简介这套中国地貌栅格数据以全国地貌分区为对象按照海拔与起伏度组合划分出低海拔平原、中海拔台地、小起伏中山、极大起伏极高山等26个类型适合GIS、地理与遥感领域的科研人员、学生及规划从业者用于地貌制图、空间分析与教学演示。压缩包共7个文件核心为tif栅格图层配套ovr金字塔、dbf属性表、tfw坐标文件、xml元数据以及docx代码表说明整体仅3.92MB轻量易用。目前已有1231人学习/下载。数据分类代码体系清晰如11代表低海拔平原、41代表小起伏低山可直接在ArcGIS等主流GIS软件中挂接属性表进行渲染与查询配合说明文档可快速理解各地貌类型含义适用于区域地貌对比、海拔与起伏度分析、专题图制作等场景也可作为教学案例数据使用帮助用户快速提取区域地貌特征。1. 别把全国地貌栅格当成普通图片它是一张带分类属性的科学数据集第一次拿到“中国地貌栅格数据.rar”的人通常会有两种反应要么双击“中国地貌.tif”却发现颜色灰蒙蒙要么用ArcGIS加载后以为数据损坏。其实这套数据的内核不是影像而是一张全国范围的地貌类型分类图每个像元存的是一个编码例如11代表低海拔平原73代表大起伏高山。这类栅格数据在土壤侵蚀、生态区划、工程选线等领域是标准底图适合有GIS基础但需要快速应用全国地貌分区的从业者。和普通遥感影像最大的区别在于它的属性信息写在旁边的代码表和VAT文件里读图时必须先理解编码体系否则后面的统计、重分类、叠加都无从下手。2. 理解中国地貌栅格数据的分类体系与文件结构2.1 从代码表说明.docx到地貌类型编码栅格值背后的分级逻辑压缩包里的“代码表说明.docx”是所有问题的起点。这套数据的编码方式是十位数字表示海拔带个位数字表示形态起伏。1到4分别为低海拔、中海拔、高海拔、极高海拔个位1到4分别对应平原、台地、丘陵、山地4到7再按起伏度分出小起伏、中起伏、大起伏、极大起伏。也就是说13是高海拔平原51是小起伏低山74是极大起伏极高山。注意原表中没有61、71因为高海拔区很少出现低起伏的低山这种空缺本身就是地貌学规律的体现。我建议在看数据之前先把这个映射表单独保存成一份CSV比如每行是“value, name, group”。这样后续用Python分析时直接读CSV就能把像元值翻译成中文省得每写一个脚本都要复制一遍dict。如果你处理过“我国地下水位栅格数据shp”或“全国城市形态栅格数据集”这类似的产品会发现它们普遍也是这种“编码外部说明”的组织方式只不过attribute表里字段名不同。下表只列前几类示意完整31类以docx为准栅格值地貌名称大类11低海拔平原平原12中海拔平原平原21低海拔台地台地31低海拔丘陵丘陵42小起伏中山山地64大起伏极高山山地表格的价值在于让你快速理解VAT里的一个数字并不只是一个ID它背后代表的是海拔带与起伏度的交叉组合。实际工作中有人会直接拿“中国地貌.tif.vat.dbf”里的Value字段和Count字段做面积估算但我更推荐先用gdalinfo检查一下这个VAT是否与tif完全同步老版本ArcGIS生成的缓存数据偶尔会因为后续裁剪操作变得不一致。2.2 tif / tfw / ovr / vat每个文件在GIS里负责干什么解压后的文件组成很有代表性主文件“中国地貌.tif”保存像元值“中国地貌.tfw”是一个六行文本的世界文件记录了左上角坐标、像元大小和旋转参数“中国地貌.tif.ovr”是金字塔文件用于缩放显示时加速“中国地貌.tif.vat.dbf”和“.vat.cpg”是栅格属性表dbf里存放Value、Count等字段cpg声明字符集“.aux.xml”则是ArcGIS自动生成的辅助元数据。很多新手会把这些小文件当成垃圾文件删掉这是错误做法。如果删掉tfwQGIS还能从tif头部读到地理参考但部分命令行工具在tif内部坐标缺失时会无法定位如果删掉ovr大尺度缩放会非常慢如果删掉vat分类图例显示就会失去颜色映射。最稳妥的方式是整套目录原样保留分析时只对主文件操作不做任何删除。注意.ovr 和 .tif 必须保持同名同目录一旦拆分GDAL 会在下次读取时尝试重建金字塔反而会拖慢你的工作流。2.3 为什么用栅格而非矢量表达地貌分区从制图学角度看地貌分区在自然界里没有刚性边界栅格比矢量更适合表达这种连续渐变的地理现象。矢量多边形在处理“低海拔向高海拔过渡”时会产生大量碎斑和锯齿而固定分辨率栅格每个像元独立分类后续与DEM、土壤类型、降雨量等栅格做像元级叠加时不需要做复杂的拓扑预处理直接矩阵运算即可。理论上讲这套栅格数据通常采用约1km的网格分辨率对于全国尺度的区域分析已经很合适但如果你研究县域或更小区域1km像元会明显抹平微地貌这时就需要寻找更精细的数据或用矢量原始资料重新采样。在选用数据前用tfw文件里的像元大小估算一下研究区内横向像元数是一个很实用的前置判断方法。例如tfw第1行是0.0083333说明像元边长约1/120度在低纬度地面距离约920米在高纬度则更短这种非等面积的特性决定了面积统计前必须先投影。3. 解压、检查与正确打开中国地貌.tif3.1 rar解压与完整性校验拿到“中国地貌栅格数据.rar”后不要直接在解压工具里双击预览tif因为RAR内的tif可能是不连续存储的预览时万一文件损坏不容易发现。我习惯先在命令行里做一次测试7z t 中国地貌栅格数据.rar 7z x 中国地貌栅格数据.rar -o./china_geomorph7z t是测试命令会逐个文件计算CRC输出“Everything is Ok”才能确定压缩包完好7z x是解压到指定目录注意-o和目录之间没有空格。如果你只有unrar也可以使用unrar t 中国地貌栅格数据.rar。测试中一旦报CRC错误说明RAR文件下载不完整或磁盘有坏道重新比对下载文件的哈希值再解压。解压完成后最好再用file命令确认tif格式file 中国地貌.tif输出若为“TIFF image data, little-endian”则正常。如果输出“data”或“HTML document”多半是解压出来的文件不对需要重新检查RAR。这行代码虽然简单但能帮你避开“数据已解压却打不开”的尴尬。3.2 用gdalinfo快速掌握栅格元数据接下来用GDAL读取元数据这比在ArcGIS界面上看更快gdalinfo 中国地貌.tif重点关注六项Size行列数、Origin左上角地理坐标、Pixel Size像元尺寸、Coordinate SystemWKT投影信息、Band 1 Minimum/Maximum像元值范围、Overviews金字塔层级。例如Origin为(73.5, 53.5)Pixel Size为(0.0083333, -0.0083333)就说明数据是经纬度网格地面像元大小约1km。如果Band 1的Maximum是74则说明这31类编码都被完整包含。如果返回值中投影为空或Origin是(0,0)说明数据丢失了地理参考需要从tfw文件手工恢复。恢复方式是读取tfw的六行数字通过gdal_translate -a_ullr命令行写入范围。例如gdal_translate -a_ullr 73.5 53.5 135.0 18.0 -a_srs EPSG:4326 中国地貌.tif 中国地貌_fixed.tif参数说明-a_ullr后面依次是左上、右下的经纬度边界-a_srs强制指定WGS84坐标系。这个方法适合tfw文件没被修改的情况但前提是你从tfw里读到的坐标是正确的。相比直接在QGIS里手动改CRS这样生成的GeoTIFF在后续命令行处理时更不容易出错。3.3 打开后黑屏、错位和NoData的正确处理在QGIS中打开如果是全黑或花屏原因是分类栅格没有默认颜色映射。此时打开图层属性-符号化选择“Paletted/Unique values”点击“Classify”按钮让QGIS读取所有唯一值再手动替换颜色。如果值很多也可以加载.qml样式文件但这里不展开。另一个高频问题是大范围全黑但局部有颜色遇到这种可以考虑是NoData值没设对。分类栅格常把海洋或无效区域设为0而最低值显示为0时也会被拉伸成黑色。用gdal_translate将0显式定义为NoDatagdal_translate -a_nodata 0 -ot Byte -co COMPRESSLZW 中国地貌.tif 中国地貌_nodata.tif-a_nodata 0设置无效值为0-ot Byte保持8位无符号整型-co COMPRESSLZW无损压缩。之后再加载黑色区域会自动变成透明或背景色后续统计也更方便。注意有些早期版本数据用了255作为NoData这时可以用gdalinfo查看原始tif的NoData参数不要直接套0。4. 实际应用地貌栅格的面积统计、重分类与制图4.1 用numpyGDAL统计各地貌类型面积拿到有效的分类栅格后最常见的需求是“低海拔平原占多少大起伏高山占多少”。使用ArcGIS的Spatial Analyst课可以完成但我更习惯用Python直接跑from osgeo import gdal import numpy as np ds gdal.Open(中国地貌_nodata.tif) band ds.GetRasterBand(1) data band.ReadAsArray() gt ds.GetGeoTransform() pixel_area abs(gt[1] * gt[5]) # 像元宽*高平方度 valid data[data 0] vals, counts np.unique(valid, return_countsTrue) for v, c in zip(vals, counts): print(v, c, round(c * pixel_area, 4))代码逻辑说明ReadAsArray把整个全国栅格读入内存约4200万整数内存占用约40MB普通电脑可以接受np.unique统计每个ID的像元数像元面积由GeoTransform里的宽度gt[1]和高度绝对值gt[5]相乘得到。注意这里是“平方度”若要公制面积先投影到等积投影如Albers再统计。输出结果后结合代码表docx就能得到类似“低海拔平原面积XX平方度”的结论。这种统计方式比在ArcGIS里打开属性表逐项导出更快而且你能看清楚哪些值出现了0次这对数据质量检查也有帮助。比如某个应该在陆地上的代码出现0次很可能说明原始分类在某个区域断档。4.2 31类合并成平原、台地、丘陵、山地四大类实际报告中通常不需要细分31类而只需要“平原/台地/丘陵/山地”。由于编码规则明确你完全可以在数组上向量化重分类而不需要使用重分类工具data data.astype(np.int16) last data % 10 res np.zeros_like(data) res[(last 1) (data 0)] 1 # 平原 res[(last 2) (data 0)] 2 # 台地 res[(last 3) (data 0)] 3 # 丘陵 res[(last 4) (data 0)] 4 # 山地 res[data 0] 0逻辑说明取个位数字last判断形态类别十位和个位组合其实没有破坏这个规则因为所有平原、台地、丘陵的个位都是1、2、3山地是4~7。(data 0)条件将NoData排除最后一行的res[data 0] 0是防御性写法确保输出值域是0~4。输出到新栅格时注意使用原GeoTransform和Projection并使用COMPRESSLZWout gdal.GetDriverByName(GTiff).Create(reclass4.tif, ds.RasterXSize, ds.RasterYSize, 1, gdal.GDT_Byte, options[COMPRESSLZW]) out.SetGeoTransform(gt) out.SetProjection(ds.GetProjection()) out.GetRasterBand(1).WriteArray(res) out.FlushCache()用这个流程替代ArcGIS的Reclassify好处是你可以随时修改合并规则。比如把“台地”和“平原”合并为“平地”只要改变last 2那一行右侧的赋值即可不用重新打开对话框。4.3 渲染成图面向汇报的地貌专题图如果你只是截图看效果可以用QGIS的“Paletted/Unique values”手动配色但如果要输出矢量出版级地图我更愿意用GDAL生成一个海拔带栅格再用伪彩色渲染。先通过gdal_calc.py提取海拔带gdal_calc.py -A 中国地貌_nodata.tif --outfilerelief.tif --calcfloor(A/10) --NoDataValue0calc中floor(A/10)将每个像元值整除10得到1~4的海拔带NoDataValue保留0。随后在QGIS中对relief.tif选择“Singleband pseudocolor”搭配一套由绿到棕的Type值渐变方案比如低海拔用绿、中海拔用黄、高海拔用橙、极高海拔用红棕图例按1~4标注为“低海拔/中海拔/高海拔/极高海拔”输出300dpiPNG即可。这样处理的好处是地貌形态平原/山地通过专题符号化表达海拔带通过颜色表达一张图既有生态含义又有地形特征。5. 进阶技巧给中国地貌栅格补上金字塔并提取像素值5.1 分类栅格的金字塔重建别让平均重采样毁掉分类值用QGIS打开全国tif时出现明显拖动卡顿通常是因为.ovr金字塔缺失或者与tif文件分离。如果数据目录里有同名.ovr请直接保留如果没有手动用gdaladdo生成gdaladdo -r nearest 中国地貌_nodata.tif 2 4 8 16参数说明-r nearest是关键分类栅格重采样只能用邻近法。若误用average缩放时像元值会被平均比如11和12变成11.5图上会出现一些不存在的类后续如果对金字塔采样而不是原始数据统计结果会受污染。2、4、8、16构建四级金字塔覆盖从轻到重的缩放级别。生成后检查tif目录是否出现.ovr文件即可。5.2 用逆仿射变换批量提取点上的地貌代码当你在野外采集了一批调查点或者在矢量上有乡镇点数据需要给每个点附加地貌类型时不必用到昂贵的空间分析扩展。直接用GDAL的仿射变换反算就可以from osgeo import gdal ds gdal.Open(中国地貌_nodata.tif) gt ds.GetGeoTransform() inv_gt gdal.InvGeoTransform(gt) def value_at(lon, lat): px, py [int(round(v)) for v in gdal.ApplyGeoTransform(inv_gt, lon, lat)] if px 0 or py 0 or px ds.RasterXSize or py ds.RasterYSize: return -1 return ds.GetRasterBand(1).ReadAsArray(px, py, 1, 1)[0][0] print(value_at(104.5, 32.5))代码说明InvGeoTransform将地图坐标转为像元坐标ApplyGeoTransform按六参数模型计算行列ReadAsArray(px, py, 1, 1)只读一个像元适合几千个点的批处理。注意越界点要先过滤否则返回错误。如果你的点文件是投影坐标那么lon, lat换成投影坐标同样适用只要坐标系与tif统一即可。每个点都对应一个4位编码拿去做后续的Logistic回归或样方统计正好与植被、土壤数据对齐。如果提取结果大面积返回0请回到第3.3节检查NoData设置不要急着把0当作“水面”写进报告。本文还有配套的精品资源点击获取 中国地貌栅格数据全解析：编码体系、文件结构与Python处理
简介这套中国地貌栅格数据以全国地貌分区为对象按照海拔与起伏度组合划分出低海拔平原、中海拔台地、小起伏中山、极大起伏极高山等26个类型适合GIS、地理与遥感领域的科研人员、学生及规划从业者用于地貌制图、空间分析与教学演示。压缩包共7个文件核心为tif栅格图层配套ovr金字塔、dbf属性表、tfw坐标文件、xml元数据以及docx代码表说明整体仅3.92MB轻量易用。目前已有1231人学习/下载。数据分类代码体系清晰如11代表低海拔平原、41代表小起伏低山可直接在ArcGIS等主流GIS软件中挂接属性表进行渲染与查询配合说明文档可快速理解各地貌类型含义适用于区域地貌对比、海拔与起伏度分析、专题图制作等场景也可作为教学案例数据使用帮助用户快速提取区域地貌特征。1. 别把全国地貌栅格当成普通图片它是一张带分类属性的科学数据集第一次拿到“中国地貌栅格数据.rar”的人通常会有两种反应要么双击“中国地貌.tif”却发现颜色灰蒙蒙要么用ArcGIS加载后以为数据损坏。其实这套数据的内核不是影像而是一张全国范围的地貌类型分类图每个像元存的是一个编码例如11代表低海拔平原73代表大起伏高山。这类栅格数据在土壤侵蚀、生态区划、工程选线等领域是标准底图适合有GIS基础但需要快速应用全国地貌分区的从业者。和普通遥感影像最大的区别在于它的属性信息写在旁边的代码表和VAT文件里读图时必须先理解编码体系否则后面的统计、重分类、叠加都无从下手。2. 理解中国地貌栅格数据的分类体系与文件结构2.1 从代码表说明.docx到地貌类型编码栅格值背后的分级逻辑压缩包里的“代码表说明.docx”是所有问题的起点。这套数据的编码方式是十位数字表示海拔带个位数字表示形态起伏。1到4分别为低海拔、中海拔、高海拔、极高海拔个位1到4分别对应平原、台地、丘陵、山地4到7再按起伏度分出小起伏、中起伏、大起伏、极大起伏。也就是说13是高海拔平原51是小起伏低山74是极大起伏极高山。注意原表中没有61、71因为高海拔区很少出现低起伏的低山这种空缺本身就是地貌学规律的体现。我建议在看数据之前先把这个映射表单独保存成一份CSV比如每行是“value, name, group”。这样后续用Python分析时直接读CSV就能把像元值翻译成中文省得每写一个脚本都要复制一遍dict。如果你处理过“我国地下水位栅格数据shp”或“全国城市形态栅格数据集”这类似的产品会发现它们普遍也是这种“编码外部说明”的组织方式只不过attribute表里字段名不同。下表只列前几类示意完整31类以docx为准栅格值地貌名称大类11低海拔平原平原12中海拔平原平原21低海拔台地台地31低海拔丘陵丘陵42小起伏中山山地64大起伏极高山山地表格的价值在于让你快速理解VAT里的一个数字并不只是一个ID它背后代表的是海拔带与起伏度的交叉组合。实际工作中有人会直接拿“中国地貌.tif.vat.dbf”里的Value字段和Count字段做面积估算但我更推荐先用gdalinfo检查一下这个VAT是否与tif完全同步老版本ArcGIS生成的缓存数据偶尔会因为后续裁剪操作变得不一致。2.2 tif / tfw / ovr / vat每个文件在GIS里负责干什么解压后的文件组成很有代表性主文件“中国地貌.tif”保存像元值“中国地貌.tfw”是一个六行文本的世界文件记录了左上角坐标、像元大小和旋转参数“中国地貌.tif.ovr”是金字塔文件用于缩放显示时加速“中国地貌.tif.vat.dbf”和“.vat.cpg”是栅格属性表dbf里存放Value、Count等字段cpg声明字符集“.aux.xml”则是ArcGIS自动生成的辅助元数据。很多新手会把这些小文件当成垃圾文件删掉这是错误做法。如果删掉tfwQGIS还能从tif头部读到地理参考但部分命令行工具在tif内部坐标缺失时会无法定位如果删掉ovr大尺度缩放会非常慢如果删掉vat分类图例显示就会失去颜色映射。最稳妥的方式是整套目录原样保留分析时只对主文件操作不做任何删除。注意.ovr 和 .tif 必须保持同名同目录一旦拆分GDAL 会在下次读取时尝试重建金字塔反而会拖慢你的工作流。2.3 为什么用栅格而非矢量表达地貌分区从制图学角度看地貌分区在自然界里没有刚性边界栅格比矢量更适合表达这种连续渐变的地理现象。矢量多边形在处理“低海拔向高海拔过渡”时会产生大量碎斑和锯齿而固定分辨率栅格每个像元独立分类后续与DEM、土壤类型、降雨量等栅格做像元级叠加时不需要做复杂的拓扑预处理直接矩阵运算即可。理论上讲这套栅格数据通常采用约1km的网格分辨率对于全国尺度的区域分析已经很合适但如果你研究县域或更小区域1km像元会明显抹平微地貌这时就需要寻找更精细的数据或用矢量原始资料重新采样。在选用数据前用tfw文件里的像元大小估算一下研究区内横向像元数是一个很实用的前置判断方法。例如tfw第1行是0.0083333说明像元边长约1/120度在低纬度地面距离约920米在高纬度则更短这种非等面积的特性决定了面积统计前必须先投影。3. 解压、检查与正确打开中国地貌.tif3.1 rar解压与完整性校验拿到“中国地貌栅格数据.rar”后不要直接在解压工具里双击预览tif因为RAR内的tif可能是不连续存储的预览时万一文件损坏不容易发现。我习惯先在命令行里做一次测试7z t 中国地貌栅格数据.rar 7z x 中国地貌栅格数据.rar -o./china_geomorph7z t是测试命令会逐个文件计算CRC输出“Everything is Ok”才能确定压缩包完好7z x是解压到指定目录注意-o和目录之间没有空格。如果你只有unrar也可以使用unrar t 中国地貌栅格数据.rar。测试中一旦报CRC错误说明RAR文件下载不完整或磁盘有坏道重新比对下载文件的哈希值再解压。解压完成后最好再用file命令确认tif格式file 中国地貌.tif输出若为“TIFF image data, little-endian”则正常。如果输出“data”或“HTML document”多半是解压出来的文件不对需要重新检查RAR。这行代码虽然简单但能帮你避开“数据已解压却打不开”的尴尬。3.2 用gdalinfo快速掌握栅格元数据接下来用GDAL读取元数据这比在ArcGIS界面上看更快gdalinfo 中国地貌.tif重点关注六项Size行列数、Origin左上角地理坐标、Pixel Size像元尺寸、Coordinate SystemWKT投影信息、Band 1 Minimum/Maximum像元值范围、Overviews金字塔层级。例如Origin为(73.5, 53.5)Pixel Size为(0.0083333, -0.0083333)就说明数据是经纬度网格地面像元大小约1km。如果Band 1的Maximum是74则说明这31类编码都被完整包含。如果返回值中投影为空或Origin是(0,0)说明数据丢失了地理参考需要从tfw文件手工恢复。恢复方式是读取tfw的六行数字通过gdal_translate -a_ullr命令行写入范围。例如gdal_translate -a_ullr 73.5 53.5 135.0 18.0 -a_srs EPSG:4326 中国地貌.tif 中国地貌_fixed.tif参数说明-a_ullr后面依次是左上、右下的经纬度边界-a_srs强制指定WGS84坐标系。这个方法适合tfw文件没被修改的情况但前提是你从tfw里读到的坐标是正确的。相比直接在QGIS里手动改CRS这样生成的GeoTIFF在后续命令行处理时更不容易出错。3.3 打开后黑屏、错位和NoData的正确处理在QGIS中打开如果是全黑或花屏原因是分类栅格没有默认颜色映射。此时打开图层属性-符号化选择“Paletted/Unique values”点击“Classify”按钮让QGIS读取所有唯一值再手动替换颜色。如果值很多也可以加载.qml样式文件但这里不展开。另一个高频问题是大范围全黑但局部有颜色遇到这种可以考虑是NoData值没设对。分类栅格常把海洋或无效区域设为0而最低值显示为0时也会被拉伸成黑色。用gdal_translate将0显式定义为NoDatagdal_translate -a_nodata 0 -ot Byte -co COMPRESSLZW 中国地貌.tif 中国地貌_nodata.tif-a_nodata 0设置无效值为0-ot Byte保持8位无符号整型-co COMPRESSLZW无损压缩。之后再加载黑色区域会自动变成透明或背景色后续统计也更方便。注意有些早期版本数据用了255作为NoData这时可以用gdalinfo查看原始tif的NoData参数不要直接套0。4. 实际应用地貌栅格的面积统计、重分类与制图4.1 用numpyGDAL统计各地貌类型面积拿到有效的分类栅格后最常见的需求是“低海拔平原占多少大起伏高山占多少”。使用ArcGIS的Spatial Analyst课可以完成但我更习惯用Python直接跑from osgeo import gdal import numpy as np ds gdal.Open(中国地貌_nodata.tif) band ds.GetRasterBand(1) data band.ReadAsArray() gt ds.GetGeoTransform() pixel_area abs(gt[1] * gt[5]) # 像元宽*高平方度 valid data[data 0] vals, counts np.unique(valid, return_countsTrue) for v, c in zip(vals, counts): print(v, c, round(c * pixel_area, 4))代码逻辑说明ReadAsArray把整个全国栅格读入内存约4200万整数内存占用约40MB普通电脑可以接受np.unique统计每个ID的像元数像元面积由GeoTransform里的宽度gt[1]和高度绝对值gt[5]相乘得到。注意这里是“平方度”若要公制面积先投影到等积投影如Albers再统计。输出结果后结合代码表docx就能得到类似“低海拔平原面积XX平方度”的结论。这种统计方式比在ArcGIS里打开属性表逐项导出更快而且你能看清楚哪些值出现了0次这对数据质量检查也有帮助。比如某个应该在陆地上的代码出现0次很可能说明原始分类在某个区域断档。4.2 31类合并成平原、台地、丘陵、山地四大类实际报告中通常不需要细分31类而只需要“平原/台地/丘陵/山地”。由于编码规则明确你完全可以在数组上向量化重分类而不需要使用重分类工具data data.astype(np.int16) last data % 10 res np.zeros_like(data) res[(last 1) (data 0)] 1 # 平原 res[(last 2) (data 0)] 2 # 台地 res[(last 3) (data 0)] 3 # 丘陵 res[(last 4) (data 0)] 4 # 山地 res[data 0] 0逻辑说明取个位数字last判断形态类别十位和个位组合其实没有破坏这个规则因为所有平原、台地、丘陵的个位都是1、2、3山地是4~7。(data 0)条件将NoData排除最后一行的res[data 0] 0是防御性写法确保输出值域是0~4。输出到新栅格时注意使用原GeoTransform和Projection并使用COMPRESSLZWout gdal.GetDriverByName(GTiff).Create(reclass4.tif, ds.RasterXSize, ds.RasterYSize, 1, gdal.GDT_Byte, options[COMPRESSLZW]) out.SetGeoTransform(gt) out.SetProjection(ds.GetProjection()) out.GetRasterBand(1).WriteArray(res) out.FlushCache()用这个流程替代ArcGIS的Reclassify好处是你可以随时修改合并规则。比如把“台地”和“平原”合并为“平地”只要改变last 2那一行右侧的赋值即可不用重新打开对话框。4.3 渲染成图面向汇报的地貌专题图如果你只是截图看效果可以用QGIS的“Paletted/Unique values”手动配色但如果要输出矢量出版级地图我更愿意用GDAL生成一个海拔带栅格再用伪彩色渲染。先通过gdal_calc.py提取海拔带gdal_calc.py -A 中国地貌_nodata.tif --outfilerelief.tif --calcfloor(A/10) --NoDataValue0calc中floor(A/10)将每个像元值整除10得到1~4的海拔带NoDataValue保留0。随后在QGIS中对relief.tif选择“Singleband pseudocolor”搭配一套由绿到棕的Type值渐变方案比如低海拔用绿、中海拔用黄、高海拔用橙、极高海拔用红棕图例按1~4标注为“低海拔/中海拔/高海拔/极高海拔”输出300dpiPNG即可。这样处理的好处是地貌形态平原/山地通过专题符号化表达海拔带通过颜色表达一张图既有生态含义又有地形特征。5. 进阶技巧给中国地貌栅格补上金字塔并提取像素值5.1 分类栅格的金字塔重建别让平均重采样毁掉分类值用QGIS打开全国tif时出现明显拖动卡顿通常是因为.ovr金字塔缺失或者与tif文件分离。如果数据目录里有同名.ovr请直接保留如果没有手动用gdaladdo生成gdaladdo -r nearest 中国地貌_nodata.tif 2 4 8 16参数说明-r nearest是关键分类栅格重采样只能用邻近法。若误用average缩放时像元值会被平均比如11和12变成11.5图上会出现一些不存在的类后续如果对金字塔采样而不是原始数据统计结果会受污染。2、4、8、16构建四级金字塔覆盖从轻到重的缩放级别。生成后检查tif目录是否出现.ovr文件即可。5.2 用逆仿射变换批量提取点上的地貌代码当你在野外采集了一批调查点或者在矢量上有乡镇点数据需要给每个点附加地貌类型时不必用到昂贵的空间分析扩展。直接用GDAL的仿射变换反算就可以from osgeo import gdal ds gdal.Open(中国地貌_nodata.tif) gt ds.GetGeoTransform() inv_gt gdal.InvGeoTransform(gt) def value_at(lon, lat): px, py [int(round(v)) for v in gdal.ApplyGeoTransform(inv_gt, lon, lat)] if px 0 or py 0 or px ds.RasterXSize or py ds.RasterYSize: return -1 return ds.GetRasterBand(1).ReadAsArray(px, py, 1, 1)[0][0] print(value_at(104.5, 32.5))代码说明InvGeoTransform将地图坐标转为像元坐标ApplyGeoTransform按六参数模型计算行列ReadAsArray(px, py, 1, 1)只读一个像元适合几千个点的批处理。注意越界点要先过滤否则返回错误。如果你的点文件是投影坐标那么lon, lat换成投影坐标同样适用只要坐标系与tif统一即可。每个点都对应一个4位编码拿去做后续的Logistic回归或样方统计正好与植被、土壤数据对齐。如果提取结果大面积返回0请回到第3.3节检查NoData设置不要急着把0当作“水面”写进报告。本文还有配套的精品资源点击获取