简介江苏省连云港市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 或浮点型 Float32GeoTIFF 存储NoData 常用-9999或0。shp 那边用ogrinfo快速查看概要不需要打开 QGISogrinfo -so LYG_boundary.shp LYG_boundary-so表示只输出图层概要信息不遍历每个要素。输出中重点看Geometry类型和Feature Count。连云港市级边界一般只有一个 Polygon 面要素或少数多面要素如果 Layer 数目和字段异常说明 shp 可能带多个图层或属性表不完整。2.3 30m DEM 常见元数据怎么读拿到 gdalinfo 输出后经验值大概能帮你判断这套数据是否正常。下表是市级 DEM 交付包最常见的几个元数据形态元数据项常见值判断方法像素类型Int16 / Float32Int16 可存负高程Nodata 阶段如果出现 0 要警惕像素尺寸0.00027778 度左右地理坐标系下约等于 30m直接用经纬度表示NoData 值-9999 或 0海岸带区域可能有无效值0 不能简单当成海平面输出范围东经 118.7~120.5 度与连云港实际纬度区间做大致对比投影WGS84 / CGCS2000判断能否与边界 shp 搭配裁剪如果 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 COMPRESSDEFLATE \ -co TILEDYES \ 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, shapesaoi.geometry, cropTrue, nodata-9999, ) meta src.meta.copy() meta.update( driverGTiff, heightout_image.shape[1], widthout_image.shape[2], transformout_transform, nodata-9999, ) with rasterio.open(OUT_PATH, w, **meta) as dst: dst.write(out_image)这段代码先把 shp 转换到栅格的坐标系然后执行带裁剪的掩膜操作最后按原影像元数据写出新 GeoTIFF。相比gdalwarprasterio 更适合把裁剪嵌在批量循环里比如把全市切成多个分幅结果依次处理。注意 rasterio 的mask带cropTrue时依然输出一个矩形范围只是范围紧贴几何体边界矩形内边界外的像素由 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 工具集里有FillQGIS 的 SAGA 工具箱里有Fill Sinks用 Python 的话可以用 whitebox 或 richdem 库。以 ArcGIS 为例参数只需要一个Z limit单位与 DEM 高程单位一致。连云港沿海区域地势平缓海拔很低填洼时 Z limit 建议先设 5 到 10 米观察填洼结果的体积变化不要一上来就用大阈值否则会把真实洼地也填平。工具常用参数适用场景gdaldem slope-p / -s坡度分级、水土保持gdaldem hillshade-z / -az / -alt地形可视化ArcGIS FillZ limit水文分析预处理SAGA Fill SinksThreshold低洼平原地区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 COMPRESSDEFLATE LYG_dem_clip.tif LYG_dem_cog.tifCOG 格式通过内置概览和按需读取机制让服务端只传输屏幕范围内的瓦片而不是整幅影像。对于全市范围的 30m DEMCOG 化之后体积虽然没变但瓦片请求速度提升非常明显。最后再用 color-relief 做一张快速目检图gdaldem color-relief LYG_dem_cog.tif color_ramp.txt LYG_dem_color.tifcolor_ramp.txt是两列文本低值和高值分别指定颜色区间例如0 34 139 34表示 0 米处为绿色200 139 69 19表示 200 米处为棕色。生成后在 QGIS 中把 color-relief 放到 hillshade 上层并设置混合模式为 Vivid Light连云港的沿海平原、中部丘陵和北部山地边界会一目了然。这张图既是验收材料也可以直接作为汇报底图的原始素材。本文还有配套的精品资源点击获取