SpringBoot2+Vue3全栈开发来访管理系统实战
2026/9/16 3:20:26
简介：辽宁省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.tif | GeoTIFF 栅格 | 高程主数据，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 84；Size 的像素行列数是否与 30 米规格匹配；NoData 值和计算出的最小最大高程是否合理。如果文件自带地理信息，gdalinfo 报告中的GeoTransform会包含左上角坐标、像元宽度和像元高度。对 WGS84 地理坐标来说，像元宽高是度数，数值量级在 0.0002777778（1/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.4e+38，说明原始制作者把填充值写在数据内，后续裁剪和坡度分析前需要先把 NoData 值设掉，否则这些无效像素会被当成真实海拔参与计算。
3. 把分块瓦片合并成辽宁省 DEM
3.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 COMPRESS=DEFLATE -co TILED=YES -co BIGTIFF=YES \ input_N41E123.tif input_N41E124.tif input_N42E123.tif ...参数作用要分清：-ot Float32强制输出浮点，避免多个整数栅格叠加时出现截断误差；COMPRESS=DEFLATE压缩率比 LZW 更稳，适合高程数据；TILED=YES会让 tif 内部按块组织读写，后续坡度计算比条带存储快很多；BIGTIFF=YES则是为输出超过 4GB 做准备，不设的话 4GB 边缘可能报错。
这个工具的问题是它只做“简单覆盖”，重叠区默认取后读入的值，不会自动羽化。如果瓦片间差异大，产物会看到明显拼接线。作为工程师我不推荐用 PS 式羽化，更好的做法是用 gdalbuildvrt 先生成虚拟目录，再用 gdalwarp 做真正的重采样，把接边权重交给重采样算法。
3.3 推荐路线：VRT 虚拟拼接 + 重置 NoData
VRT 不复制像素，只记录每块数据的空间位置和覆盖关系，因此生成几乎是瞬间的。然后一次性转换成独立 GeoTIFF：
gdalbuildvrt -resolution highest -r average Liaoning.vrt N4*.tif N3*.tif gdal_translate -ot Float32 -a_nodata -9999 \ -co COMPRESS=DEFLATE -co BIGTIFF=YES \ 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 米 | 可对照省平均海拔判断 |
| StdDev | 150~300 米 | 过低说明数据被过度平滑 |
我一般会再额外跑一次 Python 直方图，统计 -9999 的像素占比，占比超过 0.1% 时，那片区域多半是水体或云遮挡，需要后续填洼前先做插值。合并文件如果已经是单一 tif，则从下一步开始做边界裁剪就够了。
4. 用辽宁省边界裁剪，并计算坡度与山体阴影
4.1 用 gdalwarp 按省界精确裁剪
拿到“辽宁省.shp”后，不要直接拿矩形窗口裁。省界裁剪必须考虑非常不规则的海岸线和辽西、辽东山地边缘，使用 cutline 把栅格限制到面内：
gdalwarp -cutline 辽宁省.shp -crop_to_cutline \ -dstnodata -9999 -co COMPRESS=DEFLATE \ 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:32651（UTM 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 | 无 -p，UTM 坐标 |
| 坡度(百分比) | gdaldem slope | -p |
| 坡向 | gdaldem aspect | 0 为北，-9999 为平地 |
| 山体阴影 | gdaldem hillshade | -z 2 -azimuth 315 -altitude 45 |
4.3 面向水系模拟的填洼与流向
如果目的是水文模拟，坡度和坡向只是中间产品，接下来要先填洼，再算流向和汇流累计。WhiteboxTools 免费且适合批处理，三个命令串起来：
whitebox_tools -r FillDepressions -v \ --dem=liaoning_clip_utm.tif --output=liaoning_fill.tif whitebox_tools -r D8Pointer -v \ --dem=liaoning_fill.tif --output=liaoning_ptr.tif --flow_type=MFD whitebox_tools -r D8FlowAccumulation -v \ --d8_pntr=liaoning_ptr.tif --out_accum=liaoning_acc.tifFillDepressions 第一个跑，它把地形上的局部洼地填平，防止水流断在坑里。D8Pointer 里的--flow_type=MFD是多流向法，适合辽宁的丘陵漫流地形；如果做沟道提取，用单流向 D8 更稳。D8FlowAccumulation 得到的累计栅格单位是“汇入的格网数×像元面积”，要转成汇水面积就乘 900 平方米。参数表方便对照：
| 步骤 | 输入 | 输出 | 需要注意 |
|---|---|---|---|
| FillDepressions | DEM | filled.tif | 无需调阈值，默认平面填挖 |
| D8Pointer | filled.tif | ptr.tif | MFD 比单流向更平滑 |
| D8FlowAccumulation | ptr.tif | acc.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 从命令行取值，不打开 GIS：
gdallocationinfo -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 COMPRESS=DEFLATE \ -co OVERVIEWS=IGNORE_EXISTING -co RESAMPLING=NEAREST \ liaoning_clip_utm.tif liaoning_cog.tifRESAMPLING=NEAREST是为了保持高程原始值，不用双线性去修改数据。发布前再执行rio cogeo validate liaoning_cog.tif，确认没有内部偏移问题。如果你不想装 rio-cogeo，gdalinfo liaoning_cog.tif看到内部结构里含 overviews 也算基本合格。
本文还有配套的精品资源，点击获取