30m DEM数据与市级shp处理：从解压到坡度生成全流程

简介：江苏省连云港市30米分辨率数字高程数据包，是面向地理信息系统应用的基础地形栅格数据，内含全市高程模型与配套行政边界矢量文件，可服务于测绘、城乡规划、水文分析及地理教学等场景，有效解决区域高精度地形数据获取困难的问题。压缩包共12个文件，大小仅14.54MB，核心内容为高程影像与市域范围矢量边界，同时附带坐标配准、投影定义、属性表、空间索引、金字塔与元数据等辅助文件，结构完整，可直接在ArcGIS、QGIS等主流GIS软件中加载使用。目前已有591人浏览学习，数据量适中，下载后即可用于实际项目或课程实验。借助这些数据，使用者无需另行采集，即可获得30米间隔的高程矩阵和精确的市域边界，方便进行坡度坡向分析、视域分析、流域提取等操作；配套的矢量范围也能辅助快速裁剪、拼接与出图，适合作为GIS入门练习或专业研究的可靠底图。

1. 为什么 30m DEM 总要配一个市级 shp 再交付

区域规划、水利普查或通信站点选址时，“连云港市 DEM 数字高程数据 30m + 本市级范围 shp”这个组合，本质上把两件事一次解决了：高程底图不用再到处找，行政边界也省去了自己描轮廓的环节。很多工程师拿到 zip 后第一反应是右键解压、拖进 ArcMap 就出图，但真正稳妥的流程应该是先解压确认文件结构，再看投影和无效值，最后才做裁剪与衍生分析。30m 分辨率在视觉上不如 12.5m 或 5m 细腻，但它的数据量小、覆盖稳定，在坡度分级、淹没分析、片区填挖方量估算这类中尺度任务里反而是性价比最高的选择。这篇内容就按我平时处理这种交付包的习惯，把解压、坐标对齐、裁剪、坡度山体阴影生成和验收验证一条线讲完。

2. 解压先摸底：zip 里的栅格与 shp 到底谁是谁

市级 DEM 交付包通常由若干 GeoTIFF 影像和一个完整的 shp 要素类组成。shp 不是单文件，而是一组同名文件的集合，最少包含 .shp、.shx、.dbf 三个文件，正规交付还会带上 .prj、.cpg 等附属文件。先弄清楚 zip 里装了什么，比直接双击打开更重要。

2.1 先解决 zip 解压的中文文件名乱码

国内数据商打包时常用中文或拼音命名文件，文件名一旦涉及中文，跨平台解压经常出乱码。Windows 下用系统自带解压或 7-Zip 时问题不大，但在 Linux 服务器或 Python 脚本里处理时，zipfile模块默认按 cp437 解码文件名，中文名会变成锟斤拷一类的乱码。我一般先在 Python 里列出压缩包内文件，确认实际文件名：

import zipfile with zipfile.ZipFile("lianyungang_dem_30m.zip") as zf: for info in zf.infolist(): print(info.filename, info.file_size)

如果输出乱码，就手动修正编码再解压。常见做法是用cp437重新编码，再按gbk解码：

import zipfile SRC_ZIP = "lianyungang_dem_30m.zip" with zipfile.ZipFile(SRC_ZIP) as zf: for info in zf.infolist(): raw_name = info.filename try: fixed_name = raw_name.encode("cp437").decode("gbk") except UnicodeDecodeError: fixed_name = raw_name print(f"解压: {fixed_name}") with zf.open(info) as src, open(fixed_name, "wb") as dst: dst.write(src.read())

这段代码先把zipfile按 cp437 读出的文件名还原成字节，再用 GBK 解码，适用于国内 Windows 压缩软件生成的包。如果压缩包本来就用 UTF-8 文件名，UnicodeDecodeError分支会跳过，不会影响解压。注意不要用zf.extract()直接解压，因为内部仍用错误编码写文件名，落盘后依旧是乱码。

2.2 用 gdalinfo 与 ogrinfo 确认 DEM 与边界的底细

解压完成后，先对 DEM 做一次摸底。GDAL 自带的gdalinfo能把影像的尺寸、波段、坐标系、像素类型、NoData 值一次性列出来：

gdalinfo LYG_dem_30m.tif

重点看四类信息：Size决定行列数和后续算力开销；Pixel Size决定地面分辨率；Coordinate System决定能否和 shp 直接叠加；NoData Value决定后续坡度计算时无效值会不会污染统计。常见的 30m DEM 交付格式是整型 Int16 或浮点型 Float32，GeoTIFF 存储，NoData 常用-9999或0。

shp 那边用ogrinfo快速查看概要，不需要打开 QGIS：

ogrinfo -so LYG_boundary.shp LYG_boundary

-so表示只输出图层概要信息，不遍历每个要素。输出中重点看Geometry类型和Feature Count。连云港市级边界一般只有一个 Polygon 面要素或少数多面要素，如果 Layer 数目和字段异常，说明 shp 可能带多个图层或属性表不完整。

