简介面向GIS学习与科研工作者的一份全国尺度地貌栅格数据数据覆盖完整、分类体系统一采用全国地貌分区编码覆盖低海拔平原、中海拔台地、高海拔丘陵、极大起伏极高山等多种地貌类型便于开展地貌制图、区域对比与教学实训。压缩包共7个文件主体为TIFF栅格及其金字塔文件在ArcGIS等软件中可快速显示同时附TFW坐标参考、属性表DBF、配套代码表说明等便于理解分区编码含义并做后续统计。资源包仅3.92MB轻量易得已有1231人学习下载。对需要直接使用标准地貌数据完成空间分析、专题图制作或课堂案例的用户无需再从多源数据中自行拼接、分类和配准直接加载即可使用属性字段与说明文档让编码一目了然可较快投入地形地貌制图与预处理工作。1. 打开中国地貌栅格数据.rar 之前先搞清楚里面是什么在网盘、地质资料库或科研数据共享平台上拿到中国地貌栅格数据.rar第一反应通常是双击解压然后拖进 ArcGIS 或 QGIS 里拉伸渲染。结果往往有三种图像全黑、地图位置错位、或者颜色看着像地貌图却怎么也叠不到地形上。原因很简单这类压缩包基本来自全国 1:100 万或 1:25 万地貌图的栅格化成果内部可能是 GeoTIFF、ERDAS IMG、ESRI ArcGRID甚至是一套带着投影文件和属性表的栅格工程。更关键的是像元值不是高程而是地貌类型编码。真正要用起来至少要顺着格式识别、投影统一、重分类统计、制图表达这条链路走一遍。这篇内容就是把这些步骤落地成可执行的命令和参数给做地貌分析、生态评价和国土空间规划的同行一条稳妥的实现路径。2. 解压与栅格格式识别从 .rar 到可在 Python 中读取的数组拿到压缩包先别急着双击解压。很多栅格数据集是地理信息系统工程打包出来的除主文件外还带有.tfw世界文件、.prj投影文件、.dbf属性表甚至整套 ArcGRID 目录。如果只管主文件而忽略了侧车文件后续打开时坐标信息就会丢失。2.1 用命令行把压缩包拆开先看文件清单Windows 下可以优先用 7-ZipLinux/macOS 下用 unrar 或 7z 也可以但我一般习惯先用unrar l列出压缩包里的完整条目确认是单幅 TIFF 还是一个多文件栅格再决定怎么解压。# 列出压缩包内容确认文件结构和是否有隐藏文件 unrar l China_Geomorphology.rar # 解压到一个纯英文路径避免 GDAL 某些版本对乱码路径报错 7z x China_Geomorphology.rar -o~/geo_data/unrar l不产生实际文件只是预览这对判断包内是否嵌套目录和缓存文件很有用。使用7z x时注意-o后面紧跟输出目录没有空格这和其他参数的习惯不一样。解压完成后用ls -R查看文件如果发现中文文件名在软件里无法打开可以统一重命名成拼音或英文顺便把路径中的空格也去掉。GDAL 对含中文和空格的路径支持不算稳定尤其是老数据配.tfw时容易读出空投影。文件扩展名常见格式说明.tif/.tiffGeoTIFF新生数据多为此格式坐标写入文件头.tfwWorld File单独的地理配准文件必须和 tif 同目录.imgERDAS IMAGINE常见于早期上图成果带.aux侧车文件.adfESRI ArcGRID以文件夹为单位保存不能单独拖动内部文件.bil/.hdrENVI/BIL需要配合.hdr头文件解释波段和类型这张表用来判断解压后数据是否完整。特别提醒一点如果包内是 ArcGRID 格式解压后不要只拷贝其中的.adf文件而要保持整个xxx/文件夹结构和主目录同名否则 QGIS 会直接忽略这个数据集。2.2 读元数据用 gdalinfo 弄清楚投影、像素和 NoData文件解压到位后下一步不是打开渲染而是读取栅格元数据。栅格的所有关键技术参数都在元数据里包括坐标系、像素尺寸、波段数和 NoData 值。用 GDAL 自带工具最快能看清全貌gdalinfo China_Geomorphology_1m.img重点关注输出中的这一段Driver: HFA/ERDAS Imagine Size is 29184, 20416 Coordinate System is: PROJCRS[China_Lambert_Conformal_Conic, ...] Pixel Size (1000.000000000000000, -1000.000000000000000) Metadata: STATISTICS_MINIMUM1 STATISTICS_MAXIMUM18 NoData Value127Size是行列数决定后面读数组时能开多大的内存Coordinate System是真真正正决定数据往哪摆的参数NoData Value则直接关系到面积统计是否准确。很多人在统计面积时忽略 NoData把背景格子也算进地貌类型得到的结果自然偏大几倍。如果输出的Coordinate System是GCS_WGS_84说明数据是经纬度坐标不能直接用米作为单位计算面积需要先转到投影坐标。Python 侧用 Rasterio 打开会更方便因为可以直接拿到profile和boundsimport rasterio with rasterio.open(China_Geomorphology_1m.img) as src: profile src.profile # 包含驱动、波段数、数据类型、NoData crs src.crs # 坐标系可能是 None bounds src.bounds # 左、下、右、上 data src.read(1) # 读第一个波段到内存 print(profile) print(crs) print(bounds)这里直接data src.read(1)适合全国 1:100 万地貌图因为分辨率通常为 1km 或 500m整个数组只有几百万个像素内存压力小。但如果压缩包内是 30 米分辨率的派生数据不要直接read(1)建议改用src.block_shapes做分块读取否则一台 16GB 内存的机器可能会直接内存溢出。2.3 压缩包内最容易出问题的三种情况第一类缺少.prj投影文件。这是最麻烦的GDAL 会把坐标系认成未知数据加载后没有参考位置。遇到这种情况先看压缩包内有没有说明文档写明“Albers 等积投影中央经线 105E双标准纬线 25N/47N”再用手头的投影参数手动为数据赋坐标。QGIS 里也可以用“Assign Projection”工具临时指定。第二类NoData 和真实值混在一起。老数据为了节省存储背景常常不是 NoData 而是编码 0 或其他数字。此时必须显式把背景值统一成 NoData否则后面做重分类时边缘地区会被误判为某个地貌类型。第三类存在金字塔和统计缓存文件.rrd、.aux.xml。这些文件在压缩包里和主文件不在同一目录时特别容易漏掉解压后主文件一重组软件会误认为数据损坏。规避方法是在解压后保留所有文件到原目录或者直接删除缓存文件让软件重建金字塔。3. 投影与坐标系中国地貌栅格数据叠加前必做的事很多数据叠加到在线底图上会出现几十公里偏移这往往不是数据精度问题而是坐标系不一致。中国地貌栅格数据的坐标系常见有北京 54、西安 80、CGCS2000以及各种 Lambert 圆锥投影和 Albers 等积投影不能默认全都是 WGS84。3.1 先确认投影再谈处理地貌分类数据是按面积成图的1:100 万中国地貌图国家标准采用等积性质的 Albers 或双标准纬线圆锥投影因为地貌分类和面积统计必须保证面积不变形。如果原始数据是 Albers 等积投影而你直接丢进 Web 墨卡托底图渲染面积会被拉大后续统计的各类地貌面积占比失去意义。老图数字化出的数据还可能带的是“西安 1980”或“克拉索夫斯基椭球”与 CGCS2000 之间存在几十米平移。所以第一步永远是gdalinfo确认坐标系而不是直接点击重投影。怎么判断投影信息是否损毁看三个信号文件范围不在中国海域陆域边界范围内像素尺寸不是整数GDAL 输出Coordinate System is: unknown。这三个信号有一个就要先修坐标再进入后续流程。3.2 用 gdalwarp 把栅格统一到 CGCS2000 / Albers在实际项目中通用做法是把所有参与分析的栅格统一到 CGCS2000 地理坐标系EPSG:4490或 CGCS2000 高斯克吕格投影带。如果只是想和在线底图叠加用 4490 足够如果要算面积还是建议用 Albers 等积投影。下面以 CGCS2000 地理坐标为例# -r near 表示最近邻重采样保留地貌类型编码的原始值 gdalwarp -t_srs EPSG:4490 \ -r near \ -of GTiff \ China_Geomorphology_1m.img \ China_Geo_4490.tif参数说明-r near是重采样算法这里是地貌分类栅格像素值 1、2、3 是离散编码不能做平均值或整型插值否则边界会产生 1.5 这种无意义数值。-of GTiff统一输出成 GeoTIFF比.img更适合跨软件交换。如果想输出 Albers 等积投影可以用自定义投影参数gdalwarp -t_srs projaea lat_125 lat_247 lon_0105 datumCGCS2000 \ -r near \ -overwrite \ China_Geomorphology_1m.img \ China_Geo_AEA.tiflat_1和lat_2是双标准纬线lon_0是中央经线。这套参数比较接近中国制图行业对 Albers 投影的习惯设置适合全国尺度分析。如果你的研究区在某个省份可以适当调整中央经线减少经度方向的变形。3.3 按省份或流域裁切用矢量边界做 cutline拿到全中国地貌栅格后通常只需要一个研究区的范围。常见做法是直接拿行政区划或流域边界来裁剪既能压缩后续计算量也能让统计结果对应到行政单元。用gdalwarp的-cutline一条命令就能完成gdalwarp -t_srs EPSG:4490 \ -cutline 研究区边界.shp \ -crop_to_cutline \ -dstnodata 127 \ -r near \ China_Geomorphology_1m.img \ research_area.tif-cutline接收的矢量边界最好是 WGS84 或 CGCS2000 经纬度GDAL 会在内部完成坐标系转换。-crop_to_cutline让输出范围精确等于边界的范围避免四周留出大片多余区域。-dstnodata 127是经常被忽略但极其重要的参数它把裁剪后边界以外的区域设为 NoData而不是像素值 0。如果省略有些输出格式会用 0 填边后续统计时 0 被当成一个有效类型直接污染面积比例。裁剪结束后建议再复查一遍元数据gdalinfo -stats research_area.tif | grep -E Minimum|Maximum|NoData如果 Minimum 是 0 而 NoData 是 127说明背景值没处理好需要重裁或把 0 显式替换成 127。4. 地貌栅格可视化与分类统计格式和坐标都正确以后需求通常有两类一是把图渲染得清晰可读用于报告和汇报二是算出各地貌类型的面积占比。地貌栅格的像素值是分类编码不是连续高程所以默认灰度拉伸完全看不出地形语义必须给编码赋颜色表。4.1 从一片黑到能看给地貌编码配一个颜色表用 QGIS 最快捷。右键图层属性在“符号化”里把渲染类型改成“调色板/唯一值”字段选 VALUE再逐项指定颜色。如果想批量操作可以写一份 QGIS 样式文件.qml放到和栅格同名的位置下次打开就会自动加载。手边没有现成颜色表时可以用 Python 自定义palette { 1: (255, 246, 143), # 洪积台地 2: (202, 178, 214), # 冲积平原 3: (166, 206, 227), # 湖积平原 4: (141, 211, 199), # 滨海平原 5: (179, 222, 105), # 低海拔丘陵 6: (140, 150, 190), # 中海拔山地 }这组颜色参考了 Cartographic Color Brewer 的分类别配色思路相邻地貌类型使用不同色系视觉上可读性更好。实际展示时把这张字典导出成.qml或者 CSV供 QGIS 导入即可。如果需要直接在网页端展示还可以把调色板写成 JSON配合geotiff.js前端渲染。4.2 地貌类型编码与重分类把复合码归并成大形态1:100 万地貌图的编码体系并不简单一个像素值可能同时表达形态、成因和坡度等级比如“中海拔低山”可能是602之类的复合码。为了做统计分析需要把这些细类合并成一级地貌形态。最直接的方法是用 Rasterio 读取数组后做整除归并import numpy as np import rasterio with rasterio.open(research_area.tif) as src: data src.read(1) profile src.profile nodata profile.get(nodata, 127) # 假设编码规则百位为形态大类1平原 2台地 3丘陵 4山地 reclass np.where(data nodata, nodata, np.where(data 100, 1, np.floor_divide(data, 100))) reclass reclass.astype(np.int16) profile.update(dtypenp.int16) with rasterio.open(reclass_landform.tif, w, **profile) as dst: dst.write(reclass, 1)这里np.floor_divide(data, 100)正好是整除到百位把复合码变成形态大类码。不过不同分发版本的编码规则不一致操作前必须翻一遍压缩包内的图例或说明确认百位是不是形态类别。如果原始栅格已经是平铺的一级编码直接跳过这个重分类步骤不要画蛇添足。还有一个常见误用是直接用gdal_calc.py做全图计算gdal_calc.py -A input.tif --outfilelandform.tif \ --calcwhere(A127,127,floor(A/100)) \ --NoDataValue127 --typeInt16这条命令在几百 MB 的栅格上没问题但到了 GB 级大瓦片会一次性载入多数组内存很容易爆掉。所以可复用的脚本还是要写成前文的分块或者直接使用 GDAL 命令行工具配合虚拟栅格能省很多事。4.3 栅格面积统计一个像素占多少公顷中国地貌栅格数据在 Albers 等积投影下像素尺寸通常是 1000m×1000m 或 500m×500m。统计面积时不需要先矢量化再算面直接统计分类像素数乘以单像素面积就能得到可靠结果。unique, counts np.unique(reclass[reclass ! nodata], return_countsTrue) pixel_area_m2 1000 * 1000 # 根据 gdalinfo 里的 Pixel Size 调整 for cls, count in zip(unique, counts): area_sqkm count * pixel_area_m2 / 1e6 print(f地貌类型 {cls}: {area_sqkm:.2f} km²)使用这个逻辑前务必通过gdalinfo确认 Pixel Size 是等积投影下的米单位。如果还是经纬度像素实际宽度随纬度变化不能简单乘一个固定值。正确做法是先重投影到 Albers 再统计避免用近似公式和人工修正。最后的统计结果可以输出成如下结构地貌形态大类编码像素数面积(km²)平原11234512345台地267896789丘陵345674567山地42345523455如果表格要写进项目报告或论文建议在表注中标注投影坐标系、原始分辨率以及 NoData 处理方式这些信息是描述数据质量的重要部分。5. 提取等值线与沙盘叠影中国地貌栅格数据的进阶用法前面的流程已经能覆盖“把中国地貌栅格数据用起来”的主要场景。如果压缩包内是地形起伏度或高程派生栅格还能继续做等值线提取和立体叠影。即使你拿到的是纯粹的形态分类图这些技巧也可以迁移到 DEM、全国城市形态栅格数据集、地下水位栅格数据 shp 等其它栅格数据上。5.1 从起伏度栅格提取结构化等值线如果压缩包内附带的是地形起伏度栅格可以把它当作连续表面提取等值线输出给画图或者 CAD 使用。GDAL 自带gdal_contour就能直接完成gdal_contour -b 1 -a elev -i 200 -f GPKG terrain.tif contours.gpkg-b 1指定波段-a elev给输出的线图层属性字段命名-i 200表示每 200 米生成一条等高线。如果栅格是高程而不是地貌起伏度间隔可以改成 500 或 1000具体看研究区高差。输出为 GPKG 而不是 SHP是因为 GPKG 没有 Shapefile 里字段名长度 10 字符的限制线要素属性更完整加载速度也更快。5.2 山体阴影叠影与像元级交叉统计让平面分类图变得更立体常用做法是将半透明山体阴影叠加在分类栅格上。一条gdaldem命令就能生成阴影gdaldem hillshade dem.tif hillshade.tif -z 3 -alt 45 -az 315-z是垂直拉伸倍数平原地区可以给 5 增强立体感山地给 1~2避免阴影过重。-alt是太阳高度角-az是方位角315° 符合多数读图习惯。在 QGIS 中把分类栅格放到底层阴影栅格放上层混合模式选“正片叠底”再把透明度调到 40%地貌图在保持类型清晰的同时会有明显的地形起伏感。除了出图层面的技巧更实用的进阶用法是像元级交叉统计。比如拿同一研究区的城市形态栅格和地貌分类栅格做频数分析前提是两者分辨率、范围和坐标系一致不一致就先重采样。然后用np.histogram2d生成交叉表import rasterio import numpy as np with rasterio.open(landform.tif) as a, rasterio.open(city_structure.tif) as b: landform a.read(1).ravel() city b.read(1).ravel() valid_mask (landform ! a.nodata) (city ! b.nodata) cross np.histogram2d( city[valid_mask], landform[valid_mask], bins[city.max() - city.min(), landform.max() - landform.min()] )[0]交叉表每一行代表城市形态类型每一列代表地貌类型可以快速回答“低海拔台地上的建设用地占了多大比例”这类条件统计问题省去反复栅格转矢量再做空间连接的时间和精力。整个流程适合所有类似氛围的栅格数据包括地貌、土壤、地质和土地利用类型只要像素语义是一一对应的。本文还有配套的精品资源点击获取 中国地貌栅格数据从解压到分类统计完整指南
简介面向GIS学习与科研工作者的一份全国尺度地貌栅格数据数据覆盖完整、分类体系统一采用全国地貌分区编码覆盖低海拔平原、中海拔台地、高海拔丘陵、极大起伏极高山等多种地貌类型便于开展地貌制图、区域对比与教学实训。压缩包共7个文件主体为TIFF栅格及其金字塔文件在ArcGIS等软件中可快速显示同时附TFW坐标参考、属性表DBF、配套代码表说明等便于理解分区编码含义并做后续统计。资源包仅3.92MB轻量易得已有1231人学习下载。对需要直接使用标准地貌数据完成空间分析、专题图制作或课堂案例的用户无需再从多源数据中自行拼接、分类和配准直接加载即可使用属性字段与说明文档让编码一目了然可较快投入地形地貌制图与预处理工作。1. 打开中国地貌栅格数据.rar 之前先搞清楚里面是什么在网盘、地质资料库或科研数据共享平台上拿到中国地貌栅格数据.rar第一反应通常是双击解压然后拖进 ArcGIS 或 QGIS 里拉伸渲染。结果往往有三种图像全黑、地图位置错位、或者颜色看着像地貌图却怎么也叠不到地形上。原因很简单这类压缩包基本来自全国 1:100 万或 1:25 万地貌图的栅格化成果内部可能是 GeoTIFF、ERDAS IMG、ESRI ArcGRID甚至是一套带着投影文件和属性表的栅格工程。更关键的是像元值不是高程而是地貌类型编码。真正要用起来至少要顺着格式识别、投影统一、重分类统计、制图表达这条链路走一遍。这篇内容就是把这些步骤落地成可执行的命令和参数给做地貌分析、生态评价和国土空间规划的同行一条稳妥的实现路径。2. 解压与栅格格式识别从 .rar 到可在 Python 中读取的数组拿到压缩包先别急着双击解压。很多栅格数据集是地理信息系统工程打包出来的除主文件外还带有.tfw世界文件、.prj投影文件、.dbf属性表甚至整套 ArcGRID 目录。如果只管主文件而忽略了侧车文件后续打开时坐标信息就会丢失。2.1 用命令行把压缩包拆开先看文件清单Windows 下可以优先用 7-ZipLinux/macOS 下用 unrar 或 7z 也可以但我一般习惯先用unrar l列出压缩包里的完整条目确认是单幅 TIFF 还是一个多文件栅格再决定怎么解压。# 列出压缩包内容确认文件结构和是否有隐藏文件 unrar l China_Geomorphology.rar # 解压到一个纯英文路径避免 GDAL 某些版本对乱码路径报错 7z x China_Geomorphology.rar -o~/geo_data/unrar l不产生实际文件只是预览这对判断包内是否嵌套目录和缓存文件很有用。使用7z x时注意-o后面紧跟输出目录没有空格这和其他参数的习惯不一样。解压完成后用ls -R查看文件如果发现中文文件名在软件里无法打开可以统一重命名成拼音或英文顺便把路径中的空格也去掉。GDAL 对含中文和空格的路径支持不算稳定尤其是老数据配.tfw时容易读出空投影。文件扩展名常见格式说明.tif/.tiffGeoTIFF新生数据多为此格式坐标写入文件头.tfwWorld File单独的地理配准文件必须和 tif 同目录.imgERDAS IMAGINE常见于早期上图成果带.aux侧车文件.adfESRI ArcGRID以文件夹为单位保存不能单独拖动内部文件.bil/.hdrENVI/BIL需要配合.hdr头文件解释波段和类型这张表用来判断解压后数据是否完整。特别提醒一点如果包内是 ArcGRID 格式解压后不要只拷贝其中的.adf文件而要保持整个xxx/文件夹结构和主目录同名否则 QGIS 会直接忽略这个数据集。2.2 读元数据用 gdalinfo 弄清楚投影、像素和 NoData文件解压到位后下一步不是打开渲染而是读取栅格元数据。栅格的所有关键技术参数都在元数据里包括坐标系、像素尺寸、波段数和 NoData 值。用 GDAL 自带工具最快能看清全貌gdalinfo China_Geomorphology_1m.img重点关注输出中的这一段Driver: HFA/ERDAS Imagine Size is 29184, 20416 Coordinate System is: PROJCRS[China_Lambert_Conformal_Conic, ...] Pixel Size (1000.000000000000000, -1000.000000000000000) Metadata: STATISTICS_MINIMUM1 STATISTICS_MAXIMUM18 NoData Value127Size是行列数决定后面读数组时能开多大的内存Coordinate System是真真正正决定数据往哪摆的参数NoData Value则直接关系到面积统计是否准确。很多人在统计面积时忽略 NoData把背景格子也算进地貌类型得到的结果自然偏大几倍。如果输出的Coordinate System是GCS_WGS_84说明数据是经纬度坐标不能直接用米作为单位计算面积需要先转到投影坐标。Python 侧用 Rasterio 打开会更方便因为可以直接拿到profile和boundsimport rasterio with rasterio.open(China_Geomorphology_1m.img) as src: profile src.profile # 包含驱动、波段数、数据类型、NoData crs src.crs # 坐标系可能是 None bounds src.bounds # 左、下、右、上 data src.read(1) # 读第一个波段到内存 print(profile) print(crs) print(bounds)这里直接data src.read(1)适合全国 1:100 万地貌图因为分辨率通常为 1km 或 500m整个数组只有几百万个像素内存压力小。但如果压缩包内是 30 米分辨率的派生数据不要直接read(1)建议改用src.block_shapes做分块读取否则一台 16GB 内存的机器可能会直接内存溢出。2.3 压缩包内最容易出问题的三种情况第一类缺少.prj投影文件。这是最麻烦的GDAL 会把坐标系认成未知数据加载后没有参考位置。遇到这种情况先看压缩包内有没有说明文档写明“Albers 等积投影中央经线 105E双标准纬线 25N/47N”再用手头的投影参数手动为数据赋坐标。QGIS 里也可以用“Assign Projection”工具临时指定。第二类NoData 和真实值混在一起。老数据为了节省存储背景常常不是 NoData 而是编码 0 或其他数字。此时必须显式把背景值统一成 NoData否则后面做重分类时边缘地区会被误判为某个地貌类型。第三类存在金字塔和统计缓存文件.rrd、.aux.xml。这些文件在压缩包里和主文件不在同一目录时特别容易漏掉解压后主文件一重组软件会误认为数据损坏。规避方法是在解压后保留所有文件到原目录或者直接删除缓存文件让软件重建金字塔。3. 投影与坐标系中国地貌栅格数据叠加前必做的事很多数据叠加到在线底图上会出现几十公里偏移这往往不是数据精度问题而是坐标系不一致。中国地貌栅格数据的坐标系常见有北京 54、西安 80、CGCS2000以及各种 Lambert 圆锥投影和 Albers 等积投影不能默认全都是 WGS84。3.1 先确认投影再谈处理地貌分类数据是按面积成图的1:100 万中国地貌图国家标准采用等积性质的 Albers 或双标准纬线圆锥投影因为地貌分类和面积统计必须保证面积不变形。如果原始数据是 Albers 等积投影而你直接丢进 Web 墨卡托底图渲染面积会被拉大后续统计的各类地貌面积占比失去意义。老图数字化出的数据还可能带的是“西安 1980”或“克拉索夫斯基椭球”与 CGCS2000 之间存在几十米平移。所以第一步永远是gdalinfo确认坐标系而不是直接点击重投影。怎么判断投影信息是否损毁看三个信号文件范围不在中国海域陆域边界范围内像素尺寸不是整数GDAL 输出Coordinate System is: unknown。这三个信号有一个就要先修坐标再进入后续流程。3.2 用 gdalwarp 把栅格统一到 CGCS2000 / Albers在实际项目中通用做法是把所有参与分析的栅格统一到 CGCS2000 地理坐标系EPSG:4490或 CGCS2000 高斯克吕格投影带。如果只是想和在线底图叠加用 4490 足够如果要算面积还是建议用 Albers 等积投影。下面以 CGCS2000 地理坐标为例# -r near 表示最近邻重采样保留地貌类型编码的原始值 gdalwarp -t_srs EPSG:4490 \ -r near \ -of GTiff \ China_Geomorphology_1m.img \ China_Geo_4490.tif参数说明-r near是重采样算法这里是地貌分类栅格像素值 1、2、3 是离散编码不能做平均值或整型插值否则边界会产生 1.5 这种无意义数值。-of GTiff统一输出成 GeoTIFF比.img更适合跨软件交换。如果想输出 Albers 等积投影可以用自定义投影参数gdalwarp -t_srs projaea lat_125 lat_247 lon_0105 datumCGCS2000 \ -r near \ -overwrite \ China_Geomorphology_1m.img \ China_Geo_AEA.tiflat_1和lat_2是双标准纬线lon_0是中央经线。这套参数比较接近中国制图行业对 Albers 投影的习惯设置适合全国尺度分析。如果你的研究区在某个省份可以适当调整中央经线减少经度方向的变形。3.3 按省份或流域裁切用矢量边界做 cutline拿到全中国地貌栅格后通常只需要一个研究区的范围。常见做法是直接拿行政区划或流域边界来裁剪既能压缩后续计算量也能让统计结果对应到行政单元。用gdalwarp的-cutline一条命令就能完成gdalwarp -t_srs EPSG:4490 \ -cutline 研究区边界.shp \ -crop_to_cutline \ -dstnodata 127 \ -r near \ China_Geomorphology_1m.img \ research_area.tif-cutline接收的矢量边界最好是 WGS84 或 CGCS2000 经纬度GDAL 会在内部完成坐标系转换。-crop_to_cutline让输出范围精确等于边界的范围避免四周留出大片多余区域。-dstnodata 127是经常被忽略但极其重要的参数它把裁剪后边界以外的区域设为 NoData而不是像素值 0。如果省略有些输出格式会用 0 填边后续统计时 0 被当成一个有效类型直接污染面积比例。裁剪结束后建议再复查一遍元数据gdalinfo -stats research_area.tif | grep -E Minimum|Maximum|NoData如果 Minimum 是 0 而 NoData 是 127说明背景值没处理好需要重裁或把 0 显式替换成 127。4. 地貌栅格可视化与分类统计格式和坐标都正确以后需求通常有两类一是把图渲染得清晰可读用于报告和汇报二是算出各地貌类型的面积占比。地貌栅格的像素值是分类编码不是连续高程所以默认灰度拉伸完全看不出地形语义必须给编码赋颜色表。4.1 从一片黑到能看给地貌编码配一个颜色表用 QGIS 最快捷。右键图层属性在“符号化”里把渲染类型改成“调色板/唯一值”字段选 VALUE再逐项指定颜色。如果想批量操作可以写一份 QGIS 样式文件.qml放到和栅格同名的位置下次打开就会自动加载。手边没有现成颜色表时可以用 Python 自定义palette { 1: (255, 246, 143), # 洪积台地 2: (202, 178, 214), # 冲积平原 3: (166, 206, 227), # 湖积平原 4: (141, 211, 199), # 滨海平原 5: (179, 222, 105), # 低海拔丘陵 6: (140, 150, 190), # 中海拔山地 }这组颜色参考了 Cartographic Color Brewer 的分类别配色思路相邻地貌类型使用不同色系视觉上可读性更好。实际展示时把这张字典导出成.qml或者 CSV供 QGIS 导入即可。如果需要直接在网页端展示还可以把调色板写成 JSON配合geotiff.js前端渲染。4.2 地貌类型编码与重分类把复合码归并成大形态1:100 万地貌图的编码体系并不简单一个像素值可能同时表达形态、成因和坡度等级比如“中海拔低山”可能是602之类的复合码。为了做统计分析需要把这些细类合并成一级地貌形态。最直接的方法是用 Rasterio 读取数组后做整除归并import numpy as np import rasterio with rasterio.open(research_area.tif) as src: data src.read(1) profile src.profile nodata profile.get(nodata, 127) # 假设编码规则百位为形态大类1平原 2台地 3丘陵 4山地 reclass np.where(data nodata, nodata, np.where(data 100, 1, np.floor_divide(data, 100))) reclass reclass.astype(np.int16) profile.update(dtypenp.int16) with rasterio.open(reclass_landform.tif, w, **profile) as dst: dst.write(reclass, 1)这里np.floor_divide(data, 100)正好是整除到百位把复合码变成形态大类码。不过不同分发版本的编码规则不一致操作前必须翻一遍压缩包内的图例或说明确认百位是不是形态类别。如果原始栅格已经是平铺的一级编码直接跳过这个重分类步骤不要画蛇添足。还有一个常见误用是直接用gdal_calc.py做全图计算gdal_calc.py -A input.tif --outfilelandform.tif \ --calcwhere(A127,127,floor(A/100)) \ --NoDataValue127 --typeInt16这条命令在几百 MB 的栅格上没问题但到了 GB 级大瓦片会一次性载入多数组内存很容易爆掉。所以可复用的脚本还是要写成前文的分块或者直接使用 GDAL 命令行工具配合虚拟栅格能省很多事。4.3 栅格面积统计一个像素占多少公顷中国地貌栅格数据在 Albers 等积投影下像素尺寸通常是 1000m×1000m 或 500m×500m。统计面积时不需要先矢量化再算面直接统计分类像素数乘以单像素面积就能得到可靠结果。unique, counts np.unique(reclass[reclass ! nodata], return_countsTrue) pixel_area_m2 1000 * 1000 # 根据 gdalinfo 里的 Pixel Size 调整 for cls, count in zip(unique, counts): area_sqkm count * pixel_area_m2 / 1e6 print(f地貌类型 {cls}: {area_sqkm:.2f} km²)使用这个逻辑前务必通过gdalinfo确认 Pixel Size 是等积投影下的米单位。如果还是经纬度像素实际宽度随纬度变化不能简单乘一个固定值。正确做法是先重投影到 Albers 再统计避免用近似公式和人工修正。最后的统计结果可以输出成如下结构地貌形态大类编码像素数面积(km²)平原11234512345台地267896789丘陵345674567山地42345523455如果表格要写进项目报告或论文建议在表注中标注投影坐标系、原始分辨率以及 NoData 处理方式这些信息是描述数据质量的重要部分。5. 提取等值线与沙盘叠影中国地貌栅格数据的进阶用法前面的流程已经能覆盖“把中国地貌栅格数据用起来”的主要场景。如果压缩包内是地形起伏度或高程派生栅格还能继续做等值线提取和立体叠影。即使你拿到的是纯粹的形态分类图这些技巧也可以迁移到 DEM、全国城市形态栅格数据集、地下水位栅格数据 shp 等其它栅格数据上。5.1 从起伏度栅格提取结构化等值线如果压缩包内附带的是地形起伏度栅格可以把它当作连续表面提取等值线输出给画图或者 CAD 使用。GDAL 自带gdal_contour就能直接完成gdal_contour -b 1 -a elev -i 200 -f GPKG terrain.tif contours.gpkg-b 1指定波段-a elev给输出的线图层属性字段命名-i 200表示每 200 米生成一条等高线。如果栅格是高程而不是地貌起伏度间隔可以改成 500 或 1000具体看研究区高差。输出为 GPKG 而不是 SHP是因为 GPKG 没有 Shapefile 里字段名长度 10 字符的限制线要素属性更完整加载速度也更快。5.2 山体阴影叠影与像元级交叉统计让平面分类图变得更立体常用做法是将半透明山体阴影叠加在分类栅格上。一条gdaldem命令就能生成阴影gdaldem hillshade dem.tif hillshade.tif -z 3 -alt 45 -az 315-z是垂直拉伸倍数平原地区可以给 5 增强立体感山地给 1~2避免阴影过重。-alt是太阳高度角-az是方位角315° 符合多数读图习惯。在 QGIS 中把分类栅格放到底层阴影栅格放上层混合模式选“正片叠底”再把透明度调到 40%地貌图在保持类型清晰的同时会有明显的地形起伏感。除了出图层面的技巧更实用的进阶用法是像元级交叉统计。比如拿同一研究区的城市形态栅格和地貌分类栅格做频数分析前提是两者分辨率、范围和坐标系一致不一致就先重采样。然后用np.histogram2d生成交叉表import rasterio import numpy as np with rasterio.open(landform.tif) as a, rasterio.open(city_structure.tif) as b: landform a.read(1).ravel() city b.read(1).ravel() valid_mask (landform ! a.nodata) (city ! b.nodata) cross np.histogram2d( city[valid_mask], landform[valid_mask], bins[city.max() - city.min(), landform.max() - landform.min()] )[0]交叉表每一行代表城市形态类型每一列代表地貌类型可以快速回答“低海拔台地上的建设用地占了多大比例”这类条件统计问题省去反复栅格转矢量再做空间连接的时间和精力。整个流程适合所有类似氛围的栅格数据包括地貌、土壤、地质和土地利用类型只要像素语义是一一对应的。本文还有配套的精品资源点击获取