云南30米DEM数据处理：从数据识别到地形分析的完整流程

简介：云南省30米分辨率的数字高程模型数据，面向地理信息系统从业者与城市规划、环境研究等专业人士，为区域地形分析提供精细基础。压缩包共17个文件，以Shapefile格式为核心，包含几何、属性、投影定义等关键文件，辅以索引和元数据文件，可提升空间查询与读取效率，压缩后体积仅43.66MB，能在ArcGIS、QGIS等平台直接打开。目前已有1178人学习下载，数据实用性获得认可。借助该数据可进行坡度坡向计算、流域划分、洪水淹没模拟、视线分析、地形剖面提取，支撑城市扩张评估、交通路线规划、灾害风险区划、生态保护及水利工程选址等工作；30米分辨率兼顾精度与数据量，适合区域级地形研究，也能融合土地利用、人口分布等多源数据，为综合空间建模与决策支持提供可靠数据底座。

1. “各省30米分辨率DEM数据-云南.rar”到底是什么——先搞清你拿到的是哪类 DEM

拿到“各省30米分辨率DEM数据-云南.rar”，第一反应不应该是双击解压，而是先确认数据源和分辨率。30米分辨率DEM通常指 1 弧秒采样的高程格网，云南省范围内意味着横断山的深切割河谷、滇东喀斯特峰丛和滇中坝子都在同一套网格里。做地形分析、流域划分、选址评估和山地灾害评价的人，最需要的就是把全球分幅数据按行政区统一成连续可用产品。这类压缩包打开后常见是一批以图幅号命名的 GeoTIFF 或 HGT，而不是一张整图的云南 DEM。先花十分钟把坐标系、无效值和分幅范围确认清楚，能让后面的拼接和出图少踩大量坑。

2. 解开 .rar 后的数据识别：从文件名、元数据判断 30 米 DEM 的版本与坐标系

2.1 文件命名与分幅规则先告诉你版本

解压后我一般先不打开 GIS，而是先看一眼目录里的文件名。30 米级 DEM 绝大多数采用 1°×1° 分幅，这决定了解压后你会看到一批以经纬度编号命名的文件，而不是一张直接覆盖全省的大图。比如N24E100_dem.tif表示北界 24°N、东界 100°E 的那一格；N24E100.hgt则是 SRTM 原始格式的同名图幅。云南大致覆盖 N21～N29、E97～E106 这 9 个纬度带和 10 个经度带，所以一个按省整理好的 rar 里应有 80～90 个分幅文件。如果只有十几张，说明数据源已经被人为合并或按地级市切过，后续拼图方式要调整。

hgt后缀比较特殊。SRTM 的 hgt 固定是 3601×3601 的 int16 网格，1 弧秒采样，文件里不包含坐标系信息和无效值标记，GDAL 打开后把 -32768 当作空洞是常见行为。GeoTIFF 后缀则要看波段数：单波段 Int16/Float32 才是可用于高程分析的 DEM，三波段 Byte 的绝大多数是别人渲染过的晕渲图。

2.2 用 GDAL 批量读取元数据，先给数据做“体检”

拿到几十个 tif 后，第一件事是跑一组只输出关键信息的 gdalinfo，确认所有分幅的分辨率、像元大小和投影一致。命令如下：

for f in *.tif; do echo "== $f ==" gdalinfo "$f" | grep -E "Size is|Origin|Pixel Size|EPSG" done

这段循环把每个文件的尺寸、左上角坐标、像元大小和 EPSG 编号打出来。重点看 Pixel Size 的数值是否是0.00027777778,-0.00027777778。这是 1/3600 度的十进制表示，也就是 1 弧秒，约 30 米。如果出现0.00083333333之类，说明混入了 3 弧秒（约 90 米）数据，要标记出来，否则后续拼接会降低整片区域的细节。

在 Python 里同样可以批量判断投影，适合把结果写进日志或直接做异常检测：

from osgeo import gdal, osr import glob for fp in sorted(glob.glob("*.tif")): ds = gdal.Open(fp) gt = ds.GetGeoTransform() proj_wkt = ds.GetProjection() epsg = None if proj_wkt: srs = osr.SpatialReference() srs.ImportFromWkt(proj_wkt) epsg = srs.GetAuthorityCode(None) print(f"{fp:26s} {ds.RasterXSize}x{ds.RasterYSize} " f"origin=({gt[0]:.5f},{gt[3]:.5f}) " f"pixel=({gt[1]:.8f},{gt[5]:.8f}) epsg={epsg}") ds = None