2.3 30m DEM 常见元数据怎么读

拿到 gdalinfo 输出后，经验值大概能帮你判断这套数据是否正常。下表是市级 DEM 交付包最常见的几个元数据形态：

如果 DEM 是经纬度坐标的 GeoTIFF，像素尺寸会显示为度，比如(0.0002777778, -0.0002777778)，这代表每个像素约等于地表 30m。某些交付包将 DEM 直接做成 CGCS2000 3 度分带投影，像素尺寸会直接显示(30, -30)，这种数据与多数测绘成果叠加更方便。

3. 把边界 shp 转成 DEM 的裁剪条件：坐标系对齐是关键

很多人在裁剪时报错“两个图层不重叠”或裁出来是黑块，八成不是范围和边界的问题，而是坐标系没对齐。

3.1 投影坐标系与地理坐标系先统一

DEM 可能是经纬度坐标，而本市级范围 shp 来自国土或规划部门，常使用 CGCS2000 或西安 80 的高斯投影平面坐标。两者单位不同，经纬度单位是度，投影坐标单位是米，直接叠加时偏差可能达到数百米甚至上千公里。

常见的处理原则是：以 shp 的坐标系为准，将 DEM 重投影后再裁剪。先用 gdalinfo 查看 DEM 的坐标系，再用gdalsrsinfo读取 shp 的投影文件得到 EPSG 编号：

gdalsrsinfo LYG_boundary.prj -o epsg

得到 EPSG 编号后，在后续gdalwarp中用-t_srs指定。连云港所在的经度范围主要落在 3 度分带的 40 带附近，具体带号以 shp 的.prj内容为准，不要拍脑袋猜。

3.2 用 gdalwarp 一步完成裁剪与重投影

GDAL 处理 DEM 裁剪最常用的命令是gdalwarp，它可以把重投影、重采样、裁剪、压缩一次性完成：

gdalwarp -t_srs EPSG:4490 \ -cutline LYG_boundary.shp \ -crop_to_cutline \ -dstnodata -9999 \ -co COMPRESS=DEFLATE \ -co TILED=YES \ LYG_dem_30m.tif LYG_dem_clip.tif

命令中-t_srs指定目标坐标系，示例用 EPSG:4490 表示 CGCS2000 地理坐标系，实际使用时以 shp 的投影为准；-cutline传入边界 shp 路径；-crop_to_cutline让输出栅格的范围与矢量多边形完全一致，而不是只留外接矩形；-dstnodata把裁剪边缘外和原无效值统一写成 -9999；-co参数控制 GeoTIFF 的内部结构，DEFLATE 压缩可减少约一半文件体积。

这里有个容易混淆的点：-cutline只做外接矩形裁剪，多边形内部的保留逻辑仍依赖-crop_to_cutline。如果漏掉-crop_to_cutline，输出会是包含整个边界的矩形范围，边界外的像素保留为无效值，后续做坡度统计时会把这些空值一起算进去。

3.3 用 rasterio 在 Python 管线里完成相同操作

批量处理多个区块或要把裁剪流程集成到定时任务里时，我习惯用 Python 的 rasterio 和 geopandas。

import geopandas as gpd import rasterio from rasterio.mask import mask DEM_PATH = "LYG_dem_30m.tif" SHP_PATH = "LYG_boundary.shp" OUT_PATH = "LYG_dem_clip.tif" aoi = gpd.read_file(SHP_PATH) with rasterio.open(DEM_PATH) as src: # 关键：把边界数据转换到 DEM 的坐标系 aoi = aoi.to_crs(src.crs) out_image, out_transform = mask( src, shapes=aoi.geometry, crop=True, nodata=-9999, ) meta = src.meta.copy() meta.update( driver="GTiff", height=out_image.shape[1], width=out_image.shape[2], transform=out_transform, nodata=-9999, ) with rasterio.open(OUT_PATH, "w", **meta) as dst: dst.write(out_image)

这段代码先把 shp 转换到栅格的坐标系，然后执行带裁剪的掩膜操作，最后按原影像元数据写出新 GeoTIFF。相比gdalwarp，rasterio 更适合把裁剪嵌在批量循环里，比如把全市切成多个分幅结果依次处理。

注意 rasterio 的mask带crop=True时，依然输出一个矩形范围，只是范围紧贴几何体边界，矩形内边界外的像素由 nodata 填充。如果后续步骤只想要“多边形内部严格有效”，需要在坡度或坡度分析阶段再过滤一次无效值，或者改用gdalwarp的-crop_to_cutline输出带内 netCDF 外壳的掩膜版本。

4. 基于裁剪后 DEM 生成坡度、坡向与山体阴影

拿到裁剪后的 DEM，下一步通常是生成分析底图。gdaldem 是 GDAL 里处理 DEM 派生产品的工具，包含slope、aspect、hillshade、color-relief等模块，命令行一条就能出结果。

4.1 gdaldem slope 与角度/百分比输出的区别

