黔东南30m DEM数据处理全流程：从裁剪到分区统计实践

简介贵州省黔东南苗族侗族自治州30米分辨率DEM数字高程数据包内含州级范围shp矢量边界面向GIS从业者、城乡规划师、地理研究人员等可直接用于地形分析、流域模拟、灾害评估与地图制图免去自行下载拼接与裁剪的环节。压缩包共12个文件约110.95MB以tif栅格高程数据为核心配以shp边界及dbf属性表、prj坐标系定义、tfw定位信息、xml元数据等配套文件结构完整可在ArcGIS或QGIS中直接加载使用。边界文件覆盖范围略超出市界能保证接边连续性适合完整区域场景的专题分析。已有342人学习该数据包对需要黔东南州高精度地形底图与行政边界数据的用户是一份即取即用的基础地理信息资料。1. 拿到黔东南 30m DEM 数据包先想清楚要拿它做什么解压这个带地名的 30m DEM 数据包真正花时间的往往不是下载本身而是解压后的第一小时坐标系对不对、市级边界能不能对齐、NoData 是多少、shp 字段够不够用。对于做山地地形分析的工程师来说这包数据基本决定了后续坡度、坡向、山体阴影和分区统计的结果质量。30m 分辨率的数字高程模型在黔东南这样地形起伏明显的地区既有足够的宏观轮廓又不会像 12.5m 或 5m 数据那样动辄几个 GB 难以处理是区域级项目的常见选型。这篇文章围绕“贵州省黔东南苗族侗族自治州 DEM 数字高程数据 30m含市级范围 shp 文件.zip”这一包数据把从解压检查到最终产出统计表的完整路径过一遍适合准备用 DEM 做选址评估、水文分析或可视化底图的数据工程师。2. 打开 zip 先别急着分析检查 DEM 与市级 shp 的文件组成和坐标系拿到 zip 后的第一个操作不是 gdalwarp也不是打开 GIS 软件而是先看清楚压缩包里到底有什么。名称里写着“含市级范围 shp 文件”但实际压缩包内可能是多个分幅 tif、一个县域合并结果也可能是带 .tfw 世界文件的 img。数据组织方式决定你后续是直接裁剪还是需要先拼接。2.1 30m DEM 与 DSM 的差别决定你是否用错了数据数字高程模型DEM描述的是裸地表高程剔除了植被和建筑。与之对应的是数字表面模型DSM记录地表物体顶面包括树冠和房顶。你在做淹没分析、通信基站覆盖、输电线路选线时需要的是 DEM做城市天际线或森林冠层研究时才会用到 DSM。很多数据源同时发布 DEM 和 DSM如果从标题里无法确认解压后要立刻看附带的 xml 或 txt 元数据文件确认“30m”指的是空间分辨率而不是高程精度。这包数据名为“30m”最常见的来源是 SRTM 或 ALOS 的衍生版本。SRTM 原始分辨率约 1 弧秒约 30m而 ALOS 有 12.5m 版本覆盖赤道附近区域时精度更细。如果你在黔东南这种山高谷深的区域做小流域分析30m 能识别主沟道但会平滑掉部分次级冲沟如果后续要提取精细河网可以考虑再补一份 12.5m 数据做对比而不是只用单一来源。2.2 用 gdalinfo 核对分辨率、范围和 NoData 的三个关键字段解压到本地目录后先用命令行把影像元数据打出来。相比直接拖进 QGISgdalinfo 的好处是输出稳定、便于写入脚本也能在第一时间发现文件损坏或坐标系缺失的问题。cd /data/qdn_dem ls -lh gdalinfo qdn_dem_30m.tif | head -60head -60 截断输出重点看以下字段Size 显示栅格宽高Pixel Size 里第一个数如果接近 0.00027777781 弧秒说明数据是地理坐标系下的 30m 格网如果接近 30 且单位是米则是投影坐标系。坐标系段Coordinate System会写明 EPSG 编号这是后续所有重投影操作的基准。字段含义踩坑点Size栅格行列数行列数异常小说明可能只是分幅的一部分Pixel Size像元尺寸数值带 0.0002x 秒量级时别忘了转投影NoData Value无效值标记常见为 -9999 或 -32768统计前必须排除Type数据类型Float32 常见少数为 Int16影响体积与精度2.3 用 ogrinfo 查看市级 shp 字段弄清“市级范围”到底有多全DEM 栅格本身没有属性表真正的地名信息全在 shp 里。黔东南州下辖县级行政区这个“市级范围”shp 的字段设计决定了你后面能按什么粒度统计。ogrinfo -so -al qdn_city.shp输出中 Feature Count 是面要素个数Geometry Column 通常是 Polygon 或 MultiPolygon属性字段列表里如果只有 NAME 和 CODE那就是标准行政区边界如果还有 AREA、PERIMETER 或者城市等级字段说明数据经过预处理。建议顺手看一眼每个要素的范围确认 shp 的包围盒比 DEM 小。若 shp 范围明显大于州界说明文件里包含了相邻区域的面统计时要用 SQL 过滤。2.4 中文文件名在 Linux 下解压乱码的根因与处理这个 zip 的文件名带中文和括号在 Windows 下用默认资源管理器解压通常没问题但到了 Linux 服务器上用 unzip 直接解压大概率会出现乱码。原因是 zip 规范里没有强制指定文件名编码Windows 下压缩时常用 GBK而 Linux 的 unzip 默认按 UTF-8 解码头部的 flag。两种编码不对齐文件名就变成一串无效字符。unzip -O gbk 贵州省黔东南苗族侗族自治州DEM数字高程数据30m含市级范围shp文件.zip如果 unzip 版本不支持 -O 参数用 7z 也可以7z x 文件名.zip配合-mcp936指定代码页。更稳妥的做法是写一小段 Python用 zipfile 读取后按目标文件名重建压缩包把文件名统一转成 UTF-8。import zipfile, shutil src 贵州省黔东南苗族侗族自治州DEM数字高程数据30m含市级范围shp文件.zip with zipfile.ZipFile(src) as zin: names zin.namelist() for name in names: # 在 Windows 上先解压出乱码名再用原始名替换 raw name.encode(cp437).decode(gbk) print(raw)这段代码的用途是探测真实文件名。zipfile 在读取本地文件头时如果遇到非 ASCII 名会用 cp437 解码再用 gbk 重新编码就能还原中文。确认了原始名之后可以直接把单个条目解压到指定文件名避免后续所有脚本都带着乱码路径。3. 用市级 shp 把 DEM 剪出来裁剪、重投影与 NoData 处理检查完元数据下一步就是把整幅地形数据裁剪到黔东南的行政区范围内。这里有两个常见误区一是忽略坐标系差异直接裁剪结果要么报错要么输出全黑二是裁剪后不检查 NoData导致后续坡度计算把无效值当作真实高程参与运算。3.1 坐标系不一致时裁剪结果会是空的或偏移shp 和 DEM 的坐标系如果不一致gdalwarp 不会报错而是把两个范围叠在一起找交集。地理坐标系和投影坐标系的单位不同范围边界对不上时可能只留下一条像素甚至整个输出都是 NoData。黔东南区域常见的三种坐标系如下表。EPSG坐标系名称使用场景EPSG:4326WGS84 经纬度原始 DEM 常见单位是度EPSG:4490CGCS2000 经纬度国内测绘成果常见EPSG:32648 / 32649WGS84 UTM 48N / 49N投影后做距离分析更合适黔东南经度跨度约 107°E 到 109°E横跨 UTM 48 和 49 两个分带直接选单一 UTM 带会导致东侧或西侧的形变偏大。如果只是做裁剪和面积统计用 EPSG:4490 或 EPSG:4326 足够如果要做坡度和坡向必须转投影坐标系但建议用 Albers 等积投影或兰伯特等角投影来覆盖整个州域而不是死盯 UTM。3.2 gdalwarp 与 rasterio 两种方式做 cutline 裁剪确定好目标坐标系后用 gdalwarp 完成裁剪是效率最高的方式。cutline 是裁剪边界crop_to_cutline 让输出范围贴合边界-dstnodata 把边界外的区域统一写为指定值。gdalwarp \ -t_srs EPSG:4490 \ -cutline qdn_city.shp \ -crop_to_cutline \ -dstnodata -9999 \ -co COMPRESSDEFLATE \ qdn_dem_30m.tif qdn_dem_clip.tif这里没有加 -tr 重采样参数是为了保持原始像元尺寸。如果 DEM 本身是 4326 下的 1 弧秒格网转成 4490 后像元仍是 1 弧秒左右不改变数据内容。加 -co COMPRESSDEFLATE 是为了减小输出体积山地 DEM 在裁剪后通常会有大量无效值压缩率很可观。在 Python 里做同样的事情可以拿到更多控制权比如在掩膜前对 shp 做 buffer或对裁剪结果做统计。import fiona import rasterio from rasterio.mask import mask with fiona.open(qdn_city.shp, r) as shp: geoms [feature[geometry] for feature in shp] with rasterio.open(qdn_dem_30m.tif) as src: out_image, out_transform mask( src, geoms, cropTrue, nodata-9999 ) 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(qdn_dem_clip_py.tif, w, **out_meta) as dst: dst.write(out_image)mask 函数接收的 geoms 是多边形列表cropTrue 表示只保留边界内的像元。这段逻辑适合嵌入自动化流程比如按县逐个输出裁剪结果时循环传入不同的要素。需要注意如果 shp 是多面要素栅格边缘会产生锯齿状像元这是栅格化边界时的正常现象不必追求视觉平滑。3.3 裁剪后检查 NoData 与文件大小别急着用裁剪完成后别急着生成坡度先做一次数值检查。用 gdalinfo 看 NoData 值是否与原始文件一致再统计有效像元数。这里有个容易忽略的细节原始 DEM 的 NoData 可能是 -32768Int16而口头上约定用 -9999。如果你的后续脚本写死 -9999实际读取的无效值还是 -32768统计结果会带上一批异常高程最大最小值全部失真。gdalinfo -stats qdn_dem_clip.tif输出中的 STATISTICS_MINIMUM 和 STATISTICS_MAXIMUM 能快速发现问题。如果最小值是 -9999 或 -32768 而不是接近 100 的海拔值说明 NoData 混入了有效值。此时不要改原始文件而是定义一个全局常量在坡度计算和分区统计时统一传给工具。文件大小也是检查依据一幅覆盖黔东南全州的 30m Float32 裁剪结果压缩后大约在 80 到 150MB 之间如果只有几 MB说明裁剪范围或像元尺寸有问题。4. 派生地形分析坡度、坡向与山体阴影的生成与参数调整裁剪出来的 DEM 是一张数字矩阵但大部分业务问题需要的是地形因子这里坡度多少、坡向朝南还是朝北、山体阴影是怎样的。生产环境里最常见的三个派生图层是坡度、坡向和山体阴影它们分别服务于建设选址、光伏朝向和可视化底图。4.1 从 DEM 能派生哪些图层生产中最常用哪三类坡度图反映地面倾斜程度单位是度或百分比直接影响工程土方量和径流速度。坡向图把坡面朝向归为八个方向加平地用于光照分析、植被分布判断。山体阴影模拟太阳光照射下的明暗关系是做地形渲染底图的标准素材。另外还有曲率和地形粗糙度前者用于识别山谷山脊后者用于地质灾害评估但这两个因子对 30m 分辨率相对敏感12.5m 或更细的数据源表现更好。用 30m 数据做坡度分析存在一个尺度效应真实坡面上细小的陡坎会被平均到整个像元内导致坡度峰值偏低。因此在看结果时要关注相对趋势而不是绝对数值。比如同一个 30m 数据集里北部区域坡度均值大于南部这个结论是可靠的但如果拿去和实测 1m 无人机数据的坡度对比偏差会很直观。4.2 gdaldem 三行命令产出坡度、坡向与山体阴影GDAL 自带的 gdaldem 工具是生成这三类图层的标准选择不需要装额外依赖。gdaldem slope qdn_dem_clip.tif qdn_slope_deg.tif gdaldem aspect qdn_dem_clip.tif qdn_aspect.tif -zero_for_flat gdaldem hillshade qdn_dem_clip.tif qdn_hillshade.tif -z 2.0 -az 315 -alt 45第一条命令生成以度为单位的坡度取值范围 0 到 90。第二条生成坡向0 表示北90 表示东-1 表示平地加上 -zero_for_flat 参数后平地统一输出 0。第三条生成山体阴影-z 是垂直方向放大系数-az 是太阳方位角315 即西北方向-alt 是太阳高度角。这三个参数直接影响立体视觉方位角决定阴影方向高度角控制阴影长度高度角越小阴影拉得越长。参数作用常见调整-p坡度单位改为百分比道路设计常用取代默认的度-s缩放系数经纬度数据必须设为 111120否则坡度偏小-z垂直夸张系数山体阴影常用 2 到 3-az太阳方位角315 度适合常规展示-alt太阳高度角45 度阴影适中默认也是 45这里最隐蔽的坑是 -s 参数。gdaldem 计算坡度时默认假设 DEM 的水平单位与垂直单位一致。如果输入数据是 EPSG:4326 经纬度X 和 Y 单位是度而高程单位是米两者数量级完全不同算出来的坡度会小到几乎无法区分。此时必须加-s 111120将度转换为米。4.3 批量生成时的参数细节与 Python 封装当你有多个县域的裁剪 DEM 时逐条运行 gdaldem 并不合理。可以用 Python 循环调用把命令参数集中管理。import pathlib import subprocess dem_files list(pathlib.Path(out).glob(*_clip.tif)) for dem in dem_files: slope dem.with_name(dem.stem _slope.tif) subprocess.run([ gdaldem, slope, str(dem), str(slope), -s, 111120, -co, COMPRESSDEFLATE ], checkTrue)这里的 -s 111120 只对经纬度坐标的 DEM 生效。如果输入已经是投影坐标再加这个参数会把坡度值放大 111120 倍结果完全失真。建议在循环里先读取 transform判断像元单位后再组装参数。另外裁剪后如果 NoData 值没有写进配置文件gdaldem 会出现边缘异常值生成的坡度图会在边界处产生一圈高值处理时用 3.3 节的方法先统一 NoData。5. 分区统计让市级 shp 把 30m DEM 变成一张可汇报的高程统计表DEM 是栅格数据没法直接在 Excel 里汇报。真正给决策者看的往往是“某某区平均海拔多少、最大高差多少、陡坡区占多大比例”。这一步把市级 shp 作为统计单元用分区统计把栅格值聚合到矢量面上输出一张 CSV 就能完成从数据到报告的转换。5.1 为什么不用像素读高程而要用 zonal stats 做分区聚合如果你遍历某个市范围内的所有像元手动算平均值很快会遇到两个问题一是碰上边缘的无效像元需要额外过滤二是多个县域叠加分析时代码会膨胀。zonal stats 的思路是把每个多边形当作一个分组栅格像元按分组做聚合输出 min、max、mean、std 等指标。它本质上是一条 SQL 里的 GROUP BY只是分组键来自矢量边界。分区统计的可靠性取决于两个前提shp 与 DEM 坐标系一致NoData 被正确标记。如果 shp 是 GCS 而 DEM 是投影坐标可以先用 5.3 的 ogr2ogr 把 shp 也转成同一坐标系再做统计而不是依赖工具内部重投影。5.2 Python 里用 rasterstats 按市级边界输出 CSVrasterstats 是最省力的方案pip 安装后即可用。它会自动处理栅格与矢量的对齐并跳过 nodata。from rasterstats import zonal_stats import pandas as pd stats zonal_stats( qdn_city_4490.shp, qdn_dem_clip.tif, stats[min, max, mean, median, std, percentile_20, percentile_80], nodata-9999, all_touchedFalse, ) df pd.DataFrame(stats) df[city] [f[properties][NAME] for f in __import__(fiona).open(qdn_city_4490.shp)] df.to_csv(qdn_dem_by_city.csv, indexFalse, encodingutf-8-sig)zonal_stats 的第一个参数是矢量文件第二个是栅格文件stats 列表控制输出哪些统计量。nodata-9999 与前面裁剪时写入的值保持一致。all_touched 是个容易混淆的参数默认 False 表示只有像元中心落入多边形才参与统计改成 True 表示只要像元与多边形相交就计入。对于边界精确的行政 shp默认值更准确。percentile_20 和 percentile_80 比单纯的 min/max 更有业务价值因为 min 和 max 容易被个别异常像元带偏。用这两个分位数配合 mean 和 std能看出该区域高程分布是否均匀。CSV 用 utf-8-sig 编码这样直接双击打开也是中文正常显示。5.3 顺手把 shp 批量转成 GeoJSON 与 KML便于交付如果同事不熟悉 GIS却需要把边界叠加到 web 地图或导入草图大师GeoJSON 和 KML 是很常见的交付格式。ogr2ogr 一次转换一个文件批量转则用 shell 循环。for f in *.shp; do ogr2ogr -f GeoJSON ${f%.shp}.geojson $f ogr2ogr -f KML ${f%.shp}.kml $f done转换后建议检查每个 geojson 的坐标系声明。有些场景要求 WGS84 经纬度格式如果你的 shp 是 CGCS2000需要在 ogr2ogr 命令里加-t_srs EPSG:4326。KML 的字段名对中文支持一般属性里有中文长字段时偶尔会出现编码问题正规的做法是先转 GeoJSON再由前端工具生成 KML而不是直接从 shapefile 硬转。6. 重新归档 zip 前做三件收尾完整性校验、编码处理与一键复跑处理完坡度、坡向和统计表之后工作并没有结束。原始 zip 里的数据是上游提供的你的裁剪结果和分析产物也需要归档。这时候重新打一个 zip把源数据、中间产物、shp 和统计表放在一起不是保存一份备份而是让项目可复现。以下是三个值得做的收尾处理。第一校验裁剪结果的完整性。不要只靠肉眼用真实面积和像元数量交叉验证。import rasterio import geopandas as gpd with rasterio.open(qdn_dem_clip.tif) as src: arr src.read(1) valid (arr ! -9999) (arr 0) area valid.sum() * src.transform.a * src.transform.e / 1e6 print(fvalid pixel: {valid.sum()}, area: {area:.1f} km2)对比 gpd.read_file 的 shp 面积两者偏差在 1% 以内说明裁剪正常。30m 数据一个像元面积约 900 平方米计算出的公里数如果与行政区面积差异超过 2%就要回去看裁剪边界和 NoData。第二重新打包 zip 时统一编码和中国特色的文件名。Linux 下用 zip 命令压缩带中文名的文件Windows 用户打开时可能乱码。用 Python 的 zipfile 写入文件名时会自动设置 UTF-8 标志位Windows 10 以上的资源管理器可以正确处理。import zipfile, pathlib out qdn_dem_30m_final.zip files [qdn_dem_clip.tif, qdn_city.shp, qdn_dem_by_city.csv] with zipfile.ZipFile(out, w, zipfile.ZIP_DEFLATED) as zf: for f in files: zf.write(f, arcnamef)如果对数据保密有要求给压缩包加密码也在这个阶段完成。zipfile 模块支持设置密码但只对 ZIP_DEFLATED 有效默认加密算法是 ZipCrypto破解难度较低。更稳妥的是用 7z 命令生成 AES-256 加密的压缩包配合长度足够的密码比标准 zip 加密可靠得多。7z a -tzip -p -mheon qdn_dem_30m_final.7z qdn_dem_clip.tif qdn_city.shp qdn_dem_by_city.csv第三把整个流程固化成一条可复跑的脚本。数据从业者都知道两周后再看当时的处理过程常常说不清参数为什么这么定。在项目目录里放一个 build.sh注释写清裁剪坐标系、NoData 值和 gdaldem 缩放系数配合一张 checksum 清单归档才有意义。校验时用sha256sum -c checksum.txt一次性核对所有文件确认压缩包传输过程没有损坏。这套收尾做完原始的“贵州省黔东南苗族侗族自治州DEM数字高程数据30m含市级范围shp文件.zip”才算真正进入可用状态。本文还有配套的精品资源点击获取