这段代码用GetGeoTransform()读取六参数仿射变换，前两个坐标给左上角，像元尺寸用第 2 和第 6 个参数表示，负值表示数据按行向下读取。GetAuthorityCode(None)从投影定义里抽出 EPSG 编号。如果打印结果是 4326，数据和省界都是经纬度坐标，后面按省界裁剪时无需转换；如果是投影坐标，就要先把省界转成同一坐标系再裁剪，不能硬套。

2.3 不同 30 米 DEM 的来源差异：一张表看清选择和坑

网上能下载到的“30米DEM”通常来自四个主流产品：ASTER GDEM V3、SRTM V3 30m、ALOS AW3D30 和 Copernicus GLO-30。它们的标称分辨率都是 30 米级，但传感器、后处理策略和高程噪声特征差别很大。在云南做地形分析时选哪个，决定了后面参数要不要做针对性修正。

判断无效值最直接的办法是看gdalinfo -stats输出的最小值。正常云南陆地区域没有负海拔，最低点在南溪河出境处约 76 米。如果最小值出现 -9999 或 -32768，那它一定是无效值标记，而不能当成真实高程。ASTER 与 SRTM 对无效值的标记不同，统一转成同一个 nodata 值，才能在填洼、坡度计算时不把“空洞”当成“洼地”。

2.4 不要用带渲染色的 DEM 底图替代高程数据

一个常见的坑是 rar 解压后出现三个波段的彩色 tif。这个文件可能已经被人为渲染过高程渐变和山体阴影，视觉上很漂亮，但每个像元存的是 RGB 颜色而不是高程值。判断方式很简单：gdalinfo里 Band Count 是 3，Type 是 Byte。高程分析的 DEM 至少是 Int16，标准做法是 Float32。如果包里全是这种渲染图，就不要强行转格式骗自己，直接换用原始高程数据源。渲染图唯一适合的用途是当底图背景，不能用来计算坡度、填洼或提取剖面。

3. 用 Python + GDAL 把云南 30 米 DEM 裁剪拼接成可用产品

3.1 先想清楚在哪套坐标系里完成拼接

30 米 DEM 刚解压时通常是 EPSG:4326 的经纬度格网。在这个坐标系里拼接最省事，因为所有分幅的边缘严格对齐，没有投影重采样误差。我通常先在 4326 下做一个全省 VRT 和裁剪，得到一个只覆盖云南的经纬度栅格；等真正做坡度、坡向或流域面积计算前，再一次性重投影到等积投影。这样能把“拼接”和“重投影”两件事拆开，参数改起来也少踩坑。

3.2 用 gdalbuildvrt 把 90 个分幅变成虚拟文件

不要直接把 90 个 tif 用gdal_merge.py合成一个大文件。更稳的做法是先做 VRT：

gdalbuildvrt -overwrite \ -resolution user \ -tr 0.0002777777778 0.0002777777778 \ yunnan_30m.vrt *.tif

参数说明。-resolution user配合-tr把虚拟栅格的像元尺寸钉死在 1 弧秒，防止某个分幅混入异源分辨率时，VRT 自动把整张栅格迁就到更低分辨率。-tr的值写 0.0002777777778，即 1/3600 度的浮点表达，多写几位小数能减少边缘对齐时的舍入误差。如果目录里同时有 SRTM 的.hgt文件，追加*.hgt到命令末尾；注意 bash 在没有匹配文件时会原样传字符串，先确认文件存在再加。

VRT 本身不复制像素，文件只有几百 KB；后续 gdalwarp 在输出时按需求取源像元。这意味着改重采样算法、改裁切范围都不需要重建 VRT。唯一要注意的是，VRT 对随机 I/O 敏感，90 个分幅在机械硬盘上第一次生成会明显慢，放到 SSD 或本地目录会顺利得多。前面在 2.2 已经确认过所有分幅的分辨率一致，这里-tr看起来是重复劳动，却是防御性的：万一某个分幅边缘像元有 0.0000001 度的偏差，-resolution user能让整个虚拟文件对齐到同一规则网格，接边处不会出现半像元错位。

3.3 按云南省界裁剪，政区边界外的像素留着也是噪声