gdaldem slope LYG_dem_clip.tif LYG_slope_degree.tif -p -s 111120

参数-p表示输出坡度百分比，若不加则默认输出角度制；-s是垂直比例因子，表示水平单位与垂直单位的比值。当 DEM 是经纬度坐标时，水平单位是度，垂直单位是米，直接算坡度会导致结果严重偏小，必须设置-s 111120，即 1 度约等于 111120 米。如果 DEM 本身就是投影坐标且单位为米，通常不加-s。

角度和百分比在实际使用中差别很大。角度制适合人直接读，比如 30 度坡；百分比则适合做栅格计算器里的分级公式，常用于水土保持和道路选线。若两个文件都要用，建议各生成一次，不要试图从角度结果手动换算。

4.2 hillshade 参数与渲染叠加

山体阴影不是分析数据，而是可视化底图，它通过模拟光照让地形起伏肉眼可见。标准命令：

gdaldem hillshade LYG_dem_clip.tif LYG_hillshade.tif -z 2.5 -az 315 -alt 45

-z是垂直拉伸系数，平原地区 DEM 高程差小，默认值 1 会让阴影平淡，2.5 到 3 之间通常能压出地形纹理；-az是光源方位角，315 度即西北方向来光，这是地图惯例，因为西北光下地形判读符合多数人视觉习惯；-alt是太阳高度角，45 度属于中等光照，能兼顾坡面细节和整体明暗。

在 QGIS 中把 hillshade 放到底层、DEM 配色影像放上层并设置透明度约 60%，是市政规划汇报里最常见的组合呈现方式。此时可以直接看出断层、河谷和山体走向，对数据质量做初步定性判断。

4.3 水文分析前的填洼处理

如果后续要做流向、汇水面积或淹没模拟，不能直接用原始 DEM。真实地貌里存在大量由数据噪声和拼接误差造成的伪洼地，填洼是水文分析的标准前置步骤。

ArcGIS 的 Spatial Analyst 工具集里有Fill，QGIS 的 SAGA 工具箱里有Fill Sinks，用 Python 的话可以用 whitebox 或 richdem 库。以 ArcGIS 为例，参数只需要一个Z limit，单位与 DEM 高程单位一致。连云港沿海区域地势平缓，海拔很低，填洼时 Z limit 建议先设 5 到 10 米，观察填洼结果的体积变化，不要一上来就用大阈值，否则会把真实洼地也填平。

5. 验收数据与嵌入工程：几个可复现的最小验证

正式使用这套数据之前，用一组小命令做验收，能避免在后续分析到一半时才发现数据有问题。

5.1 用 gdalinfo -stats 抓无效值和取值区间

gdalinfo -stats LYG_dem_clip.tif

输出的STATISTICS_MINIMUM、STATISTICS_MAXIMUM和STATISTICS_VALID_PERCENT能直接反映数据的质量。连云港最高点大约在云台山玉女峰附近，海拔约 625 米，沿海区域大量像素应接近 0 米。如果统计结果显示最低值低到 -3000 或最高值超过 9000，基本可以判断无效值没有被正确设置，或者原始数据混入了条带噪声。

若VALID_PERCENT明显偏低，说明边界内有大片 NoData，需要回退到裁剪步骤检查-crop_to_cutline和坐标系对齐过程。

5.2 用已知地标验证高程合理性

统计值只能说明数值范围，不能证明空间位置对。我一般会挑一个市内的山体或制高点，用gdallocationinfo直接读某个经纬度对应的高程值，与公开地理数据交叉验证：

gdallocationinfo -valonly -wgs84 LYG_dem_clip.tif 119.16 34.64

-wgs84表示输入坐标是经纬度，输出结果是该点的高程值。把一个已知地点的高程与 DEM 读取值做差，差值在几米到十几米范围属于正常，如果差出数百米就要重新检查投影带号或数据源本身。

5.3 转成 COG 并生成快速预览图

验收通过后，把成果转成 Cloud Optimized GeoTIFF，后续在 QGIS、GeoServer 或 Web 地图服务里加载都会更快：

gdal_translate -of COG -co COMPRESS=DEFLATE LYG_dem_clip.tif LYG_dem_cog.tif

COG 格式通过内置概览和按需读取机制，让服务端只传输屏幕范围内的瓦片，而不是整幅影像。对于全市范围的 30m DEM，COG 化之后体积虽然没变，但瓦片请求速度提升非常明显。

最后再用 color-relief 做一张快速目检图：

gdaldem color-relief LYG_dem_cog.tif color_ramp.txt LYG_dem_color.tif

color_ramp.txt是两列文本，低值和高值分别指定颜色区间，例如0 34 139 34表示 0 米处为绿色，200 139 69 19表示 200 米处为棕色。生成后在 QGIS 中把 color-relief 放到 hillshade 上层并设置混合模式为 Vivid Light，连云港的沿海平原、中部丘陵和北部山地边界会一目了然。这张图既是验收材料，也可以直接作为汇报底图的原始素材。

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