简介：云南省保山市30米分辨率DEM数字高程数据包，面向GIS测绘、环境规划、工程勘察等领域的从业者与研究者，可用于地形分析、洪水淹没模拟、地质灾害评估、气候区划与城乡规划等典型任务。压缩包内共12个文件，以TIFF高程栅格、Shapefile区域边界、投影坐标文件、属性表及XML元数据为主，其中.tif为高程核心，.shp/.dbf/.prj等界定保山市范围并定义坐标系统，.tfw/.ovr辅助地理配准与快速渲染，整体体积约81.44MB。数据附带保山市行政范围shp文件，便于按研究区统一裁剪与叠加分析，已有399人学习下载。依托30米网格，数据能清晰呈现滇西高原河谷、山脊与陡坡等地貌细节，可直接导入ArcGIS、QGIS等平台提取坡度坡向、划分流域、开展淹没与地质灾害风险识别，也可结合土地利用、气象等专题数据，为科研与项目决策提供基础底图。
1. 这个"保山市DEM数字高程数据30m"的zip，先看内容再动手
拿到"云南省保山市DEM数字高程数据30m（含区域范围shp文件）.zip"，多数人的第一反应是解压、拖进ArcMap，然后直接做坡度分析。实际上，这个压缩包里的两个核心数据——30米分辨率的DEM栅格和保山市区域边界的shp文件——多半不是同一坐标系，甚至不是同一套制图基准。30m DEM一个像元对应地面900平方米，在保山这种起伏较大的滇西地形上，坡度计算会直接受投影变形影响；shp文件如果没有对应的投影文件（.prj），软件会默认当成WGS84经纬度，稍后图层就错开几公里。这篇文章顺着"解压→检查→配准→裁剪→提取地形因子→验证"讲一套能落地的工作流，适合在使用DEM和shp文件做地形分析的GIS工程师、测绘与规划从业者，也适合刚接触开源GIS的学生。
2. 检查保山市DEM与shp文件：坐标系、范围和像元深度
在动手裁剪之前，先用工具看清楚数据的底细。常见做法是用GDAL的命令行工具读栅格元数据，用Python或QGIS读shp范围。这里用gdalinfo和ogrinfo就能完成。
2.1 用gdalinfo和Python读出一个DEM的元数据
假设解压后得到dem_baoshan_30m.tif，先运行：
gdalinfo dem_baoshan_30m.tif输出中重点看Size is 2400, 1800这样的尺寸、Pixel Size = (0.00027777778, -0.00027777778)这样的像元大小，以及Coordinate System is:这一行。像元大小是0.00027777778度，约等于30米（1度约110km，0.00027778*110km≈30m）。但如果坐标系统不是经纬度而是UTM，那么像元大小就会直接显示30,30。很多从数据平台下载的30m DEM默认使用WGS84经纬度，直接计算坡度和面积时会产生约1.5~2%的变形，在保山这种北纬25度左右的区域，东西方向被拉伸的比例更明显。
接下来用Python检查shp文件范围，这里借助fiona或gdal的ogr接口。先安装依赖：
pip install fiona shapely然后跑一个脚本：
import fiona with fiona.open('baoshan.shp', 'r') as src: print(src.crs) print(src.bounds) for feature in src: geom = feature['geometry'] # 只打印第一个要素的边界 print(geom['type'], feature['properties']) break这段代码会输出shp文件的坐标系、最小外接矩形和第一个要素的属性。对比DEM的范围，确认二者是否大致重合。如果DEM的范围是98.5, 24.8, 100.5, 26.2，而shp的范围是237235, 2745678, 375345, 2987654，那很明显一个是经纬度，一个是投影坐标系，必须做配准。
2.2 shp文件不是单一文件：shp/shx/dbf/prj的作用
一个完整的shp文件至少包含四个文件。.shp存放几何坐标，.shx是空间索引，.dbf存放属性表，.prj是投影信息。很多下载站的包只给.shp和.dbf，没有.prj，这时导入QGIS或ArcMap时会被当作未知坐标系，软件通常默认当作WGS84经纬度。正确的做法是先打开.prj文件（纯文本），确认它的坐标系描述；如果没有.prj，用ogrinfo -al -so baoshan.shp查看坐标值范围也能大致判断：X在98~100、Y在24~26之间的通常是经纬度；X有6位数、Y有7位数的是高斯投影或UTM。
提示：在ArcMap中导入没有
.prj的shp文件时，会弹出坐标系警告。建议先用此方法确认原始坐标系再定义，不要直接选WGS84。
2.3 坐标系配准的决策矩阵：经纬度还是投影
30m DEM到底用经纬度还是投影坐标系，取决于后续分析。计算坡度、坡向、山体阴影这类地形因子时，GDAL算法内部会按像元大小计算梯度，如果DEM是经纬度，在北纬25度处东西向一个像元约27.5米，南北向约30.8米，差值会造成坡度和坡向的系统偏差。因此我一般会先把DEM投影到UTM 47N（EPSG:32647）或CGCS2000 / 3-degree Gauss-Kruger zone 40（EPSG:4541或相关），再做分析。shp文件如果是经纬度，同样重投影到一致的投影坐标系。
下面是一个用GDAL把DEM和shp统一到UTM 47N的例子（先处理shp）：
ogr2ogr -t_srs EPSG:32647 baoshan_utm47.shp baoshan.shp然后处理DEM：
gdalwarp -t_srs EPSG:32647 -tr 30 30 -r bilinear dem_baoshan_30m.tif dem_baoshan_utm30.tif-tr 30 30让输出像元尺寸保持30米，-r bilinear是双线性重采样，适合连续高程数据。处理后再用gdalinfo确认Pixel Size = (30, -30)。这样DEM和shp就处在了同一个坐标系，裁剪才能对得上。
| 分析任务 | 建议投影 | 理由 |
|---|---|---|
| 坡度/坡向/山体阴影 | UTM 或 高斯-克吕格 | 像元正方形，角度计算无畸变 |
| 面积/距离统计 | UTM 或 等面积投影 | 避免纬度方向的拉伸 |
| 与经纬度栅格叠加 | WGS84 经纬度 | 保持网格对齐 |
| 出图展示 | 按项目要求 | 无强约束 |
3. 用保山市shp边界裁剪30m DEM：QGIS与gdalwarp两个做法
投影统一后，下一步是用shp边界把DEM裁剪到保山市区域。注意，这里讲的不是简单地按矩形范围裁剪，而是按shp的多边形裁剪，边缘要贴合行政边界。有两种常见做法：QGIS界面操作适合临时验证，gdalwarp命令行适合批量处理多块分幅数据。
3.1 在QGIS中把shp裁剪DEM的最小操作
打开QGIS，把dem_baoshan_utm30.tif和baoshan_utm47.shp拖入图层面板。确认图层顺序是shp在上。然后从菜单选择"栅格工具→裁剪栅格"（或Processing Toolbox里搜Clip raster by mask layer）。
参数设置时注意：
- 输入图层选择DEM，掩膜图层选择shp；
Target extent需要选择当前图层范围（保山市shp）或手工输入；Output resolved或保持原分辨率勾选，避免输出像元变成大网格；- 输出格式选GeoTIFF，文件名为
dem_baoshan_clip.tif。
执行后放大查看边缘，如果裁剪后的DEM与shp边界仍有缝隙，说明两个图层的坐标系或空间基准不严格一致。常见的坑是：shp虽然带.prj，但坐标系是CGCS2000，而DEM是WGS84。两者在1~2米的差距在30米分辨率下通常不明显，但在边界上会露出锯齿状空白。解决办法是用gdalwarp强制将DEM重投影到shp的坐标系。
3.2 gdalwarp裁剪命令与关键参数
命令行方式更适合写成脚本复用。核心命令：
gdalwarp -t_srs EPSG:32647 -cutline baoshan_utm47.shp -crop_to_cutline -dstnodata -32768 -of GTiff dem_baoshan_utm30.tif dem_baoshan_final.tif说明：-cutline指定shp矢量路径；-crop_to_cutline表示让输出范围贴合多边形，而不是使用shp的外包矩形；-dstnodata -32768把裁剪后多边形外部的区域设为NoData，避免黑色0值参与后面的坡度计算；-of GTiff输出GeoTIFF。这里有一个容易犯的错：只加-cutline不加-crop_to_cutline，命令会以shp的外包矩形裁剪，保山边界外的像元会保留为0值干扰后续统计。
裁剪完成后，用gdalinfo dem_baoshan_final.tif | grep -E "Size|NoData|Pixel"验证。正常时，NoData Value=-32768，Size的宽高比例与shp边界的外包矩形一致。
3.3 不同裁剪工具的优劣势对比
| 工具 | 优势 | 劣势 | 适用场景 |
|---|---|---|---|
| QGIS 裁剪栅格工具 | 可视化，简单 | 每次一座，无法批量；参数藏在GUI里 | 一次性任务 |
| gdalwarp | 可批量、可精确控制重采样和投影 | 需要记忆参数 | 批量处理或脚本化 |
| ArcGIS Extract by Mask | 与ArcGIS其他工具无缝衔接 | 授权限制，命令行不易传参 | 已有Desktop环境的用户 |
如果要批量裁剪多个shp（比如把保山市每个县区的shp都裁一份），可以用一个bash循环：
for shp in /path/to/counties/*.shp; do name=$(basename "$shp" .shp) gdalwarp -cutline "$shp" -crop_to_cutline \ -tr 30 30 output_${name}.tif dem_baoshan_utm30.tif done注意-tr 30 30在批量情况下必须加，否则裁剪后的像元尺寸会根据输出范围自动调整，可能与原来的30m不一致。
4. 从裁剪后的DEM提取坡度、坡向与山体阴影：gdaldem参数设置
DEM最主要的用途是提取地形因子。常见的是用gdaldem命令，它可以根据DEM生成坡度、坡向、山体阴影和彩色高程图。参数选择直接影响结果，下面逐个讲清楚。
4.1 坡度是度数还是百分比？这里按实际需求选
坡度常见有两种输出：度数和百分比。度数适合水土保持和生态学，百分比适合工程与道路设计。比如保山市算平均坡度，度数直观；但在评估坡地建设时，百分比坡度（如15%以下）更常被引用。gdaldem的-p选项输出百分比，不写该选项则默认输出度数。
命令示例：
gdaldem slope dem_baoshan_final.tif baoshan_slope_deg.tif gdaldem slope -p dem_baoshan_final.tif baoshan_slope_percent.tif坡向输出是0到360的角度，0表示北、90表示东，方向代表水流去向。一般在山地分析中很少直接使用，需重分类为8个方位。山体阴影命令按光源方向生成灰度栅格，适合做底图：
gdaldem hillshade -az 315 -alt 45 -z 2 dem_baoshan_final.tif baoshan_hs.tif参数-az 315表示光源方位角315度（西北），-alt 45是太阳高度角45度，-z 2是高程放大倍数。保山位于横断山脉余脉，地形起伏大，-z取值可以降低到1.5，不然阴影会过深。
提示：如果DEM存在很多小沟壑，坡向输出会显得非常噪。先对DEM做一次中值滤波或焦点统计平滑，再算坡向，结果更干净。
4.2 用gdaldem生成三套地形因子
一次完整的分析通常同时生成坡度、坡向、山体阴影。可以放在脚本里执行：
# 生成坡度（度） gdaldem slope dem_baoshan_final.tif baoshan_slope.tif # 生成坡向 gdaldem aspect dem_baoshan_final.tif baoshan_aspect.tif # 生成山体阴影 gdaldem hillshade -az 315 -alt 45 -z 1.8 dem_baoshan_final.tif baoshan_hillshade.tif三个输出GeoTIFF的坐标系和范围与输入DEM一致。接着可以做一个检查：用gdalinfo -stats baoshan_slope.tif | grep "Minimum="查看坡度的最小值。如果最小值为负数，说明DEM存在边缘NoData被计算为0；如果最大值超过90，说明DEM有异常高程值。
4.3 叠加shp和DEM做可视化出图
生成的坡度图层要与shp边界一起展示，以说明数据来源。QGIS中，把baoshan_hillshade.tif放到底层，baoshan_slope.tif用半透明渐变叠加上去，然后设置shp边界的透明度为70%。合适的比例是山体阴影用灰度，坡度用彩虹色带，shp边界用白色线宽0.7。这样在保山地图上能一眼看出坝区与山区的过渡。
如果要把叠加结果保存为图片，用QGIS的"导出地图"功能，注意DPI至少150，避免截图模糊。也可以直接用Python脚本从GeoTIFF合成PNG，但那需要额外的matplotlib代码，对一般工程来说不划算。
5. 验证DEM数据可用性的三个快速方法，以及和不完美DEM共存
拿到别人的DEM数据，最怕里面有空洞、负值或接边错位。这里提供三个验证方法，最后一个技巧能帮你在不重新下载的情况下把数据救回来。
5.1 用直方图和基础统计识别DEM空洞
先运行：
gdalinfo -hist -stats dem_baoshan_final.tif查看Minimum、Maximum和直方图。高程值范围如果在保山地区，一般不会为负（保山市最低海拔约500米，最高在4500米左右）。如果Minimum是-32768，说明NoData没被正确读取。这时候用gdal_calc.py把NoData替换为DEM有效范围内的最小值或者在坡度和山体阴影计算时忽略它。更稳妥的验证是直接看空洞分布：用Python读入数组，统计0值或NoData像元所占比例：
from osgeo import gdal import numpy as np ds = gdal.Open('dem_baoshan_final.tif') band = ds.GetRasterBand(1) nodata = band.GetNoDataValue() arr = band.ReadAsArray() if nodata is not None: invalid = (arr == nodata).sum() else: invalid = (arr <= -9999).sum() print('invalid pixels:', invalid, 'of', arr.size, 'ratio:', invalid / arr.size)如果ratio大于0.1，说明裁剪时保留了大量边界外区域，建议重新用-crop_to_cutline裁剪。
5.2 和ALOS 12.5米等高分辨率DEM比较
在关键区域，可以用ALOS 12.5米DEM做交叉验证。ALOS 12.5米是JAXA提供的全球公开数据，比30m更精细。下载范围覆盖保山市的分块后，先重采样到30m，再与现有DEM做差值：
gdalwarp -tr 30 30 alos_baoshan_12m.tif alos_baoshan_30m.tif gdal_calc.py -A dem_baoshan_final.tif -B alos_baoshan_30m.tif --outfile=diff.tif --calc="A-B"运行后用gdalinfo -stats diff.tif查看差值均值和标准差。正常情况下，30m DEM与ALOS之间的均方差在5米以内，如果标准差超过15米，说明数据存在系统性误差或投影问题。
5.3 一个小技巧：用shp投影同步DEM，避免重复折腾
处理多个DEM分块时，最烦的是每个tif的投影都不一致。一个省事的小技巧是把shp的投影文件直接复制给DEM（仅当shp的坐标系是准确投影坐标系时）。在Linux或Windows的GDAL环境里，用gdalsrsinfo导出shp的WKT，再用gdal_edit.py写入DEM：
gdalsrsinfo -o wkt baoshan_utm47.shp > baoshan.prj gdal_edit.py -a_srs baoshan.prj dem_baoshan_final.tifgdal_edit.py只修改元数据，不重采样，速度快且不损失原有像元值。它适用于DEM本身坐标位置正确，仅仅缺失.prj或投影信息错误的情况。如果实际坐标已经有偏移，gdal_edit.py救不了，还需要用gdal_translate -a_ullr重新设定四点坐标。
本文还有配套的精品资源，点击获取