VRT 生成后，下一步就是把范围裁到云南边界。用 rasterio 的 mask 函数比 gdalwarp 更直观，适合在 Python 工作流里调试：

import geopandas as gpd import rasterio from rasterio.mask import mask boundary = gpd.read_file("yunnan.shp") if boundary.crs is None: raise RuntimeError("省界缺少坐标系定义，请先补充 .prj 或重新导出") if boundary.crs.to_string() != "EPSG:4326": boundary = boundary.to_crs("EPSG:4326") with rasterio.open("yunnan_30m.vrt") as src: if src.crs.to_string() != "EPSG:4326": raise RuntimeError("DEM 不是 EPSG:4326，请先统一坐标系") shapes = [g for g in boundary.geometry] dem, out_transform = mask( src, shapes, crop=True, all_touched=False, invert=False, nodata=-9999, ) meta = src.meta.copy() meta.update( driver="GTiff", height=dem.shape[1], width=dem.shape[2], transform=out_transform, nodata=-9999, compress="deflate", tiled=True, bigtiff="IF_SAFER", ) with rasterio.open("yunnan_30m_cut.tif", "w", **meta) as dst: dst.write(dem)

这段代码先读省界 shp，并强制把省界转换到与 DEM 相同的 EPSG:4326。mask()的核心参数含义是：crop=True让输出范围收缩到多边形外接矩形；all_touched=False表示只保留像素中心落在云南界内的栅格，结果是锯齿感边界但无效值极少；如果你希望边界外保留一圈用于填洼，改True即可。invert=False指保留多边形内部，不是反向裁剪。meta.update里把压缩和瓦片存储打开，后续处理大文件时访问效率差异明显；bigtiff="IF_SAFER"会在文件超过 4GB 时自动升级。

命令行方案在脚本化交付里更常用，参数和 rasterio 是等效的：

gdalwarp -overwrite \ -cutline yunnan.shp -crop_to_cutline \ -dstnodata -9999 \ -co COMPRESS=DEFLATE -co TILED=YES -co BIGTIFF=IF_SAFER \ yunnan_30m.vrt yunnan_30m_cut.tif

-crop_to_cutline是决定成败的参数。若不写，gdalwarp 默认保留整个 VRT 范围，云南外的空白会一起被裁成无效值大背景，文件虚胖几倍。早期 GDAL 写法是-cropToCutline，如果你的 gdalwarp 版本较老，换成驼峰拼写。-dstnodata统一写 -9999，不管源数据是 -32768 还是空值，输出里只有一个无效值标记，后面所有算法只认这一个。

3.4 重投影到 Albers 等积投影，30 米才真正是 30 米

云南横跨约 9 度经度，如果坚持用 UTM，一个省会落在 47N 和 48N 两个带，全省拼接需要两次投影再接边，误差和复杂度都上来了。全省尺度的地形分析建议一次性转到 Albers 等积投影，像元用米为单位：

gdalwarp -overwrite \ -t_srs "+proj=aea +lat_1=22 +lat_2=28 +lat_0=0 +lon_0=102 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs" \ -tr 30 30 -r bilinear -dstnodata -9999 -ot Float32 \ -co COMPRESS=DEFLATE -co TILED=YES -co BIGTIFF=IF_SAFER \ yunnan_30m_cut.tif yunnan_30m_aea.tif

+lat_1=22和+lat_2=28是双标准纬线，覆盖云南的主要纬度带；+lon_0=102中央经线接近云南地理中线，能最大限度减小距离变形。-tr 30 30把输出像元固定为 30 米边长。输出用 Float32 而不是 Int16，是因为双线性插值会产生小数高程，int16 会截断成整数，斜坡细节直接丢失。这套参数做出来的 Albers 文件，坡度、坡向和面积计算都不再受纬度影响。

4. 云南 30 米 DEM 的质量验证：无效值、控制点与坐标系的隐形误差

4.1 无效值清洗：不要让 -9999 参与坡度计算

3.3 输出的文件还带着边界一圈无效值。在填洼工具里，-9999 会被当成一个极深的洼地，直接把整个云南省界外区域和边界像素都连成一个巨型坑。所以质量验证第一步不是看平均高程，而是检查统计信息：

gdalinfo -stats yunnan_30m_aea.tif

