尧图网页设计 
简介一份覆盖中国全境的DEM栅格数据集分辨率为1km源自SRTM 30米分辨率原始数据经ArcGIS镶嵌拼接与重采样制作完成数据完整、边界准确可作为多种空间分析的统一底图。面向地理信息科学、测绘工程等专业学生以及需要全国尺度地形底图的科研人员可用于课程设计、地形因子提取、宏观区域分析和多源数据叠加等场景。数据采用WGS-84坐标系TIFF格式无需转换即可在ArcGIS中直接打开方便坡度坡向提取、可视域分析、栅格计算及制图出图等常见操作。压缩包共7个文件、26.51MB核心为ChinaDEM.tif栅格主数据同时收录tfw坐标配准文件、ovr金字塔文件、xml元数据及dbf属性表等辅助文件便于快速加载和二次处理。目前已有693人学习下载既能支撑教学演示也可作为科研前期数据准备数据精度经过重采样优化整体实用性强。1. 中国范围的DEM栅格数据从数据源选型到地形分析落地做全国性或区域性地形分析时第一道坎不是算法而是数据本身。中国范围的DEM栅格数据看似随处可下但不同来源的精度、坐标系、覆盖范围差异极大——SRTM在新疆某些区域存在空洞ASTER GDEM在国内高山区域噪声明显ALOS 12.5米数据精度高却动辄几十GB。更麻烦的是从地理空间数据云下载的DEM通常是分幅存储的坐标参考系为WGS84而国内项目往往要求CGCS2000或高斯投影直接拿来用轻则投影变形重则分析结果完全偏离。这篇文章就想把“获取中国范围DEM栅格数据到完成地形指标计算”这条链路拆开来讲覆盖数据源横向对比、预处理参数、地形指标计算和大范围数据的性能瓶颈给出一套可以直接照做的方案。适合的读者群很明确GIS开发、遥感应用工程师、水文或规划领域的分析人员以及做三维可视化或DSM生成DEM相关工作的朋友。无论你只需要某个省份的裁剪DEM数据还是想拼出全国一张图下面的思路都适用。2. DEM数据源选型与中国区域覆盖精度、坐标系与空洞问题2.1 三大主流DEM数据源在中国区域的真实表现覆盖中国范围的全球DEM数据源中SRTM、ASTER GDEM和ALOS AW3D30是实际项目里最常被用到的。SRTM由NASA和NGA联合获取在中国区域主要覆盖北纬60°以南国内绝大部分领土都在范围内分辨率分30米和90米两档其中90米版本在国内用得最广文件组织为分幅TIF单幅约5°×5°。ASTER GDEM由日本METI和美国NASA基于光学立体像对生成分辨率30米覆盖纬度83°N到83°S中国全域可覆盖但光学影像受云影响较大在青藏高原、天山等高山区域经常出现异常凹陷或尖峰业内俗称“麻点”噪声。ALOS AW3D30是JAXA基于PRISM立体传感器生成的全球DSM标称30米实际在绝大多数区域能体现更细致的地形纹理。而热词里常被提起的ALOS 12.5米DEM数据其实来源于ALOS PALSAR的12.5米分辨率产品精度确实高出其他全球产品一个量级但中国区域需要分块下载且原始数据为DSM——包含地表覆盖物高度不是剔除了植被和建筑物的纯地面DEM。基于这个差别做洪水淹没模拟或耕地平整度分析前需要先明确自己到底要DSM还是DEM。用DSM生成DEM的做法通常是对原始DSM做滤波剔除植被冠层和建筑物高度这在城区和密林区域效果有限但在地形起伏主导的区域表现尚可。2.2 坐标系差异与国内常用投影的选择下载到的全球DEM原始数据坐标系基本都是WGS84地理坐标系单位为度。很多初学者直接用原始数据做坡度计算或面积统计得到的结果毫无意义——坡度在经纬度坐标系下X和Y方向长度不等得到的坡度和实际物理坡度存在系统性偏差。处理中国范围的DEM栅格数据时第一步应该是确认目标坐标系。国内生产项目常用的有三类CGCS2000地理坐标系代号EPSG:4490适合全国范围展示和简单分析CGCS2000 / 3-degree Gauss-Kruger zone代号EPSG:4524到4559按经度分带适合省级或地市级项目UTM分区代号EPSG:32649到32652覆盖中国从东到西的分区如果上下游工具链是国际通用软件UTM反而更省事。需要特别注意CGCS2000高斯投影分带规则和WGS84 UTM分带规则不同高斯3度带从东经75°起每隔3°一带中央经线为75、78、81……直到135UTM则从西经180°起每隔6°一带中央经线为75、81、87……285。对于范围跨多个分带的全国数据正确做法是保留地理坐标做全国分析只在局部高精度分析时才投影到对应分带。某些省级项目直接把数据投影到CGCS2000 / Gauss-Kruger zone 20中央经线117°处理全省范围投影变形会导致东西跨度大的省份边缘区域精度损失。2.3 数据空洞和异常值的识别的快速检查法中国区域DEM一个老生常谈的痛点是SRTM在陡峭地形的数据空洞以及ASTER GDEM在湖泊、积雪区域的异常值。拿到数据后先做概览检查别直接进入分析。# 用 GDAL 查看原始文件的元数据确认覆盖范围和 nodata 设置 gdalinfo ASTGTM_N30E090_dem.tif # 统计每个波段的数值分布重点关注 min/max 是否异常 gdalinfo -stats ASTGTM_N30E090_dem.tif # 用 Python 快速检查 nodata 像素的占比判断空洞是否严重 python3 -c from osgeo import gdal import numpy as np ds gdal.Open(ASTGTM_N30E090_dem.tif) band ds.GetRasterBand(1) arr band.ReadAsArray() nodata band.GetNoDataValue() if nodata is not None: mask (arr nodata) print(fnodata 占比: {mask.sum() / arr.size:.2%}) else: print(未设置 nodata检查极端值: min{np.nanmin(arr)}, max{np.nanmax(arr)}) 上面的逻辑是先用gdalinfo查看栅格的整体元数据包括尺寸、波段、坐标系和nodata标记再用-stats触发GDAL计算统计信息最后用Python原生读取数组统计空洞比例。如果nodata占比在0.1%以下通常可以直接用邻域插值补洞如果超过1%比如某些山地区域SRTM的V2版本仍有残留空洞我一般会换数据源而不是强行插值因为空洞区域往往地形破碎插值结果完全不可信。值得留意的是不同下载站点的SRTM版本可能有差异——USGS的版本通常已经做过空洞填充而某些镜像站提供的是原始版本这也是同一区域不同数据看起来差异很大的一个重要原因。3. 中国范围DEM拼接、裁剪与投影转换的实际操作3.1 多幅DEM栅格数据的高效拼接国内下载DEM的习惯是按分幅拉取一次项目可能要下载几十甚至上百幅。GDAL提供两条拼接路径gdal_merge.py适合幅数少、无需过度优化的场景gdalbuildvrt先生成虚拟栅格再转GeoTIFF则在大幅数据时优势明显因为VRT不复制像素数据只记录映射关系。拼接前需要先确认所有分幅的坐标系、分辨率和nodata设置是否一致不一致时gdal_merge.py虽然会自动重投影但速度极慢且插值算法默认是最近邻精度不可控。# 先用 gdalinfo 快速列出所有文件的坐标系检查一致性 for f in *.tif; do gdalinfo $f | grep -E ([0-9]{4,6}); done # 推荐方式先建 VRT再转为正式 GeoTIFF gdalbuildvrt -resolution highest -r bilinear all_dem.vrt *.tif # 将 VRT 烘焙成带压缩的 GeoTIFF适合后续频繁读写 gdal_translate -co COMPRESSDEFLATE -co BIGTIFFYES -co TILEDYES all_dem.vrt china_dem.tif这里解释一下几个关键参数-resolution highest保证各分幅采样到一致分辨率选最高而非最低是为了避免后续分析时重采样信息损失-r bilinear使用双线性插值处理边缘接缝比默认的最近邻平滑但要注意在陡峭地形可能引入轻微模糊BIGTIFFYES强制使用BigTIFF格式因为中国全境的30米分辨率DEM拼接后轻松超过4GB传统TIFF格式根本存不下TILEDYES开启内部瓦片化后续按需读取小区域时性能提升明显。3.2 按行政边界或矢量范围裁剪DEM数据拿到拼接后的全国DEM通常需要按行政区或研究区裁剪。gdalwarp配合行政区shp文件是最主流的做法但有几个容易被忽略的细节。# 按 shp 范围裁剪并同时完成投影转换 gdalwarp -cutline china_dem.shp -crop_to_cutline -dstnodata -9999 \ -t_srs EPSG:4524 -tr 30 30 -r bilinear -co COMPRESSDEFLATE \ input_china_dem.tif output_province_dem.tif-crop_to_cutline以shp的外边界对齐裁剪结果的范围和矢量边界完全一致-dstnodata -9999显式设置输出nodata值防止原始文件里不同nodata值混入结果-t_srs EPSG:4524把WGS84地理坐标转为CGCS2000高斯投影-tr 30 30同时重新采样到30米分辨率。一个常见误区是只加-cutline不加-crop_to_cutline这样输出范围仍是全图矩形只是矩形内边界外的区域被设为nodata文件体积不会减小后续处理还得靠掩膜操作。3.3 不同分辨率DEM之间的重采样策略ALOS 12.5米DEM要参与一个设计分辨率30米的分析流程或者把30米DEM重采样到90米这是DEM栅格数据预处理的典型需求。重采样算法选型有讲究最邻近插值只在分类栅格或坡度等级数据中使用避免引入不存在的值双线性插值通用首选地形平滑且计算快适合一般分析三次卷积插值需要更精细化地形时用典型场景是水文分析中的流向提取但会产生超出原始值范围的插值结果在洼地填充阶段可能引入人工洼地聚合方法从高分辨率聚合成低分辨率时如12.5米到30米使用gdalwarp的-r average做像素聚合得到的结果比任何插值都更物理。# 将 12.5 米 ALOS DEM 聚合重采样到 30 米使用平均聚合 gdalwarp -tr 30 30 -r average -co COMPRESSDEFLATE alos_12_5m.tif alos_30m.tif这个场景我踩过一次坑用双线性插值把12.5米DEM降到30米后生成的坡度图在山脊线附近有规律的条带状伪影检查后发现是重采样时栅格对齐偏移导致的。解决方法是先用gdal_align或gdalwarp -te把网格对齐到目标数据的边界再做重采样。4. 基于中国DOM的典型地形分析与DSM生成DEM的流程4.1 用GDAL计算坡度、坡向与曲率完成预处理的DEM栅格数据下一步通常是坡度、坡向和曲率的地形指标计算。GDAL的gdaldem工具支持slope、aspect、hillshade、TRI地形粗糙度指数、TPI地形位置指数等常用地形指标但国内项目里真正的高频需求是TWI地形湿度指数或汇流累积量这类复合指标需要依赖WhiteboxTools或TauDEM不过基础的坡度坡向计算用gdaldem完全足够。# 坡度单位度 gdaldem slope input_dem.tif output_slope.tif -p -s 111120 # 坡向按罗盘方位角输出 gdaldem aspect input_dem.tif output_aspect.tif -zero_for_flat # 地形位置指数 TPI用于划分山脊、山谷与平地 gdaldem TPI input_dem.tif output_tpi.tif上面-p表示输出坡度为百分比还是度不加法默认输出度-s 111120是指定水平和垂直方向的比例关系这里把1度约等于111120米换算传入原因是当输入DEM坐标系是地理坐标时GDAL需要知道X和Y方向的真实尺度才能算出物理意义上的坡度。如果DEM已经投影到高斯或UTM坐标系-s参数不是必须的GDAL能直接从投影定义中获取单位。-zero_for_flat把平坦区域的坡向设为0度否则平坦区域因为微小噪声会产生随机坡向影响后续分析。4.2 计算TWI地形湿度指数与汇流累积量水文分析是国内DEM应用的高频场景从原始DEM到汇流累积量通常需要完成填洼、流向、累积三步。填洼这一步在水文分析里的地位被很多初学者低估——原始DEM中存在的伪洼地尤其是ASTER GDEM的麻点噪声如果不处理流向会被截断汇流累积量结果严重失真。常用的流程是使用WhiteboxTools因为它既可以命令行操作又天然支持中文区域大数据量处理。# 填洼提升阈值法避免过度填洼破坏真实地形 whitebox_tools -r FillDepressions -p --deminput_dem.tif --outputdem_filled.tif --fix_flatstrue # D8 流向算法 whitebox_tools -r D8Pointer -p --demdem_filled.tif --outputd8pointer.tif # D8 汇流累积量 whitebox_tools -r D8FlowAccumulation -p --demdem_filled.tif --outputd8accum.tif --out_typesca实际项目中比较棘手的是填洼参数的选择。FillDepressions默认填充到最大洼地深度这对真实存在的封闭盆地如塔里木盆地局部区域会过度修改地形。我一般先做一次快速统计看洼地深度分布然后在--max_depth参数里限制最大填充深度否则计算结果里那些本来是真实盆地的区域完全被淹没在平坦面里。4.3 用DSM生成DEM植被与建筑物剔除ALOS 12.5米数据本质是DSM之前提到过。那在国内实践中“DSM生成DEM”到底是个什么流程业内通常有两种理解一种是对原始点云或者DSM栅格做形态学滤波剔除植被和建筑物另一种是用DSM减去对应区域的地物高度模型nDSM恢复地面高程。第一种更适合大范围因为不需要额外的高精度地物高度数据。# 使用 GDAL 的形态学开运算做初步平滑剔除细微的植被冠层噪声 gdal_fillnodata.py -md 10 -si 1 input_dsm.tif output_dem_smooth.tif # 对城市区域做差值分析提取建筑物轮廓DEM与DSM之差 gdal_calc.py -A input_dsm.tif -B input_dem.tif --outfileheight_diff.tif \ --calc(A-B) --NoDataValue-9999形态学开运算的思路是先腐蚀再膨胀能把小于结构元素尺寸的凸起如树冠、小建筑削平同时保留大型地形起伏。城市区域的中高层建筑物宽度远大于结构元素形态学滤波效果有限更好的做法是用现有DSM和DEM相减提取建筑高度再对建筑区域做插值替换。国内一些数据集就专门提供“我国地下水位栅格数据shp”之类的分层产品这也说明DEM相关数据在国内正在向综合数据产品发展。5. 全国范围栅格数据在内存与磁盘上的性能优化5.1 大数据量DEM的分块读写策略拼接完成的全国30米DEM文件动辄50100GB读取全量进内存不仅没必要而且几乎不可行。处理这种数据的正确姿势是用Block块方式读取配合GDAL的ReadAsArray指定窗口。from osgeo import gdal import numpy as np ds gdal.Open(china_dem_30m.tif) band ds.GetRasterBand(1) # 获取栅格尺寸 xsize band.XSize ysize band.YSize # 分块大小根据内存余量调整一般 2048x2048 较为合适 block_size 2048 for y in range(0, ysize, block_size): rows min(block_size, ysize - y) for x in range(0, xsize, block_size): cols min(block_size, xsize - x) data band.ReadAsArray(x, y, cols, rows) # 这里做具体的分析逻辑比如计算坡度、提取阈值等 print(f处理块: ({x}, {y})-({xcols}, {yrows}))Python的range遍历块的逻辑保证每个循环取到的rows和cols不会超出栅格边界这是分块处理容易犯的越界错误的标准规避方式。需要留意的是ReadAsArray默认返回浮点数组对整型DEM可以在读取时指定band.ReadAsArray(x, y, cols, rows, buf_obj...)复用缓冲区减少内存分配次数。对于需要做邻域分析的算法如坡度计算分块时每边多读一个radius像素halo处理完再丢弃这就是带缓冲的分块策略。5.2 压缩格式选型DEFLATE与LZMA的取舍中国范围DEM栅格数据的压缩选型对磁盘占用和运行速度影响很大。GeoTIFF内部压缩有DEFLATE和LZMA两种主流选择。DEFLATE在速度与压缩率之间平衡最好绝大多数场景推荐使用LZMA压缩率高出约15%20%但压缩和解压时CPU消耗翻倍。对于需要反复读取的底图数据我建议用DEFLATE搭配PREDICTOR2参数后者对浮点栅格启用水平差分预测能显著提高压缩率且解码开销很低。gdal_translate -co COMPRESSDEFLATE -co PREDICTOR2 -co TILEDYES \ -co BLOCKXSIZE512 -co BLOCKYSIZE512 input.tif output_deflate.tifBLOCKXSIZE和BLOCKYSIZE都是512表示瓦片尺寸512×512像素。GDAL默认的块大小是128×128但国内范围大图的随机读取场景下512×512能减少I/O次数实测在机械硬盘上全图拼接速度提升尤其明显。还有一个隐形坑如果原始DEM带有内部金字塔overview要先确认金字塔也是压缩的否则文件压缩了但金字塔占用的空间仍然巨大。5.3 DEM缓存金字塔与可视化提速大范围DEM在Web端展示或桌面端交互浏览时直接渲染全分辨率数据必然卡顿。标准做法是构建内部概览金字塔overviews。# 为已有的 GeoTIFF 添加金字塔4级每级缩小1/4面积 gdaladdo -r average china_dem.tif 2 4 8 16 32-r average用平均法重采样生成金字塔级别的数据适合DEM这类连续表面数据。如果不加-r默认使用最近邻山脊线会有明显锯齿。检查金字塔是否生效可以用gdalinfo china_dem.tif | grep Overviews确认没有丢失。实际操作中我发现许多从地理空间数据云下载的DEM已经自带金字塔直接执行gdaladdo会追加而不是覆盖反而浪费空间。正确做法是先gdalinfo查看已有金字塔级别再补充缺少的级别。6. 融合多源数据填补空洞的一个可复现技巧不同源DEM数据在中国不同区域有各自的优劣一个实用技巧是用区域权重融合SRTM和ALOS数据在平坦地区信任ALOS 12.5米的高分辨率细节在陡峭山区信任SRTM的填充完整性。具体做法是以地形坡度作为权重依据坡度小的区域给ALOS更大的权重坡度大的区域给SRTM更大权重。# 第一步计算坡度的绝对值作为权重参考的倒数 gdaldem slope alos_12_5m.tif slope_alos.tif gdaldem slope srtm_30m.tif slope_srtm.tif # 第二步用 gdal_calc 做加权融合 gdal_calc.py -A alos_12_5m.tif -B srtm_30m.tif -C slope_alos.tif --outfilefused_dem.tif \ --calcA * (C 15) B * (C 15) --NoDataValue-9999这段代码的逻辑是凡是ALOS数据的坡度大于15度的地方保留ALOS原始值坡度小于等于15度的平坦区域取SRTM的值。之所以在陡峭区域保留ALOS是因为ALOS的立体成像在陡峭地形虽然存在部分噪声但整体纹理细节保留远好于SRTM的C波段数据而平坦区域SRTM不会出现空洞数值更稳定。用gdal_calc的布尔表达式向量化操作可以处理大范围栅格而不需要分块循环。验证融合结果的正确性可以用两种手段一是对融合前后坡度直方图做对比看是否出现明显的双峰或中空分布异常二是提取融合后DEM沿某条已知河流断面的高程和实测水面比降对比较。这两个验证动作比任何抽象评价指标都能快速暴露融合参数不当的问题。本文还有配套的精品资源点击获取 中国范围DEM栅格数据全指南：从数据源选型到地形分析落地
简介一份覆盖中国全境的DEM栅格数据集分辨率为1km源自SRTM 30米分辨率原始数据经ArcGIS镶嵌拼接与重采样制作完成数据完整、边界准确可作为多种空间分析的统一底图。面向地理信息科学、测绘工程等专业学生以及需要全国尺度地形底图的科研人员可用于课程设计、地形因子提取、宏观区域分析和多源数据叠加等场景。数据采用WGS-84坐标系TIFF格式无需转换即可在ArcGIS中直接打开方便坡度坡向提取、可视域分析、栅格计算及制图出图等常见操作。压缩包共7个文件、26.51MB核心为ChinaDEM.tif栅格主数据同时收录tfw坐标配准文件、ovr金字塔文件、xml元数据及dbf属性表等辅助文件便于快速加载和二次处理。目前已有693人学习下载既能支撑教学演示也可作为科研前期数据准备数据精度经过重采样优化整体实用性强。1. 中国范围的DEM栅格数据从数据源选型到地形分析落地做全国性或区域性地形分析时第一道坎不是算法而是数据本身。中国范围的DEM栅格数据看似随处可下但不同来源的精度、坐标系、覆盖范围差异极大——SRTM在新疆某些区域存在空洞ASTER GDEM在国内高山区域噪声明显ALOS 12.5米数据精度高却动辄几十GB。更麻烦的是从地理空间数据云下载的DEM通常是分幅存储的坐标参考系为WGS84而国内项目往往要求CGCS2000或高斯投影直接拿来用轻则投影变形重则分析结果完全偏离。这篇文章就想把“获取中国范围DEM栅格数据到完成地形指标计算”这条链路拆开来讲覆盖数据源横向对比、预处理参数、地形指标计算和大范围数据的性能瓶颈给出一套可以直接照做的方案。适合的读者群很明确GIS开发、遥感应用工程师、水文或规划领域的分析人员以及做三维可视化或DSM生成DEM相关工作的朋友。无论你只需要某个省份的裁剪DEM数据还是想拼出全国一张图下面的思路都适用。2. DEM数据源选型与中国区域覆盖精度、坐标系与空洞问题2.1 三大主流DEM数据源在中国区域的真实表现覆盖中国范围的全球DEM数据源中SRTM、ASTER GDEM和ALOS AW3D30是实际项目里最常被用到的。SRTM由NASA和NGA联合获取在中国区域主要覆盖北纬60°以南国内绝大部分领土都在范围内分辨率分30米和90米两档其中90米版本在国内用得最广文件组织为分幅TIF单幅约5°×5°。ASTER GDEM由日本METI和美国NASA基于光学立体像对生成分辨率30米覆盖纬度83°N到83°S中国全域可覆盖但光学影像受云影响较大在青藏高原、天山等高山区域经常出现异常凹陷或尖峰业内俗称“麻点”噪声。ALOS AW3D30是JAXA基于PRISM立体传感器生成的全球DSM标称30米实际在绝大多数区域能体现更细致的地形纹理。而热词里常被提起的ALOS 12.5米DEM数据其实来源于ALOS PALSAR的12.5米分辨率产品精度确实高出其他全球产品一个量级但中国区域需要分块下载且原始数据为DSM——包含地表覆盖物高度不是剔除了植被和建筑物的纯地面DEM。基于这个差别做洪水淹没模拟或耕地平整度分析前需要先明确自己到底要DSM还是DEM。用DSM生成DEM的做法通常是对原始DSM做滤波剔除植被冠层和建筑物高度这在城区和密林区域效果有限但在地形起伏主导的区域表现尚可。2.2 坐标系差异与国内常用投影的选择下载到的全球DEM原始数据坐标系基本都是WGS84地理坐标系单位为度。很多初学者直接用原始数据做坡度计算或面积统计得到的结果毫无意义——坡度在经纬度坐标系下X和Y方向长度不等得到的坡度和实际物理坡度存在系统性偏差。处理中国范围的DEM栅格数据时第一步应该是确认目标坐标系。国内生产项目常用的有三类CGCS2000地理坐标系代号EPSG:4490适合全国范围展示和简单分析CGCS2000 / 3-degree Gauss-Kruger zone代号EPSG:4524到4559按经度分带适合省级或地市级项目UTM分区代号EPSG:32649到32652覆盖中国从东到西的分区如果上下游工具链是国际通用软件UTM反而更省事。需要特别注意CGCS2000高斯投影分带规则和WGS84 UTM分带规则不同高斯3度带从东经75°起每隔3°一带中央经线为75、78、81……直到135UTM则从西经180°起每隔6°一带中央经线为75、81、87……285。对于范围跨多个分带的全国数据正确做法是保留地理坐标做全国分析只在局部高精度分析时才投影到对应分带。某些省级项目直接把数据投影到CGCS2000 / Gauss-Kruger zone 20中央经线117°处理全省范围投影变形会导致东西跨度大的省份边缘区域精度损失。2.3 数据空洞和异常值的识别的快速检查法中国区域DEM一个老生常谈的痛点是SRTM在陡峭地形的数据空洞以及ASTER GDEM在湖泊、积雪区域的异常值。拿到数据后先做概览检查别直接进入分析。# 用 GDAL 查看原始文件的元数据确认覆盖范围和 nodata 设置 gdalinfo ASTGTM_N30E090_dem.tif # 统计每个波段的数值分布重点关注 min/max 是否异常 gdalinfo -stats ASTGTM_N30E090_dem.tif # 用 Python 快速检查 nodata 像素的占比判断空洞是否严重 python3 -c from osgeo import gdal import numpy as np ds gdal.Open(ASTGTM_N30E090_dem.tif) band ds.GetRasterBand(1) arr band.ReadAsArray() nodata band.GetNoDataValue() if nodata is not None: mask (arr nodata) print(fnodata 占比: {mask.sum() / arr.size:.2%}) else: print(未设置 nodata检查极端值: min{np.nanmin(arr)}, max{np.nanmax(arr)}) 上面的逻辑是先用gdalinfo查看栅格的整体元数据包括尺寸、波段、坐标系和nodata标记再用-stats触发GDAL计算统计信息最后用Python原生读取数组统计空洞比例。如果nodata占比在0.1%以下通常可以直接用邻域插值补洞如果超过1%比如某些山地区域SRTM的V2版本仍有残留空洞我一般会换数据源而不是强行插值因为空洞区域往往地形破碎插值结果完全不可信。值得留意的是不同下载站点的SRTM版本可能有差异——USGS的版本通常已经做过空洞填充而某些镜像站提供的是原始版本这也是同一区域不同数据看起来差异很大的一个重要原因。3. 中国范围DEM拼接、裁剪与投影转换的实际操作3.1 多幅DEM栅格数据的高效拼接国内下载DEM的习惯是按分幅拉取一次项目可能要下载几十甚至上百幅。GDAL提供两条拼接路径gdal_merge.py适合幅数少、无需过度优化的场景gdalbuildvrt先生成虚拟栅格再转GeoTIFF则在大幅数据时优势明显因为VRT不复制像素数据只记录映射关系。拼接前需要先确认所有分幅的坐标系、分辨率和nodata设置是否一致不一致时gdal_merge.py虽然会自动重投影但速度极慢且插值算法默认是最近邻精度不可控。# 先用 gdalinfo 快速列出所有文件的坐标系检查一致性 for f in *.tif; do gdalinfo $f | grep -E ([0-9]{4,6}); done # 推荐方式先建 VRT再转为正式 GeoTIFF gdalbuildvrt -resolution highest -r bilinear all_dem.vrt *.tif # 将 VRT 烘焙成带压缩的 GeoTIFF适合后续频繁读写 gdal_translate -co COMPRESSDEFLATE -co BIGTIFFYES -co TILEDYES all_dem.vrt china_dem.tif这里解释一下几个关键参数-resolution highest保证各分幅采样到一致分辨率选最高而非最低是为了避免后续分析时重采样信息损失-r bilinear使用双线性插值处理边缘接缝比默认的最近邻平滑但要注意在陡峭地形可能引入轻微模糊BIGTIFFYES强制使用BigTIFF格式因为中国全境的30米分辨率DEM拼接后轻松超过4GB传统TIFF格式根本存不下TILEDYES开启内部瓦片化后续按需读取小区域时性能提升明显。3.2 按行政边界或矢量范围裁剪DEM数据拿到拼接后的全国DEM通常需要按行政区或研究区裁剪。gdalwarp配合行政区shp文件是最主流的做法但有几个容易被忽略的细节。# 按 shp 范围裁剪并同时完成投影转换 gdalwarp -cutline china_dem.shp -crop_to_cutline -dstnodata -9999 \ -t_srs EPSG:4524 -tr 30 30 -r bilinear -co COMPRESSDEFLATE \ input_china_dem.tif output_province_dem.tif-crop_to_cutline以shp的外边界对齐裁剪结果的范围和矢量边界完全一致-dstnodata -9999显式设置输出nodata值防止原始文件里不同nodata值混入结果-t_srs EPSG:4524把WGS84地理坐标转为CGCS2000高斯投影-tr 30 30同时重新采样到30米分辨率。一个常见误区是只加-cutline不加-crop_to_cutline这样输出范围仍是全图矩形只是矩形内边界外的区域被设为nodata文件体积不会减小后续处理还得靠掩膜操作。3.3 不同分辨率DEM之间的重采样策略ALOS 12.5米DEM要参与一个设计分辨率30米的分析流程或者把30米DEM重采样到90米这是DEM栅格数据预处理的典型需求。重采样算法选型有讲究最邻近插值只在分类栅格或坡度等级数据中使用避免引入不存在的值双线性插值通用首选地形平滑且计算快适合一般分析三次卷积插值需要更精细化地形时用典型场景是水文分析中的流向提取但会产生超出原始值范围的插值结果在洼地填充阶段可能引入人工洼地聚合方法从高分辨率聚合成低分辨率时如12.5米到30米使用gdalwarp的-r average做像素聚合得到的结果比任何插值都更物理。# 将 12.5 米 ALOS DEM 聚合重采样到 30 米使用平均聚合 gdalwarp -tr 30 30 -r average -co COMPRESSDEFLATE alos_12_5m.tif alos_30m.tif这个场景我踩过一次坑用双线性插值把12.5米DEM降到30米后生成的坡度图在山脊线附近有规律的条带状伪影检查后发现是重采样时栅格对齐偏移导致的。解决方法是先用gdal_align或gdalwarp -te把网格对齐到目标数据的边界再做重采样。4. 基于中国DOM的典型地形分析与DSM生成DEM的流程4.1 用GDAL计算坡度、坡向与曲率完成预处理的DEM栅格数据下一步通常是坡度、坡向和曲率的地形指标计算。GDAL的gdaldem工具支持slope、aspect、hillshade、TRI地形粗糙度指数、TPI地形位置指数等常用地形指标但国内项目里真正的高频需求是TWI地形湿度指数或汇流累积量这类复合指标需要依赖WhiteboxTools或TauDEM不过基础的坡度坡向计算用gdaldem完全足够。# 坡度单位度 gdaldem slope input_dem.tif output_slope.tif -p -s 111120 # 坡向按罗盘方位角输出 gdaldem aspect input_dem.tif output_aspect.tif -zero_for_flat # 地形位置指数 TPI用于划分山脊、山谷与平地 gdaldem TPI input_dem.tif output_tpi.tif上面-p表示输出坡度为百分比还是度不加法默认输出度-s 111120是指定水平和垂直方向的比例关系这里把1度约等于111120米换算传入原因是当输入DEM坐标系是地理坐标时GDAL需要知道X和Y方向的真实尺度才能算出物理意义上的坡度。如果DEM已经投影到高斯或UTM坐标系-s参数不是必须的GDAL能直接从投影定义中获取单位。-zero_for_flat把平坦区域的坡向设为0度否则平坦区域因为微小噪声会产生随机坡向影响后续分析。4.2 计算TWI地形湿度指数与汇流累积量水文分析是国内DEM应用的高频场景从原始DEM到汇流累积量通常需要完成填洼、流向、累积三步。填洼这一步在水文分析里的地位被很多初学者低估——原始DEM中存在的伪洼地尤其是ASTER GDEM的麻点噪声如果不处理流向会被截断汇流累积量结果严重失真。常用的流程是使用WhiteboxTools因为它既可以命令行操作又天然支持中文区域大数据量处理。# 填洼提升阈值法避免过度填洼破坏真实地形 whitebox_tools -r FillDepressions -p --deminput_dem.tif --outputdem_filled.tif --fix_flatstrue # D8 流向算法 whitebox_tools -r D8Pointer -p --demdem_filled.tif --outputd8pointer.tif # D8 汇流累积量 whitebox_tools -r D8FlowAccumulation -p --demdem_filled.tif --outputd8accum.tif --out_typesca实际项目中比较棘手的是填洼参数的选择。FillDepressions默认填充到最大洼地深度这对真实存在的封闭盆地如塔里木盆地局部区域会过度修改地形。我一般先做一次快速统计看洼地深度分布然后在--max_depth参数里限制最大填充深度否则计算结果里那些本来是真实盆地的区域完全被淹没在平坦面里。4.3 用DSM生成DEM植被与建筑物剔除ALOS 12.5米数据本质是DSM之前提到过。那在国内实践中“DSM生成DEM”到底是个什么流程业内通常有两种理解一种是对原始点云或者DSM栅格做形态学滤波剔除植被和建筑物另一种是用DSM减去对应区域的地物高度模型nDSM恢复地面高程。第一种更适合大范围因为不需要额外的高精度地物高度数据。# 使用 GDAL 的形态学开运算做初步平滑剔除细微的植被冠层噪声 gdal_fillnodata.py -md 10 -si 1 input_dsm.tif output_dem_smooth.tif # 对城市区域做差值分析提取建筑物轮廓DEM与DSM之差 gdal_calc.py -A input_dsm.tif -B input_dem.tif --outfileheight_diff.tif \ --calc(A-B) --NoDataValue-9999形态学开运算的思路是先腐蚀再膨胀能把小于结构元素尺寸的凸起如树冠、小建筑削平同时保留大型地形起伏。城市区域的中高层建筑物宽度远大于结构元素形态学滤波效果有限更好的做法是用现有DSM和DEM相减提取建筑高度再对建筑区域做插值替换。国内一些数据集就专门提供“我国地下水位栅格数据shp”之类的分层产品这也说明DEM相关数据在国内正在向综合数据产品发展。5. 全国范围栅格数据在内存与磁盘上的性能优化5.1 大数据量DEM的分块读写策略拼接完成的全国30米DEM文件动辄50100GB读取全量进内存不仅没必要而且几乎不可行。处理这种数据的正确姿势是用Block块方式读取配合GDAL的ReadAsArray指定窗口。from osgeo import gdal import numpy as np ds gdal.Open(china_dem_30m.tif) band ds.GetRasterBand(1) # 获取栅格尺寸 xsize band.XSize ysize band.YSize # 分块大小根据内存余量调整一般 2048x2048 较为合适 block_size 2048 for y in range(0, ysize, block_size): rows min(block_size, ysize - y) for x in range(0, xsize, block_size): cols min(block_size, xsize - x) data band.ReadAsArray(x, y, cols, rows) # 这里做具体的分析逻辑比如计算坡度、提取阈值等 print(f处理块: ({x}, {y})-({xcols}, {yrows}))Python的range遍历块的逻辑保证每个循环取到的rows和cols不会超出栅格边界这是分块处理容易犯的越界错误的标准规避方式。需要留意的是ReadAsArray默认返回浮点数组对整型DEM可以在读取时指定band.ReadAsArray(x, y, cols, rows, buf_obj...)复用缓冲区减少内存分配次数。对于需要做邻域分析的算法如坡度计算分块时每边多读一个radius像素halo处理完再丢弃这就是带缓冲的分块策略。5.2 压缩格式选型DEFLATE与LZMA的取舍中国范围DEM栅格数据的压缩选型对磁盘占用和运行速度影响很大。GeoTIFF内部压缩有DEFLATE和LZMA两种主流选择。DEFLATE在速度与压缩率之间平衡最好绝大多数场景推荐使用LZMA压缩率高出约15%20%但压缩和解压时CPU消耗翻倍。对于需要反复读取的底图数据我建议用DEFLATE搭配PREDICTOR2参数后者对浮点栅格启用水平差分预测能显著提高压缩率且解码开销很低。gdal_translate -co COMPRESSDEFLATE -co PREDICTOR2 -co TILEDYES \ -co BLOCKXSIZE512 -co BLOCKYSIZE512 input.tif output_deflate.tifBLOCKXSIZE和BLOCKYSIZE都是512表示瓦片尺寸512×512像素。GDAL默认的块大小是128×128但国内范围大图的随机读取场景下512×512能减少I/O次数实测在机械硬盘上全图拼接速度提升尤其明显。还有一个隐形坑如果原始DEM带有内部金字塔overview要先确认金字塔也是压缩的否则文件压缩了但金字塔占用的空间仍然巨大。5.3 DEM缓存金字塔与可视化提速大范围DEM在Web端展示或桌面端交互浏览时直接渲染全分辨率数据必然卡顿。标准做法是构建内部概览金字塔overviews。# 为已有的 GeoTIFF 添加金字塔4级每级缩小1/4面积 gdaladdo -r average china_dem.tif 2 4 8 16 32-r average用平均法重采样生成金字塔级别的数据适合DEM这类连续表面数据。如果不加-r默认使用最近邻山脊线会有明显锯齿。检查金字塔是否生效可以用gdalinfo china_dem.tif | grep Overviews确认没有丢失。实际操作中我发现许多从地理空间数据云下载的DEM已经自带金字塔直接执行gdaladdo会追加而不是覆盖反而浪费空间。正确做法是先gdalinfo查看已有金字塔级别再补充缺少的级别。6. 融合多源数据填补空洞的一个可复现技巧不同源DEM数据在中国不同区域有各自的优劣一个实用技巧是用区域权重融合SRTM和ALOS数据在平坦地区信任ALOS 12.5米的高分辨率细节在陡峭山区信任SRTM的填充完整性。具体做法是以地形坡度作为权重依据坡度小的区域给ALOS更大的权重坡度大的区域给SRTM更大权重。# 第一步计算坡度的绝对值作为权重参考的倒数 gdaldem slope alos_12_5m.tif slope_alos.tif gdaldem slope srtm_30m.tif slope_srtm.tif # 第二步用 gdal_calc 做加权融合 gdal_calc.py -A alos_12_5m.tif -B srtm_30m.tif -C slope_alos.tif --outfilefused_dem.tif \ --calcA * (C 15) B * (C 15) --NoDataValue-9999这段代码的逻辑是凡是ALOS数据的坡度大于15度的地方保留ALOS原始值坡度小于等于15度的平坦区域取SRTM的值。之所以在陡峭区域保留ALOS是因为ALOS的立体成像在陡峭地形虽然存在部分噪声但整体纹理细节保留远好于SRTM的C波段数据而平坦区域SRTM不会出现空洞数值更稳定。用gdal_calc的布尔表达式向量化操作可以处理大范围栅格而不需要分块循环。验证融合结果的正确性可以用两种手段一是对融合前后坡度直方图做对比看是否出现明显的双峰或中空分布异常二是提取融合后DEM沿某条已知河流断面的高程和实测水面比降对比较。这两个验证动作比任何抽象评价指标都能快速暴露融合参数不当的问题。本文还有配套的精品资源点击获取