30米DEM数据处理全流程：裁剪、坡度坡向与填洼实战

简介福建省南平市三十米分辨率数字高程模型数据包内含区域范围矢量文件面向地理信息科学学习者、城市规划与环境评估从业者适用于地势分析、环境规划、城市设计、灾害评估等场景。三十米精度可支撑中小尺度地形研究覆盖南平市全域及周边便于区域综合分析。包体共十二个文件以高程栅格和范围矢量为主辅以投影坐标、属性数据、空间索引、地理配准、元数据等配套文件整体约九十七兆字节结构完整可直接导入主流地理信息系统平台使用。已有二百九十八人浏览学习适合作为教学实训或科研项目的真实数据源。依托该数据集使用者可练习栅格配准、投影转换、地形因子提取等操作也能借助范围矢量完成区域裁剪统计深入理解数字高程模型数据组织与矢量文件格式的协作关系。1. 拿到“福建省南平市 30m DEM”压缩包先别急着看地形那份“福建省南平市DEM数字高程数据30m含区域范围shp文件.zip”解压后基本就是两类东西一个 30m 分辨率的 DEM 栅格和一个对应区域的边界 shp。很多人第一反应是赶紧打开 DEM 看山形我的习惯是先开 shp、读元数据、核对坐标系。因为 30m DEM 在分发过程中经常被重投影、被截成外接矩形包内自带 shp 就是给分析范围兜底用的。对大多数项目来说这个包能解决两件事地形底图来源和区域边界统一。像汇水区划分、坡度分级、房址初筛、通信覆盖模拟这类市域尺度的分析30m 分辨率是性价比最高的起点。但如果你打算拿它做宅基地或地块内部的高程推算至少要再叠加更高精度控制点否则坡度、填洼的结果会翻车。下面直接按我做项目的顺序写怎么验证数据、怎么用 shp 裁剪、怎么算坡度坡向和填洼以及哪些坑是必踩的。2. 先做数据体检30m DEM 的坐标系、边界和像元信息怎么查拿到压缩包以后不要急着往 GIS 里拖。先花三分钟做一次数据体检确认 DEM 的坐标系、空间范围、像元尺寸、无效值再确认 shp 的坐标系和范围。原因是一半以上的后续问题——裁剪位置偏移、坡度结果异常、填洼范围错误——都出在这一步而不是后面的处理环节。2.1 30m 分辨率背后的体量一个像元代表多大地面30m 分辨率意味着每个像元大致对应地面 30 米乘以 30 米也就是约 900 平方米。南平市是福建省面积较大的地级市面积差不多在 2.6 万平方千米量级。简单算一下2.6 万平方千米除以 900 平方米大约是 2900 万个像元。按 Float32 类型估算未压缩的 GeoTIFF 大约需要上百 MB如果带压缩且转成合适的整型通道实际落地通常在几十 MB 量级。这个体量对单机 QGIS 或 Python 都很友好不需要上分布式集群。30m 数据能做什么不能做什么要先有边界感。市域尺度的坡度分级、汇水范围勾画、灾害风险初步评估它够用道路选线和场地平整的精细断面它会显得粗糙。另一个容易忽略的点是30m DEM 对山谷、山脊的刻画已经能体现地形趋势但小尺度的微地形起伏会被平均掉。所以项目需求如果写着“精确到地块”就要另找更高精度数据不要再在这个包上较劲。2.2 用 rasterio 和 geopandas 读取元数据一个脚本做完体检数据体检我一般直接用 Python一次性把 DEM 和 shp 的元数据都打出来。下面这个脚本是固定开头import rasterio import geopandas as gpd dem_path 福建南平DEM30m.tif shp_path 南平市区域范围.shp with rasterio.open(dem_path) as ds: print(CRS:, ds.crs) print(像元尺寸:, ds.res) print(行列数:, ds.width, x, ds.height) print(范围:, ds.bounds) print(NoData:, ds.nodata) print(数据类型:, ds.dtypes[0]) gdf gpd.read_file(shp_path, encodingutf-8) print(shp CRS:, gdf.crs) print(shp 要素数量:, len(gdf)) print(shp 范围:, gdf.total_bounds)这个脚本解决四个问题。第一确认 DEM 的坐标系看它到底是投影坐标还是地理坐标分辨率输出是米还是度第二确认 NoData 值是什么后面裁剪和统计都要用它来区分无效像元第三看一共几个面要素区域 shp 可能是一个整体面也可能是多个区县面第四把 shp 外接矩形和 DEM 范围放一起对比初步判断两者是不是真的叠得上。特别注意ds.res这个输出。如果看到类似(0.00027, 0.00027)这种小数点后面一串的数值说明 DEM 还是地理坐标系横向和纵向单位是度而不是米。这种数据后续算坡度一定要先重投影否则水平距离和高程单位对不上坡度会算出一堆伪 89 度这个坑后面专门讲。2.3 当 shp 和 DEM 坐标系不一致用 to_crs 统一后再比较范围shp 文件和 DEM 坐标系不一致是常态。常见组合是 DEM 用 UTM 投影而 shp 用 CGCS2000 或 WGS84 地理坐标。直接肉眼对比两个外接矩形没有意义要先把 shp 转到 DEM 的坐标系里再判断重叠关系。from rasterio.warp import transform_bounds with rasterio.open(dem_path) as ds: dem_crs ds.crs dem_bounds ds.bounds if gdf.crs is None: gdf gdf.set_crs(EPSG:4326) shp_in_dem_crs gdf.to_crs(dem_crs) print(shp 转到 DEM 坐标系后的范围:, shp_in_dem_crs.total_bounds) print(DEM 范围:, dem_bounds)这里的逻辑是先把 shp 强制指定一个合理坐标系如果没有.prj文件gdf.crs会显示 None再通过to_crs(dem_crs)投影到 DEM 坐标系最后比较两个范围。如果 shp 在 DEM 坐标系下的范围远大于 DEM 范围说明前者覆盖的是更大区域不能直接拿来做细节统计如果两者南辕北辙多半是 shp 的原始坐标系被猜错了。提示如果gdf.crs是 None不要急着猜 EPSG先打开 shp 配套的.prj文件看内容没有.prj文件时尽量从数据来源处确认坐标系猜错了后续所有裁剪结果都会跟着错。3. 用区域边界 shp 裁剪 DEM从整幅影像到市域范围3.1 为什么剪不剪差别很大外接矩形、镶嵌边与 NoData有相当一部分人拿到 30m DEM 后直接加载到 QGIS 里就开始做坡度不裁剪。如果数据源本身只覆盖南平市范围问题不大但如果原始 tif 是更大范围的瓦片裁剪出来的问题就来了。不裁剪意味着整个外接矩形都会参与后续栅格计算包括南平市边界之外的邻区数据以及大量边缘 NoData。做过统计的人会看到海拔最小值、最大值、平均值全部被污染边界上的无效像元可能被当成 0 米或某个 NoData 值参与计算。包内单独放一个区域范围 shp本质就是让你拿着它做掩膜裁剪。裁剪不只是切掉边缘更关键的是把“不属于分析范围”的像元全部标记成 NoData让后续坡度、坡向、填洼只在有效范围内计算。否则边界处会出现一圈奇怪的陡坡因为紧挨着 NoData 的像元在计算领域窗口时缺邻居结果会被拉高。3.2 gdalwarp 快速裁剪cutline 参数是核心命令行裁剪我用得最多的是 gdalwarp它自带 cutline 支持一条命令同时完成裁剪和 NoData 设置。gdalwarp -cutline 南平市区域范围.shp -crop_to_cutline \ -dstnodata -9999 -overwrite \ -co COMPRESSDEFLATE -co TILEDYES \ 福建南平DEM30m.tif 南平市DEM_cut.tif这里几个参数都是必须的。-cutline指定边界 shp-crop_to_cutline让输出范围严格贴合 shp 边界而不是外接矩形-dstnodata -9999把所有无效像元统一写成 -9999避免原数据里 NoData 是 0 或者 65535 之类的值干扰后续判断。-co COMPRESSDEFLATE做无损压缩-co TILEDYES开启分块存储后续读像元时性能更好。如果不加-crop_to_cutline输出还是会按外接矩形给一整块区域边界外的像元全填 NoData。表面上也能用但后续坡度坡向计算会把边界外 NoData 当作领域参与运算导致边界一圈出现明显错误条带。所以这条建议加上。另一个注意点是 shp 自带.prj时 gdalwarp 会自动转换到 DEM 坐标系如果 shp 没有投影信息会默认用地理坐标系猜测结果大概率错位。3.3 用 Python 做掩膜裁剪顺便统一 NoData如果后续还要按多个区县分别出成果我建议直接用 rasterio 的 mask 接口在 Python 里把裁剪和 NoData 设置一起搞定方便批量循环。import rasterio from rasterio.mask import mask import geopandas as gpd gdf gpd.read_file(南平市区域范围.shp, encodingutf-8) if gdf.crs is None: gdf gdf.set_crs(EPSG:4326) with rasterio.open(福建南平DEM30m.tif) as src: gdf_proj gdf.to_crs(src.crs) geoms [geom for geom in gdf_proj.geometry if geom is not None] out_image, out_transform mask( src, geoms, cropTrue, nodata-9999, all_touchedTrue ) out_meta src.meta.copy() out_meta.update({ driver: GTiff, height: out_image.shape[1], width: out_image.shape[2], transform: out_transform, nodata: -9999, }) with rasterio.open(南平市DEM_cut_py.tif, w, **out_meta) as dst: dst.write(out_image)mask是 rasterio 的核心裁剪函数传进去的是面几何列表。cropTrue表示按几何裁剪而不是把几何外的像元全设 NoDatanodata-9999强制输出 NoDataall_touchedTrue会把与边界相交的所有像元都保留下来避免贴边像元被剔除后产生一圈空洞。out_meta是在原元数据基础上更新宽、高、变换信息写文件时保持坐标系不变。注意gdf.to_crs(src.crs)这步不能省shp 和 DEM 坐标系不一致时几何必须转成和栅格相同的坐标参考裁剪位置才准确。3.4 裁剪后立即检查NoData 占比和海拔统计裁剪完成不代表结果正确。我每次做完裁剪都会补一个检查脚本统计 NoData 占比和有效高程范围。import numpy as np import rasterio with rasterio.open(南平市DEM_cut_py.tif) as ds: arr ds.read(1) nodata ds.nodata if ds.nodata is not None else -9999 valid arr[arr ! nodata] nodata_ratio (arr nodata).sum() / arr.size print(有效像元:, valid.size) print(NoData 占比: {:.2%}.format(nodata_ratio)) print(海拔范围:, valid.min(), valid.max()) print(平均海拔:, valid.mean())NoData 占比超过 3% 就要警惕。要么是 shp 和 DEM 范围没有真正叠上要么是原始 DEM 在边界处存在空洞。海拔范围也要看合不合理南平这类闽北山区常见海拔从几十米到一千多米都有如果范围异常小多半是 NoData 设置把有效值也一起覆盖了。这一步发现问题后面还有后悔药吃直接拿去算坡度很难回头定位是哪一步出了问题。4. 把 DEM 变成地形因子坡度、坡向、填洼和出图流程4.1 先重投影再算坡度水平单位必须和高程单位对齐这是 DEM 处理里最容易踩的暗坑。原始 DEM 如果是经纬度坐标系横向单位是度纵向高程单位是米坡度计算公式里水平距离和垂直距离量纲不统一算出来的坡度几乎全部偏大丘陵都能算出 80 多度。所以算坡度之前先确认ds.res的输出。如果看到的是度而不是米必须先重投影。常见做法是转到一个本地投影坐标系。南平市大约处于东经 118 度附近UTM 50N 是一个可用的选择EPSG 编码是 32650。如果你使用的是 CGCS2000 框架的数据选对应的 CGCS2000 3 度带高斯投影也可以原理相同只是参考框架和带号不同。gdalwarp -t_srs EPSG:32650 -r bilinear -dstnodata -9999 \ -co COMPRESSDEFLATE -overwrite \ 南平市DEM_cut_py.tif 南平市DEM_prj.tif-t_srs指定目标坐标系-r bilinear采用双线性重采样。高程是连续表面双线性插值比最近邻更平滑有效避免地形台阶感。完成后可以再打印一次ds.res确认分辨率输出已经是米。如果显示 30 左右说明投影关系正确。4.2 gdaldem 计算坡度坡向两个命令拿到基础地形因子重投影之后直接用 gdaldem 算坡度坡向。gdaldem slope 南平市DEM_prj.tif 南平市坡度.tif gdaldem aspect 南平市DEM_prj.tif 南平市坡向.tif坡度输出的单位是度范围 0 到 900 表示平地90 表示垂直崖壁。坡向输出范围 0 到 3600 度表示正北90 度表示正东顺时针方向。拿到这两个文件后我一般会在 QGIS 里做一次可视化检查坡度图用白到红的单色渐变坡向图用分类色带。如果坡度图上大片区域显示 80 度以上回到 4.1 检查重投影这步如果坡向图出现大量 -9999 像素说明 NoData 没处理好。坡度坡向结果能不能直接用于工程计算取决于单位。项目需要百分比坡度时gdaldem 可以加-p参数输出百分比需要平均坡度时建议先把 NoData 掩膜统计出来再计算有效像元平均值别用原始栅格直接做全局统计。4.3 填洼汇水分析前必做的一步DEM 里存在大量伪洼地可能是插值误差、噪声或者真实地形中细小的凹陷。直接做流向分析水流会陷在洼地里出不来汇水路径断成一截一截。所以水文分析前先填洼是固定流程。常见做法是用 WhiteboxTools 的 FillDepressions命令比较直接whitebox_tools --runFillDepressions \ --input南平市DEM_prj.tif \ --output南平市DEM_filled.tif填洼把每个洼地填到它的最低溢出点让地表水流可以连续往外排。输出结果是浮点型 DEM只用于后续流向分析和汇水累积不要拿它再去算坡度和坡向。填洼后的 DEM 地形已经被人为修改过坡度值会明显偏缓用在需要真实地形表面的场景里会失真。4.4 山体阴影和配色让成果能放进报告原始 DEM 直接显示是灰蒙蒙一片看不出地形起伏。报告里通常叠加山体阴影提升立体感。gdaldem hillshade 南平市DEM_prj.tif 南平市山体阴影.tif -az 315 -alt 45 -z 1.5-az 315是光照方位角默认西北光适合表达山地-alt 45是太阳高度角高度角越低阴影越长-z 1.5是垂直放大系数地形起伏不够时可以适当调大。出图时在 QGIS 里把彩色 DEM 放底层山体阴影放上层设置混合模式为正片叠底再调整透明度就能得到一张既有高程分层又带立体感的底图。5. 避坑手册30m DEM 项目里最容易翻车的四个环节5.1 裁剪结果出现黑边边界像被狗啃现象裁剪出来的 tif 加载到 GIS 里四周有一圈黑色或者边界处锯齿状空缺统计面积时明显不对。原因一是 gdalwarp 时漏了-crop_to_cutline输出还是外接矩形边界外全是 NoData二是 Python 的mask没开all_touchedTrue贴边像元被当成越界剔除三是原始 DEM 的 NoData 值是 0 或者 255没有在裁剪时重设成负数显示时被当成黑色。解决命令行补上-crop_to_cutlinePython 端补上all_touchedTrue同时统一用-dstnodata -9999重设无效值。裁剪完成后跑一遍 3.4 的检查脚本NoData 占比和有效范围都在预期内再往下走。5.2 坡度结果大片接近 89 度和实际地形完全不符现象坡度图大面积高值丘陵地带也显示 80 度以上完全不是印象里的地形。原因DEM 还是地理坐标系横向单位是度而高程单位是米。坡度算法用邻域像元的水平距离做分母水平距离被当成了 1 度等于 1 米结果自然全部失真。这种情况在我接手过的数据里出现过不止一次几乎成了定式。解决算坡度之前先检查ds.res如果输出是度按 4.1 的命令重投影到 EPSG:32650 或当地投影坐标系。重投影后重新打印分辨率确认变成米级数值再做坡度。不要用gdaldem slope -s 111120这种方式硬拉比例尺那只是把错误结果做了一次单位换算不如直接换投影坐标系干净。5.3 shp 和 DEM 明明在一个区域裁剪结果却是全空或明显偏移现象两个文件在 GIS 里叠加看起来在同一区域但 gdalwarp 跑完输出要么全黑要么整体偏移一大段距离。原因shp 缺失.prj文件gdf.crs为 None程序按 EPSG:4326 猜但 shp 实际是 CGCS2000 或某个投影坐标或者 shp 有投影信息但 DEM 的坐标系定义不完整裁剪时几何转换发生了位置漂移。解决先用 2.2 的脚本把 shp 的 CRS 和范围打出来再和 DEM 范围做一次to_crs后的对比。给 shp 设置坐标系时不要靠猜优先找数据来源说明。裁剪命令里确保 shp 通过to_crs(src.crs)转到 DEM 坐标系再传给 mask。5.4 填洼把整个流域填成高原等高线全部失真现象填洼后的 DEM 在平缓地区也变成大片平地等高线扭曲后续坡度分析结果偏缓。原因FillDepressions 默认把所有洼地都填到最低溢出点不管它是伪洼地还是真实存在的天然洼地。真实地形里有湿地、湖泊、封闭盆地直接全填会把真实地形特征抹掉。解决填洼只用于水文分析的流向计算不要拿填洼结果去算坡度、坡向或做等高线。需要保留部分地形特征时可以改用带深度限制的填洼算法比如最大填充深度阈值把小于阈值的洼地填平大于阈值的天然洼地保留。处理前用 3.4 的统计脚本看一眼高程直方图确定一个合理的阈值范围而不是无脑默认参数。6. 让 30m DEM 进入生产管线虚拟镶嵌与批量裁剪6.1 先用 VRT 做虚拟镶嵌不要着急 merge当数据不止一个压缩包比如周边区域也拿到了同样规格的 30m DEM需要合成一块大范围地形底图时我一般不直接 merge而是先构建 VRT。VRT 是一个虚拟栅格目录不真正拷贝像元只是记录每个瓦片的范围和分辨率几十幅 tif 构建 VRT 几乎是瞬间完成。gdalbuildvrt 南平市_周边_虚拟.vrt 瓦片1.tif 瓦片2.tif 瓦片3.tif gdalwarp -cutline 南平市区域范围.shp -crop_to_cutline \ -t_srs EPSG:32650 -dstnodata -9999 \ -co COMPRESSDEFLATE \ 南平市_周边_虚拟.vrt 南平市DEM_final.tif最后一步才真正输出目标文件同时完成拼接、裁剪、重投影、压缩四个动作。这样做的好处是中间不产生冗余的大文件磁盘占用小如果某个瓦片有问题只需要重新 build 一次 VRT不用重跑整个 merge 流程。6.2 按 shp 属性字段批量输出各区县子集项目里经常需要按区县分别提交地形数据包。手动逐幅裁剪很费时间我一般在 shp 属性里找一个代表区县名称的字段用 groupby 循环批量生成。import geopandas as gpd import rasterio from rasterio.mask import mask gdf gpd.read_file(南平市区域范围.shp, encodingutf-8) key 区县名 # 按实际 shp 字段名调整 with rasterio.open(南平市DEM_prj.tif) as src: for name, sub in gdf.groupby(key): geoms list(sub.geometry) out_image, out_transform mask( src, geoms, cropTrue, nodata-9999, all_touchedTrue ) meta src.meta.copy() meta.update({ height: out_image.shape[1], width: out_image.shape[2], transform: out_transform, nodata: -9999, }) with rasterio.open(f南平市_{name}.tif, w, **meta) as dst: dst.write(out_image)这个脚本按区县名逐幅裁剪并输出独立文件文件名里带上区县名交付时不需要额外写说明。需要注意groupby前确认字段值没有重名或空值否则输出文件会互相覆盖。批量处理多幅 DEM 时建议在循环里也加上 3.4 的检查逻辑发现 NoData 异常立刻打印文件名避免最后交付时才发现某幅图是坏的。我自己的固定流程现在很简单解包、读元数据、统一 CRS、裁剪、检查 NoData、重投影、算坡度坡向。这个顺序看着繁琐但能挡住绝大多数返工。早年间某个模拟项目就因为少做了 NoData 统计拿着未裁剪的影像直接做汇水区结果边界外大量无效像元被算进低洼区域整个汇水路径全偏了。数据本身不复杂复杂的是每一步都默认它是对的。希望帮到你。本文还有配套的精品资源点击获取