简介辽宁省30米分辨率数字高程模型DEM是一份基于ASTER GDEM V3数据拼接生成的栅格地形数据集覆盖全省范围以GeoTiff格式存储采用WGS84坐标系。适合GIS研究人员、城乡规划师、环境与防灾工程人员用于地形分析、坡度坡向计算、流域水文模拟及地质灾害评估等场景。压缩包共包含10个文件总大小约190.28MB。其中主数据文件为30米分辨率GeoTiff栅格辅以tfw地理配准文件、prj投影信息文件以及一组shp/shx/dbf/sbn/sbx格式的辽宁省边界矢量文件便于在ArcGIS、QGIS等软件中直接叠加与裁剪分析。目前已有341人学习下载。数据来自ASTER GDEM V3经过多源遥感数据融合处理具有较高的可靠性边界矢量与栅格配套提供用户解压后可直接用于区域地形提取、可视性分析、淹没模拟等多种定量研究是开展辽宁省地理空间分析的基础数据支撑。1. 这份辽宁 DEM 到底能做什么拿到“辽宁DEM.rar”这个包时别被里面一堆 .sbx、.tfw、.aux.xml 后缀吓住。压缩包本质是两个 GIS 数据一个覆盖辽宁省的 30 米分辨率高程栅格LiaoNing_DEM_30m_ASTGTMV003.tif一套辽宁省边界矢量.shp 系列。它解决的是“没现成地形数据做分析”的问题——水文模拟、坡度计算、视域分析、规划选址这些活儿都需要先有 DEM 打好底子。我会按拿到资源后的处理顺序讲清楚坐标检查、拼接、裁剪、坡度计算和验证。适合 GIS 数据处理工程师、规划口的技术人员5 年以上的人也能在有坑的地方省时间。接下来直接从文件结构讲起。2. ASTER GDEM V3 的数据组织与坐标检查2.1 30 米分辨率到底意味着什么数字高程模型DEM里每个像素的高程值并不是一个点而是一个 30 米×30 米网格内的平均海拔。ASTER GDEM V3 由 NASA 与日本 JAXA 基于 ASTER 卫星立体像对生成覆盖全球陆地区域2019 年 8 月发布的这个包对应 V3 版本。需要说明的是它本质是“数字表面模型”DSM记录的是地表物体顶部的高度在森林覆盖区会比真实地面高出一截。辽宁东部山区天然林较多做地形切割深度时要预留这一误差不要把它当厘米级控制网成果。2.2 拆开压缩包看文件GeoTIFF 与 Shapefile 各自的职责压缩包内“辽宁DEM”目录下除了主文件 LiaoNing_DEM_30m_ASTGTMV003.tif还有一组同名辅助文件和一套辽宁省边界矢量。一个常见错误是以为 .tfw、.aux.xml、.vat.dbf 都必须留在原目录才能用实际上 GDAL 项目真正会读的只有 .tif 内嵌的地理信息和边界文件的 .shp/.shx/.dbf/.prj。文件类型处理中的角色LiaoNing_DEM_30m_ASTGTMV003.tifGeoTIFF 栅格高程主数据30m 分辨率.tif.aux.xml栅格辅助元数据记录统计值、色彩解释缺失不影响读取.tfw世界文件不带内嵌地理信息时才用本场景可忽略.tif.vat.dbf栅格属性表对连续高程栅格没有分析意义可忽略辽宁省.shp/.shx/.dbf/.prj矢量边界及坐标系定义用于裁剪、掩膜、面积统计辽宁省.sbn/.sbx空间索引老式 ArcGIS 生成QGIS/GDAL 无需这里特别提一下 .sbn/.sbx。看到压缩包里只有 .sbn 没有 .sbx或者两个都有都不影响后续操作。它们只是 ArcGIS 的二进制空间索引GDAL 读取 Shapefile 时从不依赖这对文件。你甚至可以只拷 .shp、.shx、.dbf、.prj 四个文件去其他环境。2.3 用 GDAL 读取元数据确认坐标系与范围在命令行先跑一次 gdalinfo是拿到数据的第一步。Windows、Linux 和 macOS 下 GDAL 均可一键安装这里以系统命令为最小复现路径gdalinfo LiaoNing_DEM_30m_ASTGTMV003.tif输出里重点看三处Coordinate System 是否为 WGS 84Size 的像素行列数是否与 30 米规格匹配NoData 值和计算出的最小最大高程是否合理。如果文件自带地理信息gdalinfo 报告中的GeoTransform会包含左上角坐标、像元宽度和像元高度。对 WGS84 地理坐标来说像元宽高是度数数值量级在 0.00027777781/3600左右对应约 30 米。如果希望把检查集成到自动化批处理中用 Python 更合适from osgeo import gdal ds gdal.Open(LiaoNing_DEM_30m_ASTGTMV003.tif) print(ds.GetProjection()) print(ds.RasterXSize, ds.RasterYSize) print(ds.GetGeoTransform()) band ds.GetRasterBand(1) print(band.GetNoDataValue(), band.ComputeRasterMinMax())这段代码先打开 tif再打印投影、行列数、六参数仿射变换以及第一波段的 NoData 和最小/最大值。GetGeoTransform返回六元组分别是左上角 x、东西方向像元分辨率、旋转项、左上角 y、旋转项、南北方向像元分辨率最后一项通常为负表示行方向从北到南。若输出里最小最大值出现 -9999 或 -3.4e38说明原始制作者把填充值写在数据内后续裁剪和坡度分析前需要先把 NoData 值设掉否则这些无效像素会被当成真实海拔参与计算。3. 把分块瓦片合并成辽宁省 DEM3.1 为什么原始数据要分块ASTER GDEM V3 全球数据下载时按 1°×1° 经纬度分幅每个文件 3601×3601 像素。辽宁省位于东经约 118.8°~125.8°、北纬约 38.7°~43.5°之间用地理坐标跨度计算至少要覆盖 7 列×5 行共 35 个瓦片实际边界图幅还不止这些。只要有多个分幅就会产生两个问题相邻瓦片之间有 1 行 1 列重叠不同图幅因成像时间、立体匹配条件不同同一位置的像素值存在系统性偏差。分开处理时山脊线两侧会出现“对不齐”的错位所以要先合并再让边界裁剪按统一坐标网格走。3.2 用 gdal_merge 做第一次拼接GDAL 里最直接的拼装工具是 gdal_merge.py。假设你已经把需要的瓦片放在同一个目录最简单的命令是gdal_merge.py -o Liaoning_merged.tif -ot Float32 \ -co COMPRESSDEFLATE -co TILEDYES -co BIGTIFFYES \ input_N41E123.tif input_N41E124.tif input_N42E123.tif ...参数作用要分清-ot Float32强制输出浮点避免多个整数栅格叠加时出现截断误差COMPRESSDEFLATE压缩率比 LZW 更稳适合高程数据TILEDYES会让 tif 内部按块组织读写后续坡度计算比条带存储快很多BIGTIFFYES则是为输出超过 4GB 做准备不设的话 4GB 边缘可能报错。这个工具的问题是它只做“简单覆盖”重叠区默认取后读入的值不会自动羽化。如果瓦片间差异大产物会看到明显拼接线。作为工程师我不推荐用 PS 式羽化更好的做法是用 gdalbuildvrt 先生成虚拟目录再用 gdalwarp 做真正的重采样把接边权重交给重采样算法。3.3 推荐路线VRT 虚拟拼接 重置 NoDataVRT 不复制像素只记录每块数据的空间位置和覆盖关系因此生成几乎是瞬间的。然后一次性转换成独立 GeoTIFFgdalbuildvrt -resolution highest -r average Liaoning.vrt N4*.tif N3*.tif gdal_translate -ot Float32 -a_nodata -9999 \ -co COMPRESSDEFLATE -co BIGTIFFYES \ Liaoning.vrt Liaoning_final.tif这里-resolution highest会让 VRT 以所有输入文件中最高分辨率作为输出分辨率-r average指定重叠区取像素平均这是最简单也最有效的接边平滑方式。如果相邻瓦片的高程差超过 20 米平均后会留下低对比的平滑条带但至少不会出现台阶状断错。-a_nodata -9999把原来的 0 或空值统一成 -9999给后续 gdalwarp 提供准确的判定条件。合并完成后用 gdalinfo -stats 检查整体数值分布gdalinfo -stats -hist Liaoning_final.tif对比输出中的 Minimum、Maximum、StdDev可以快速识别异常。以一个正常辽宁区域 DEM 为例大致会出现下面的量级检查项合理范围异常时处理最小高程0 附近沿海若为 -9999确认 NoData 设置最大高程1300 米以下若超过 2000 米查是否为残余噪声平均值200~500 米可对照省平均海拔判断StdDev150~300 米过低说明数据被过度平滑我一般会再额外跑一次 Python 直方图统计 -9999 的像素占比占比超过 0.1% 时那片区域多半是水体或云遮挡需要后续填洼前先做插值。合并文件如果已经是单一 tif则从下一步开始做边界裁剪就够了。4. 用辽宁省边界裁剪并计算坡度与山体阴影4.1 用 gdalwarp 按省界精确裁剪拿到“辽宁省.shp”后不要直接拿矩形窗口裁。省界裁剪必须考虑非常不规则的海岸线和辽西、辽东山地边缘使用 cutline 把栅格限制到面内gdalwarp -cutline 辽宁省.shp -crop_to_cutline \ -dstnodata -9999 -co COMPRESSDEFLATE \ Liaoning_final.tif Liaoning_clip.tif如果 shp 的投影与 tif 不一致GDAL 会自动在读取时做矢量重投影所以很多时候不需要手动先转换投影。但要注意-cutline指定的 shp 必须有正确的 .prj否则默认 WGS84可能与实际相差很大。若你希望边界外不保留多余的 0 值-crop_to_cutline是必须的不加它的话只按 cutline 做 alpha 掩膜输出范围仍是原矩形。4.2 坡度、坡向、山体阴影参数选型坡度计算必须处理“水平单位”问题。LiaoNing_DEM_30m_ASTGTMV003.tif 是 WGS84 地理坐标经度方向一个像元约 99 公里/度纬度方向约 111 公里/度如果不缩放GDAL 会直接把度当米得到的坡度会整体偏大几十倍。做法有两种一是先用 gdalwarp 重投影到 EPSG:32651UTM 51N单位为米再计算坡度二是在 gdaldem 里直接指定比例尺度。对辽宁地形分析我建议重投影因为后面所有基于距离的水文、可视域分析都会受益gdalwarp -t_srs EPSG:32651 -r bilinear -tr 30 30 \ -overwrite Liaoning_clip.tif Liaoning_clip_utm.tif gdaldem slope Liaoning_clip_utm.tif liaoning_slope.tif -p -compute_edges-tr 30 30强制输出 30 米格网-r bilinear用双线性内插避免过高频噪声。gdaldem slope输出默认是度加-p则输出百分比坡度。水文上常用百分比土建上常用度看下游用途决定。-compute_edges让边缘行也能参与计算避免裁掉最外一圈。坡向和山体阴影同样在 UTM 投影下执行gdaldem aspect Liaoning_clip_utm.tif liaoning_aspect.tif gdaldem hillshade Liaoning_clip_utm.tif liaoning_hs.tif \ -z 2.0 -azimuth 315 -altitude 45hillshade 参数里-z是垂直夸大系数。平原地区用 2 或 3 能突出微地形山地用 1 即可太大会让阴影过爆。-azimuth 315是太阳方位角西北光对东北中国的丘陵沟壑辨认度最好-altitude 45是太阳高度角45 度比较均衡想要更强的立体感改用 30 度。把这组参数表保存到项目文档里便于复现输出产品使用工具关键参数坡度(度)gdaldem slope无 -pUTM 坐标坡度(百分比)gdaldem slope-p坡向gdaldem aspect0 为北-9999 为平地山体阴影gdaldem hillshade-z 2 -azimuth 315 -altitude 454.3 面向水系模拟的填洼与流向如果目的是水文模拟坡度和坡向只是中间产品接下来要先填洼再算流向和汇流累计。WhiteboxTools 免费且适合批处理三个命令串起来whitebox_tools -r FillDepressions -v \ --demliaoning_clip_utm.tif --outputliaoning_fill.tif whitebox_tools -r D8Pointer -v \ --demliaoning_fill.tif --outputliaoning_ptr.tif --flow_typeMFD whitebox_tools -r D8FlowAccumulation -v \ --d8_pntrliaoning_ptr.tif --out_accumliaoning_acc.tifFillDepressions 第一个跑它把地形上的局部洼地填平防止水流断在坑里。D8Pointer 里的--flow_typeMFD是多流向法适合辽宁的丘陵漫流地形如果做沟道提取用单流向 D8 更稳。D8FlowAccumulation 得到的累计栅格单位是“汇入的格网数×像元面积”要转成汇水面积就乘 900 平方米。参数表方便对照步骤输入输出需要注意FillDepressionsDEMfilled.tif无需调阈值默认平面填挖D8Pointerfilled.tifptr.tifMFD 比单流向更平滑D8FlowAccumulationptr.tifacc.tif输出值是格网计数如果你拿到的是作者拼好的单一 tif可以直接从这个阶段开始前面的 gdal_merge 可以跳过。最后再强调一句填洼需要 DEM 水平单位是米如果保持 WGS84 坐标系的度填洼结果完全不可用。5. 数据可用性验证与导出技巧5.1 用 QGIS 叠加山体阴影在 QGIS 中把liaoning_hs.tif拖入再把liaoning_slope.tif或原始 DEM 放在其上。渲染时把上层色彩带设为半透明混合模式选“叠加”透明度设为 40% 左右。这个技巧能同时看清高程变化和地形骨架比单独看一个栅格信息量大得多。山体阴影被拉伸到 0~255不需要再做色带。5.2 用 gdallocationinfo 快速比对高程拿到已知高程点时可以用 gdallocationinfo 从命令行取值不打开 GISgdallocationinfo -valonly -geoloc liaoning_clip_utm.tif 512000 4580000这里给的是 UTM 投影坐标如果手头只有一个经纬度点123.5, 41.5先得投到 EPSG:32651。ASTER GDEM V3 垂直中误差一般在 5~10 米所以差个 3 米以内都算正常。如果出现 -9999说明那个点正好落在海洋或 NoData 区不能参与精度统计。5.3 发布前导出为 COG如果要把最终 DEM 放到内网地图服务或对象存储上建议转成 Cloud Optimized GeoTIFF。COG 的瓦片组织方式让服务器只读取用户请求的那部分web 端显示会明显变快gdal_translate -of COG -co COMPRESSDEFLATE \ -co OVERVIEWSIGNORE_EXISTING -co RESAMPLINGNEAREST \ liaoning_clip_utm.tif liaoning_cog.tifRESAMPLINGNEAREST是为了保持高程原始值不用双线性去修改数据。发布前再执行rio cogeo validate liaoning_cog.tif确认没有内部偏移问题。如果你不想装 rio-cogeogdalinfo liaoning_cog.tif看到内部结构里含 overviews 也算基本合格。本文还有配套的精品资源点击获取