如果STATISTICS_MINIMUM等于 -9999，说明无效值仍在图像内部；如果最小值接近 50～80 米，说明裁剪和转投影已经把无效值清干净。检查后仍然要显式把全部负值视为无效值处理，因为这个高程场在云南不可能出现负值。用一次 numpy 清洗：

import rasterio import numpy as np with rasterio.open("yunnan_30m_aea.tif") as src: data = src.read(1) meta = src.meta.copy() data[data < 0] = -9999 with rasterio.open("yunnan_30m_clean.tif", "w", **meta) as dst: dst.write(data, 1)

这段判断对云南安全，对其它沿海省份不一定安全，因为海平面以下区域可能有真实负高程。写成通用脚本时应把data < 0改成data == src.nodata，并额外检查src.nodata是否为空。

4.2 用控制点比对实测高程，算 RMSE 而不是肉眼对比

质量验证最可靠的是拿已知高程控制点与 DEM 采样值对比。常见做法是每个县选 3～5 个点，覆盖河谷、城区、山顶，记录经纬度和实测高程，再写一个 Python 循环采样。如果读取的是 4.1 清洗后的 Albers 文件，需要先把经纬度转成投影坐标：

import rasterio import numpy as np from pyproj import Transformer # (lon, lat, measured_height) points = [ (100.21, 25.35, 1915.2), (101.50, 24.80, 1240.6), ] with rasterio.open("yunnan_30m_clean.tif") as src: transformer = Transformer.from_crs("EPSG:4326", src.crs, always_xy=True) diffs = [] for lon, lat, meas in points: x, y = transformer.transform(lon, lat) row, col = src.index(x, y) window = ((row, row + 1), (col, col + 1)) val = src.read(1, window=window)[0, 0] if val is None or val == src.nodata: continue diff = meas - val diffs.append(diff) print(f"({lon},{lat}) 实测 {meas:.1f} m，DEM {val:.1f} m，差 {diff:.1f} m") rmse = np.sqrt(np.mean(np.square(diffs))) print(f"RMSE = {rmse:.2f} m")

Transformer.from_crs(..., always_xy=True)保证经纬度按 (x, y) 即 (lon, lat) 顺序传入，避免方向反了；src.index(x, y)把投影坐标换算成行列号，window参数只读取单个像元，避免把整张栅格读进内存。RMSE 少于 10 米对 30 米级 DEM 属于可接受范围；如果出现 30 米以上的系统性偏差，先检查坐标系转换而不是抱怨数据源。

4.3 坐标系“隐形误差”：为什么坡度在不同投影下会差 5%

同一份 DEM 在 EPSG:4326 下直接算坡度，和在 Albers 投影下算坡度，结果可能相差 5%～15%。原因是经纬度格网横纵几何比例随纬度变化，云南北部约 29°N 处一个像元的东西向实际宽度约 27 米，南北向约 30 米，而坡度算法默认把两个方向视作等长。等高线密集的横断山区尤其敏感。

所以第 3 章先做 4326 拼接、后转 Albers 不是多此一举。它让拼接基于原始采样，让后续所有地形因子计算在一个几何均匀的空间里进行。如果你在 GIS 界面里直接用经纬度 DEM 叠加在线底图，看起来没差，但坡度、坡长和汇水面积都会被放大或缩小。

5. 用 30 米 DEM 提取地形因子：坡度、坡向、填洼与流域分析的参数选择

5.1 gdaldem 三连：坡度、坡向、山体阴影一次生成

有了 Albers 投影的 30 米格网，地形因子可以用 GDAL 自带工具生成，不需要安装重型 GIS。三行命令：

gdaldem slope yunnan_30m_clean.tif yunnan_slope_deg.tif gdaldem aspect yunnan_30m_clean.tif yunnan_aspect.tif gdaldem hillshade yunnan_30m_clean.tif yunnan_hillshade.tif -z 2.0 -az 315 -alt 45

gdaldem slope输出的是每个像元的坡度角，单位度，平坦为 0，垂直为 90。因为输入是 Albers 投影，水平单位是米，高程单位也是米，GDAL 不再需要-s缩放参数；如果你拿 4326 文件直接跑，必须加-s 111120，否则坡度会小得离谱。gdaldem aspect输出坡向角度，0 度指向正北，顺时针递增，平坦区域返回 nodata。hillshade里的-z 2.0是垂直放大倍数，云南高山峡谷用默认 1.0 会显得平，2.0 让阴影像有立体感；-az 315光线来自西北，-alt 45是太阳高度角 45 度，这是制图最常见的组合。

