简介：云南省普洱市30米分辨率数字高程模型（DEM）数据包，面向GIS开发、地理分析与规划学习者，提供精细地形栅格及配套行政边界矢量文件。DEM通过等间隔海拔值描述地表形态，本数据分辨率30米，可支撑中小尺度地形研究。压缩包共12个文件，核心包括TIFF高程栅格（普洱市DEM.tif及其金字塔.ovr、世界文件.tfw、辅助元数据.xml）和Shapefile边界矢量（普洱市范围.shp及其属性.dbf、投影.prj、索引.shx/.sbn/.sbx），整体约211.48MB，格式规范完整。数据可直接导入QGIS、ArcGIS等平台，用于地形渲染、坡度坡向提取、淹没模拟或城市选址等场景。目前已有526人学习使用。该数据将高程信息与普洱市行政边界无缝整合，省去自行配准和裁剪步骤；投影文件和元数据便于坐标系确认与数据溯源，适合开展区域地理研究、教学演示及国土空间辅助分析。
1. 这个DEM压缩包不是“地图”，是能算账的地形底图
下载到本地的是一个叫“云南省普洱市DEM数字高程数据30m（含区域范围shp文件）.zip”的压缩包，里面是一张30米分辨率的普洱市数字高程模型（DEM），外加一个界定普洱市范围的shp面文件。很多人第一反应是把它当“地形图”看，实际上，DEM和shp组合的价值在于“分析”：拿到它，你可以把普洱市的地形从周边切出来，算坡度、坡向、山体阴影、地形起伏度，给道路选线、汇水区划分、光伏选址或者灾害评估做底图数据。适合谁用？接到区域性制图任务、做GIS课程设计，或者正准备开展地形相关分析项目的从业者。下面我会按处理这类数据最常见、最可靠的路径，把整个流程从解压到出成果走一遍，重点说清楚哪些参数不改就会翻车。
2. 先拆包看透这两样东西：30m高程的存储方式和shp的边界角色
拿到zip后不要急着拖进ArcGIS，先把它当成一个“待验收的数据包”来处理。理解它的内部结构，能省掉后面至少一半的排查时间。这章我会拆开讲DEM文件的格式本质、shp文件的真实组成，以及最关键的坐标基准问题。
2.1 30m到底指什么：不是精度，是栅格像元尺寸
标题里的“30m”指的是DEM栅格数据的空间分辨率——每个像元在地面上对应的实际边长约为30米×30米。这是数字高程模型的行业惯例表达，不代表高程测量精度达到30米。实际上，这类开放获取的30米DEM产品，高程绝对误差通常在±5米到±15米之间，具体看原始数据源和处理工艺。
一个像元代表一个地面范围，像元内的高程值是一个统计或采样结果。栅格文件里存储的就是一个二维矩阵，每个格子存一个高程值。你在ArcGIS或QGIS里把它拉伸成彩色图像来显示，那只是“可视化”，不是数据本身。理解这点后，后面所有派生分析才好理解：坡度、坡向、山体阴影，都是在这个矩阵上做邻域运算。
这类数据常见的包装格式有GeoTIFF（.tif）、IMG（Erdas Imagine格式）、以及部分国内的Grid格式。标题里既然叫“DEM数字高程数据”，大概率是GeoTIFF封装，因为它能内嵌坐标系信息和无效值标识，跨平台兼容性最好。打开一看如果是.img也不用慌，QGIS和ArcGIS都原生支持。注意一点：30m分辨率在经纬度坐标下表现为近似0.00027778度（即1/3600度，也叫1弧秒），在投影坐标下才是整齐的30米整数格网。这两个坐标系下的“30m”文件，实际处理流程差异不大，但涉及距离计算和面积计算时，必须统一到投影坐标系。
2.2 shp不是一个文件，是一组文件：少一个都打不开
“区域范围shp文件”这个说法在业内很常见，但严格来说shp是一组文件的统称。一个完整的Esri Shapefile至少包含三件套：主文件（.shp，存储几何）、索引文件（.shx，存储几何索引）、属性表（.dbf，存储属性字段）。如果带了坐标系定义，还会有.prj文件；在Windows系统下通常还有.cpg（字符编码）和.sbn/.sbx等可选文件。
我们拿到的是zip压缩包，解压后要检查是否同时包含上述几个文件。很多人只把.shp拖进GIS，结果图层不显示图形，其实是因为.shx或.dbf缺失。这里有个判断技巧：如果压缩包里的.shp有几十MB以上，它多半是完整的；如果只有几KB，要么只是边界简化为几条线，要么就是文件缺胳膊少腿。
打开shp后，用属性表看看有多少要素、有没有字段。区域范围文件通常是面要素，也就是一个多边形，表示普洱市的行政边界或研究区域边界。它是矢量数据，不是栅格。它的作用是“火柴盒”——决定DEM裁剪范围，而不是给你提供地形信息。另外，检查一下shp的几何类型：有“Polygon”（面）和“MultiPolygon”（多面）之分，如果普洱市由多个地块组成，shp里会有多个要素。后续做掩膜裁剪时，这种多要素边界用法和单要素没有本质区别，但导出时要留意坐标系是否统一。
2.3 坐标基准是第一关：WGS84还是CGCS2000，UTM分带怎么选
这是整个流程里最容易翻车的地方。普洱市位于云南省南部，地理范围大约在东经100度到102.3度之间，北纬22度到24.5度左右。这个经度范围恰好跨越了UTM 47N和48N两个分带。如果DEM和shp使用不同的坐标基准或投影带，叠加后会出现肉眼可见的整体偏移，轻则几百米，重则几公里。
先做个“体检”：在QGIS中分别加载DEM和shp，右键图层属性查看坐标系信息。常用的30米DEM产品，比如SRTM V3和ASTER GDEM V3，通常采用WGS84椭球基准，以经纬度坐标存储；国内发布的数据则可能采用CGCS2000基准。两者在高程基准上差异不大，在平面位置上差异小则几十厘米、大则几米，对30米分辨率分析来说通常可接受，但区域网格系统不同，还是要先归一化。
如果shp自带的.prj文件描述是“GCS_WGS_1984”，而DEM也是WGS84经纬度，那么两者可以直接叠加，裁剪前不需要做投影。如果后面要算面积、坡度坡向、距离，就要把数据统一投影到UTM。普洱市的中央经线区域在100度到102度之间，更靠近99度或105度的中央经线都行，我一般推荐用WGS 84 / UTM zone 47N（EPSG:32647）作为默认工作坐标系，因为普洱大部分区域落在47N内。若你的研究区横跨两带，那就要整景处理并考虑分带，或者改用自定义中央经线的投影。
提示：先记录原始坐标系，再决定是否投影。裁剪和投影这两个操作不要同时做，否则一旦结果不对，你分辨不清是哪一步出了问题。
2.4 解压后立刻执行的“四查”清单
有些坑在一开始就能避掉。我处理这类数据包的固定习惯是，解压后不急于做分析，先按下面四个检查走一遍。检查项目都在十分钟内能完成，但对后续流程是“后悔药”。
| 检查项 | 操作 | 期望结果 |
|---|---|---|
| 文件完整度 | 列出目录，检查.shp/.shx/.dbf/.prj是否齐全；DEM的.tif能否预览 | 无缺失文件，预览能显示地形起伏 |
| 坐标系属性 | QGIS属性信息查看COORDINATE REFERENCE SYSTEM | 两个数据都有明确坐标系，且是同一基准 |
| 无效值（nodata） | 用栅格统计查看最小值和最大值比例 | 无效值数量占比合理，且不是大面积黑洞 |
| 数据格式 | 查看像元深度（8bit/16bit/32bit float）和压缩类型 | 16bit整数或32bit浮点都可用，压缩类型不阻塞读写 |
重点说下无效值。DEM产品在填补数据空洞时，常用-9999、-32767或0来表示“无数据”。如果你不处理这些值，裁剪后统计平均值会严重失真，还会在坡度计算时产生诡异的大数值边界。在GDAL或ArcGIS里，读取时就要识别nodata设置。很多“处理结果全是黑色空洞”的求助帖，根源就是没有处理nodata值。
3. 把shp和DEM叠起来：用QGIS做行政区范围裁剪
这章进入动手环节。用QGIS配合命令行工具，把原始DEM裁剪成普洱市边界范围内的数据，并讲清楚“掩膜裁剪”和“范围裁剪”这两种容易混淆的操作到底有何不同。
3.1 加载图层并做最基础的“空间对齐”检查
在QGIS中，把zip解压后的.tif和.shp直接拖入图层区。如果DEM是大范围影像（比如覆盖全省），加载后屏幕看到的是一整块灰度或彩色的地形；shp会以线框形式叠加在上面。这时先用“识别”工具点击shp边界和DEM某个同名地物点，看两者位置是否大致吻合。普洱市的地形特征很明显，无量山、哀牢山和澜沧江，shp边界应该与DEM上的河谷走向有对应关系。
如果发现shp边界完全不在DEM可显示的区域内，基本可以判定坐标系不一致或投影带偏差。Zoom to layer之后，数字显示范围完全对不上，那就先放弃视觉上的“看起来一样”，回看第2.3节做坐标系归一化。日常经验中，这里用到的最多问题是“shp显示出来了但裁剪结果是空白”，原因大多是shp本身偏移，而不是裁剪工具有问题。
这一步完成后，按Ctrl+Shift+C复制当前画布范围做记录，方便后续验证裁剪结果的空间范围是否正确。也可以把当前范围的左上角/右下角坐标存下来，作为裁剪结果“是否切对了”的客观凭证。
3.2 裁剪栅格 vs 掩膜裁剪：精确度完全不同
QGIS菜单栏Raster菜单下的“Extraction”包含两个常用工具：“Clip Raster by Extent”和“Clip Raster by Mask Layer”。新手最容易用错的就是它们。
“Clip Raster by Extent”是按外围矩形范围裁剪，不对shp的真实边界做任何处理。也就是说，即使shp是圆形的、不规则形状的，最终裁剪出的DEM仍是矩形区域。它在处理大范围文件时速度快，适合先粗裁到普洱市附近，再进一步细裁。
“Clip Raster by Mask Layer”才是真正按shp面的形状进行裁剪。原理是先把shp矢量栅格化为一个0和1的掩膜，落在shp内的像元保留，shp外设为nodata，最后边缘用锐化算法形成贴合边界的栅格。这个操作对计算资源要求更高，但结果才是“正正好好切出一块普洱市”。
实际操作参数值得注意：在Clip Raster by Mask Layer对话框里，有个“Crop the extent of the raster to the extent of the mask layer”选项，勾选后能避免输出一大片空值区域；缺点是处理边缘像元时会以原始像元位置为准，可能产生半像元偏移。如果后续要做像元对齐分析，建议在最后输出时用“Align Rasters”工具重采样一次。
3.3 用Python（rasterio）验证裁剪结果：代码贴出来照着跑
图形界面适合单次操作，但“验证结果”和“批量处理”最好用脚本。下面这段Python代码用rasterio做一次掩膜裁剪验证，顺便打印裁剪前后的范围、像元数量和数据有效性：
import rasterio from rasterio.mask import mask from shapely.geometry import shape import geopandas as gpd import matplotlib.pyplot as plt # 第一步：读入shp并提取几何对象 gdf = gpd.read_file("puer_admin.shp") # 如果shp是多个要素，先合并为一个几何体 geometry = gdf.union_all() print("shp边界范围:", gdf.total_bounds) # 第二步：读入DEM并用shp裁剪 with rasterio.open("puer_dem_30m.tif") as src: print("原始坐标系:", src.crs) print("原始范围:", src.bounds) print("原始nodata:", src.nodata) out_image, out_transform = mask( src, [geometry], crop=True, # crop=True时输出范围紧贴掩膜边界 nodata=-9999, # 保持和原始一致的nodata设置 filled=True # 裁剪区域外的像元填为nodata ) # 检查裁剪后有多少有效像元 valid_pixels = (out_image[0] != -9999).sum() print("裁剪后有效像元数量:", valid_pixels) print("裁剪后范围:", rasterio.transform.array_bounds(out_image.shape[1], out_image.shape[0], out_transform)) # 第三步：简单可视化，确认边界贴合 plt.imshow(out_image[0], cmap="terrain") plt.colorbar(label="Elevation (m)") plt.title("Clipped DEM for Pu'er") plt.show() # 第四步：把裁剪结果写回新的GeoTIFF with rasterio.open( "puer_dem_clipped.tif", "w", driver="GTiff", height=out_image.shape[1], width=out_image.shape[2], count=1, dtype=out_image.dtype, crs=src.crs, transform=out_transform, nodata=-9999, ) as dst: dst.write(out_image)上面这段代码分成四步：读shp合并几何，读DEM按掩膜裁剪，打印裁剪后参数，写回新文件。初学者容易漏掉“union_all”：如果shp包含多个面要素，直接传给mask会得到多个输出，后面写文件时还得逐带处理。这里先合并成一个几何，输出自然就是一个单波段文件。
参数说明：crop=True是最关键的设置，不加它输出范围会和原始DEM一样大，只是边界外的像元被置为nodata，文件依然巨大；filled=True负责把边界外的像元填成nodata，否则二维数组里区域外是随机数或空洞。最后一个可选参数“double_precision”，如果你的shp比较精细，可以调高栅格化精度，但处理时间会增加。
3.4 裁剪时要不要顺带重采样：能不动就不动
裁剪过程中，对话框中通常有“Resolution”或“Output size”设置。很多人习惯顺手把分辨率改大，比如从30米改成10米，以为能获得更精细的地形。这是一个很大的认知误区：DEM本身是30米分辨率数据，超采样并不能增加信息量，只是把每个像元复制成多个小像元，输出文件变大了几倍，后续计算速度还慢。
相反，如果把分辨率调低（比如聚合为60米），就要明确知道你已经做了重采样和均值化，这会改变坡度、坡向计算结果。除非在后处理中需要与60米分辨率数据集对齐，否则我建议裁剪时保持原始分辨率不变。裁剪本身不是数据改变的时机。
4. 30m DEM能派上的真实用场：坡度、坡向、阴影和河网
从DEM提取派生地形因子是绝大多数项目核心需求。普洱市多山，坡度分析对于道路规划、竹林种植区划和地质灾害隐患识别都有价值。这章给出几个最常用参数的提取方法和参数边界。
4.1 算坡度：度（Degree）还是百分率（Percent）？
坡度是最常用高程派生参数。QGIS的Raster菜单→Analysis→Slope工具，底层算法是用Horn方法（也叫格网邻域法），在每个像元周围3×3窗口内拟合平面。无论界面语言怎么显示，关键参数是输出单位：Degree表示水平面与斜面夹角的角度，0到90度；Percent表示垂直落差除以水平距离的百分比，平地为0，45度坡等于100%。
普洱市地形切割深，坡度差异大。如果你做对外的坡耕地统计，我一般推荐用Degree输出，语义直观；如果你做工程上的土方平衡或勘探设计，用Percent更贴近工程习惯。部分工具默认输出是Percent，不留意就会拿到一个最大上千的图，实际坡度却远没有那么恐怖。
额外要提的是“Slope与高程分辨率的关系”：30米分辨率的DEM算出的坡度值系统偏小，因为大尺度栅格会平滑局部微地形。手里只有30米数据时，计算结果用于区域尺度规划是可靠的，但不要把它当成精确到具体地块的坡度量算。
4.2 等高线生成：间距怎么选，平滑用什么方法
从DEM提取等高线，多数人用r.contour（GRASS）或ArcGIS的Contour工具。间距选多少，没有标准答案，但有选择逻辑：做全普洱市的宏观地形图，20米或50米等高距不至于线条密集成一团黑；做县域地形图，10米可读性尚可；如果针对一个小山谷地块做详细分析，5米等高距才有意义。
参数上还需要指定“Base contour”和“Equidistance”。以QGIS的Contour工具为例，最终结果会自动投影，但你要记住“等高线”是矢量线数据，由栅格算法追踪生成。算法在每个像元内判断高程是否穿过给定值，并用折线连接。如果DEM有nodata空洞，在空洞边缘等高线会突然断开，生成一片扇状的短线。处理办法是先对DEM做空洞填充，或者把nodata区域用“0”的高程值作为mask剔除掉。
平滑是后期一个容易被滥用的环节。如果你的等高线锯齿感很强，可以用GRASS工具r.contour.step配合v.generalize做平滑。注意：过度平滑会让等高线失去真实感，甚至跨越山脊线画出违背地形的形状线。实操时我一般只做“容忍度1个像元”的平滑，防止地形细节被抹掉。
4.3 山体阴影（Hillshade）：两个参数最容易影响出图观感
山体阴影图不是地形分析的核心成果，但是做PPT汇报和制图底图时，它决定整体观感。QGIS的Hillshade工具或GDAL hillshade，参数有“Azimuth”（太阳方位角）和“Altitude”（太阳高度角）。默认值通常为315度和45度，但这组参数最适合展示中低纬度地区的宏观地形。
普洱市地形起伏大，使用默认值时，东北—西南走向的山脊线会特别黑，而南坡特别亮，人眼注意力容易被高亮区域带走。如果想做平面底图，建议把Altitude改成35度左右，Azimuth改成320度，阴影层次更柔和；如果要突出沟谷水系的走向，把Azimuth设置成与主要河流垂直的方向，这样河谷会产生明显阴影，视觉上河流脉络清晰。
在GDAL中，Hillshade还有“ZFactor”参数，这是高程夸大系数。30米分辨率的DEM在1200dpi打印时，起伏感可能不足，适当用ZFactor=1.5或2，能让山地显得更有立体感，但也会放大噪声。定量分析项目，不要加ZFactor，保持1.0。
4.4 地形湿度指数与河网提取：填洼（Fill Sinks）参数别乱改
在普洱这类多山区域做水文分析，基于DEM提取水流方向是常见需求。标准流程是：填洼（Fill Sinks）→流向（D8）→累积流量→河网提取。填洼这个步骤争议很大。经典算法把地形中的所有洼地填平，但真实地形中确实存在封闭洼地，盲目填平会改变真实水文路径。
QGIS中填洼工具的参数“Maximum elevation difference to fill”，默认值0表示填平所有负地形；当你用30米DEM时，我通常限制最大填洼深度。普洱市岩溶地貌区，地下河发育，地表多漏斗状洼地，全填的话会把岩溶地貌的水文特征抹掉。保守做法是填洼深度限制在DEM高程精度范围内，比如30米DEM填洼限值设为10米，只清除噪声级洼地，保留真实地形。
提取河网的“阈值”参数也需要经验。累计流量大于阈值的像元被认为是河道，阈值越小，河网越密集。对30米DEM，我一般在面积约1万平方千米区域用1000至5000像元的阈值，得到的河网密实度与1:25万地形图基本吻合。不要随便用“默认值100”，那会得到几乎每个像元都是河道的荒谬结果。
5. DEM和shp联合使用中的常见问题排查：五个真实踩坑记录
由于DEM和shp来自不同生产流程，甚至存在多年时间差，合并使用时的坑非常集中。下面五条是我反复遇到的典型问题，每条都按“现象→原因→解决”的记录方式列出，可以直接对号入座。
5.1 现象：裁剪结果大面积空白，只有边界处有一圈数据
原因：shp的坐标系和DEM不一致。常见于shp是CGCS2000 / 3-degree Gauss-Kruger投影，DEM是WGS84经纬度。叠加后边界偏移几十到几百千米，掩膜全部覆盖在高程为nodata的区域上。
解决：在QGIS中先确认shp的CRS信息。右键shp图层，查看Source的CRS描述。如果与DEM不同，用“Vector→Data Management→Reproject Layer”将shp转成与DEM一致的坐标基准。转换后重新检查总范围，再执行裁剪。
注意：Reproject Layer时勾选“Selected features only”会只转当前选中要素，容易漏掉多要素边界。除非确认研究对象只有一个要素，否则保持不选中。
5.2 现象：DEM加载后整个画面几乎纯黑或一片灰，看不出地形
原因：界面自动拉伸策略不合适。如果数据是16bit表示的DEM，值域可能从-100到4000米；显示模块按最低最高值线性拉伸，导致大范围中段高程（例如900到1500米）被压缩成接近黑色。
解决：在QGIS的图层样式里，把渲染类型改为“单波段假彩色”或“悬崖峭壁”，并调整色带区间；最好的方式是打开“Min / Max Value Settings”选择“Min max”，并勾选“Mean ± standard deviation × 2”。ArcGIS里同理用拉伸类型“Standard Deviations”，n值设为2。不要以为数据坏了，通常只是显示问题。
5.3 现象：shp文件解压后拖进软件，提示“无效数据源”或只显示坐标表
原因：shp依赖的.shx、.dbf、.prj文件缺失或路径中文导致读取失败。Windows环境下的中文路径或特殊字符（比如括号、空格）有时也让底层GDAL读取异常。
解决：把解压出来的所有文件放在纯英文路径下（例如D:\gis_data\puer），并且确保目录下没有与.shp同名但扩展名零散的其他版本。如果提示“缺少.shx”，用QGIS的数据源管理器手动选择.shp文件（而非整个文件夹），有时能强制重新生成索引。实在缺失.shx，可用ogr2ogr重写一份完整shp：
ogr2ogr -overwrite puer_fixed.shp puer_broken.shp这条命令会读取破旧文件并重建.shx/.dbf/.prj。只要.shp主体没损坏，一般能抢救过来。
5.4 现象：裁剪后的DEM边缘出现白色细线/锯齿状空洞
原因：矢量边界与栅格像元边缘不对齐，尤其当shp精度高于栅格分辨率时，掩膜算法切割出的边界像元会有一半是nodata，显示时形成一圈齿状白线。
解决：这是栅格数据的固有属性，不是错误。后续使用中不要直接拿这种边界数据叠加其他栅格，更不要做邻域运算。解决办法是做一个“Buffer”（距离为1个像元宽度的缓冲区）后再裁剪，把边缘部分包含进来；或者在制图输出时对nodata区域设置透明色，视觉边缘便干净了。最推荐的方式是后续分析时使用“边界内缩一个像元”的掩膜分析，这能避免边缘像元的不完整高程参与坡度计算。
5.5 现象：DEM范围统计出负高程或异常高程，比如-3000多米
原因：一是无效值（nodata）没被标识，二是数据源本身包含SRTM水域掩膜值（例如水域高程为0被记录为负数）。精度较低的数据还会在陡峭峡谷或雷达阴影区域产生假低值。
解决：先用GDAL或rasterio检查数据统计信息，确认最小值的物理合理性。普洱市最低海拔大约为300米（澜沧江河谷），如果统计结果低于这个值，优先怀疑nodata。用如下命令修复nodata标记：
gdalwarp -dstnodata -9999 -of GTiff puer_original.tif puer_fixed_nodata.tif再重新计算统计信息，如果仍有负值，可用“栅格计算器”把所有小于等于0的像元赋为同一高程或设为nodata。这里要特别强调：不要小看这个步骤，直接用带假高坡的DEM算坡度，会在河谷沿线生出一串虚假的陡坡值。
6. 把这一套流程固化成批处理脚本：从单次操作到可复用成果
当你验证过单次裁剪、单次坡度提取的结果后，就该把它做成脚本，因为普洱市这种级别的区域分析，往往会涉及多块图幅、多个指标或者多次参数试验。这章分享我常用的批处理思路，以及如何把输出文件整理得井井有条。
6.1 批处理核心脚本：一个Python函数解决裁剪和派生计算
下面代码做的是一次性输入DEM与shp，输出裁剪后的DEM、坡度图和高程统计表。适合作为模板改造，直接拷贝到Jupyter或命令行运行：
import rasterio import rasterio.mask import numpy as np import geopandas as gpd from rasterio.warp import reproject, Resampling def process_dem_with_shp(dem_path, shp_path, output_dir): # 读取shp并转成统一坐标系 gdf = gpd.read_file(shp_path) with rasterio.open(dem_path) as src: # 如果shp的坐标系与DEM不一致，先重投影到DEM的CRS if gdf.crs != src.crs: gdf = gdf.to_crs(src.crs) geom = gdf.union_all() # 裁剪 out_img, out_transform = rasterio.mask.mask( src, [geom], crop=True, nodata=-9999 ) meta = src.meta.copy() meta.update({ "driver": "GTiff", "height": out_img.shape[1], "width": out_img.shape[2], "transform": out_transform, "nodata": -9999 }) # 写裁剪后的DEM dem_out = f"{output_dir}/dem_clipped.tif" with rasterio.open(dem_out, "w", **meta) as dst: dst.write(out_img) # 同时输出一份高程统计csv valid = out_img[0][out_img[0] != -9999] if valid.size > 0: stats = { "min": float(valid.min()), "max": float(valid.max()), "mean": float(valid.mean()), "std": float(valid.std()) } import csv with open(f"{output_dir}/elevation_stats.csv", "w") as f: w = csv.DictWriter(f, fieldnames=stats.keys()) w.writeheader() w.writerow(stats) print(f"处理完成: {dem_out}") return dem_out # 调用示例 process_dem_with_shp("puer_dem_30m.tif", "puer_boundary.shp", "./output")这段代码把三个最常用的动作封装成一个函数：读取shp并做坐标对齐、按掩膜裁剪、输出高程统计。核心优势在于，如果后续你再拿到相邻区域的DEM，只要换路径就能跑出一致的格式结果，避免不同批次处理的手感差异。
参数说明：to_crs逻辑上先比较再转换，能减少不必要的重复采样。另外stats输出CSV对后续填报告或对比数据非常有用，不需要再重新打开GIS软件去查属性。如果你的输出目录不存在，建议在函数里加上os.makedirs，否则会因路径错误中断任务。
6.2 批量处理多幅DEM时的文件命名与中间目录设计
如果你要处理的是整个普洱市分幅的多块DEM，建议按“原始数据”和“处理结果”两个目录来管理。原始目录里保留官方命名；处理结果目录用统一后缀。例如，原始文件“puer_dem_30m.tif”处理完输出“puer_dem_clipped.tif”，坡度输出“puer_dem_slope.tif”。这样你在项目交接给同事时，不需要额外说明书就能看出文件层级。
批处理时，脚本里的输出路径不要用中文命名。这在Linux服务器上尤其实用，在Windows上也能减少编码出错概率。QGIS自带的Batch Processing界面也可以完成同样的任务，但对参数控制的灵活性不如脚本；当你的参数超过三组时，用脚本更靠谱。
6.3 验证成果是否满足出图要求：三张图叠一个报告
最后一步，把裁剪的DEM、坡度图和shp边界同时加载进QGIS，叠加卫星影像或在线地形图，检查三点：边界是否贴实、坡度极值是否位于实际陡崖区域、阴影图有没有常识性错误（比如大面积黑色出现在阳坡）。这三关都过了，这批数据就可以放心交出去。
我个人的习惯是不管多忙，都会在交付前导出一次高原程区域分布图，打印出来对着纸质地形图扫一遍。普洱市的山地地形特征在纸质地形图上的表现是硬指标，没有异常突变一般问题不大。这套“先体检，再裁剪，后验证”的流程，我用了很多年，踩过的坑都在前几章列全了。DEM和shp的组合处理是一门实践手艺，经验比参数重要，但你只要把坐标基准、nodata和重采样这三道关把好，大多数项目都能顺利跑通。希望帮到你。
本文还有配套的精品资源，点击获取