5.2 填洼：河网提取前必须做，但别让“真实洼地”消失

30 米 DEM 在做流向提取前，通常要先填洼。直接用原始 DEM 跑 D8 流向，任何一个比周围低 1 米的像元都会成为断流点，河网里出现大量不应该存在的孤立河段。填充后的 DEM 才能让水沿最大坡度方向一路流到出口。

常见做法是使用 WhiteboxTools 的 FillDepressions，命令行在脚本里很稳定：

whitebox_tools --run=FillDepressions \ --dem=yunnan_30m_clean.tif \ --out=yunnan_fill.tif --fix_flats=true

--fix_flats=true会把填洼后产生的大片平坦区做微地貌调整，让 D8 流向在平坦区也有明确方向。这个参数在云南坝区尤其重要，比如昆明盆地边缘的缓坡区，不处理会在河网提取时生成平行水流。

5.3 云南特殊地形下的参数修正

3×3 中值滤波是另一个常用预处理。横断山区的 30 米 DEM 在陡峭峡谷底部会出现混合像元，一个像元同时覆盖河面和崖壁，高程被平均后形成短小的锯齿。用scipy.ndimage.median_filter把单点尖刺压掉，再做坡度：

from scipy.ndimage import median_filter import rasterio with rasterio.open("yunnan_30m_clean.tif") as src: dem = src.read(1) meta = src.meta.copy() dem_f = median_filter(dem, size=3) with rasterio.open("yunnan_30m_med3.tif", "w", **meta) as dst: dst.write(dem_f, 1)

size=3的窗口只影响单像素级异常，对真正的地形轮廓几乎无改变。若做滑坡评价，我倾向于用滤波后的 DEM 算坡度；若做填挖方和剖面分析，用原始 DEM，因为中值滤波会把 3 米宽的岩墙抹平。预处理参数没有通用默认值，取决于你要算的是形态还是水文。

6. 把云南 30 米 DEM 发布成服务前，先做这 3 个验证

6.1 验证 1：裁剪后无效值占比超标等于没裁干净

服务发布前，我通常先算一下云南裁剪结果里 nodata 的占比。边界锯齿会留一圈无值像素，占比超过 1% 就说明裁剪参数选择有问题。

import rasterio with rasterio.open("yunnan_30m_cut.tif") as src: a = src.read(1) ratio = (a == src.nodata).mean() * 100 print(f"nodata 占比: {ratio:.2f}%")

如果比例高于 1%，回到第 3 章把all_touched改成False；如果低于 0.1%，说明边界是紧贴省界的干净产品，可以直接进入发布流程。

6.2 验证 2：COG 转换必须保留统计值

发布到切片服务前把 GeoTIFF 转为 Cloud Optimized GeoTIFF，能显著减少全景拖动时的加载延时，具体命令如下：

gdal_translate yunnan_30m_clean.tif yunnan_30m_cog.tif \ -of COG -co COMPRESS=DEFLATE -co BLOCKSIZE=512 \ -co OVERVIEW_RESAMPLING=NEAREST -co BIGTIFF=YES

转完后用gdalinfo -stats对比原始文件的 min/max/mean。COG 内部概览金字塔如果用了双线性重采样，均值会有轻微变化但 min/max 通常不变。若 min 从 76 变成 75 以下，说明重采样时插值到了边界外的无效值，应改OVERVIEW_RESAMPLING=NEAREST或先清掉边界再转。

6.3 验证 3：检查边界锯齿和叠底偏移

最后一步是把 COG 与云南省界的矢量叠加到同一个 GIS 窗口，缩放至 1:10000 检查道路和县界两侧的高程突变。30 米 DEM 的边界锯齿通常小于一个像元，但如果坐标参考基准不一致，云南西部高山区会出现明显的水平偏移。发现偏移后别手动平移 DEM，先把省界重新转成 EPSG:4326，再走一遍第 3 章的裁剪流程。

瓦片服务上线后，把 COG 的 mean 和 stddev 写进服务元数据，前端用mean +/- 2*stddev做拉伸范围，这样不同批次数据在浏览器里渲染出的明暗层次一致，接边处不会出现一眼能看出的亮度断层。

本文还有配套的精品资源，点击获取