ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

福州市30米DEM数据处理:从GDAL解压、坡度分析到边界裁剪全流程

福州市30米DEM数据处理:从GDAL解压、坡度分析到边界裁剪全流程 简介福州市30米分辨率DEM数字高程数据是一套面向GIS分析、城市规划、环境评估、灾害模拟等场景的基础地形数据集。压缩包共12个文件约35.83MB核心为TIFF格式的高程栅格同时包含福州市行政范围Shapefile.shp及配套的.dbf属性表、.prj投影文件、.sbn/.sbx空间索引、.tfw坐标配准信息与多份xml元数据可在ArcGIS、QGIS等平台中直接读取使用。该数据以30米间距的网格记录每个像元的海拔值覆盖福州市域并略向外扩展为研究区域地形起伏、径流走向、坡度坡向及通视条件提供了量化底图。目前已有1300人学习下载。借助附带的范围文件用户可快速裁剪、统计或叠加其他图层适用于区域规划、地质灾害风险评估、三维地形展示及水利模拟等专业工作是福建省内地理信息项目落地时可靠的基础数据支撑。1. 福州市30m DEM数据包能做什么为什么优先选它30米分辨率的DEM听起来没有激光雷达“高级”但做区域尺度的地形分析它反而是最常用的底图。福州市这份压缩包把DEM栅格和市范围shp打包成套下载解压后直接拖进ArcGIS、QGIS不用再去遥感平台花几小时拼接和投影。它可以用于福州地区的洪水淹没模拟、道路选线、景观视廊分析也可以作为地形分类的基础输入。和ALOS 12.5米那种精细数据相比30m的数据量小跑水文模型和可视域分析更快对行政边界是否贴合影响不大。对沉浸多年的GIS开发来说它的价值在于边界shp和tif已经对齐省去了自己勾范围坐标的步骤对教学场景则是一套完整的文件格式教材。30米一个像元沟谷形态和城市街廓能大致辨认但细节不如5米数据如果你只做宏观研判这份数据足够。接下来按“文件结构、数据体检、地形分析、边界裁剪”四个步骤往下拆。2. 拆开zip包从tif、tfw到shp每个文件都别丢很多人解压后只拖走福州市DEM.tif剩下的shp和tfw当垃圾等下次打开发现DEM“漂移”了才回来找。zip里其实分两组一组是tif栅格及辅助文件一组是边界矢量。下面逐个说明为什么这些文件都值得保留。2.1 主数据文件福州市DEM.tif和它的辅助文件福州市DEM.tif是GeoTIFF高程值以像元为单位存储每个像元的数值代表海拔高度单位通常是米。30m分辨率意味着地面一个像元对应30米×30米的范围。要确认tif底层分辨率和坐标系最直接的是用GDAL自带的gdalinfogdalinfo 福州市DEM.tif | head -n 40 # head 用于截断输出避免元数据太长刷屏输出中Size is后面是列数和行数Pixel Size 括号里是x方向、y方向的分辨率Coordinate System is:后面是坐标系描述。如果看到EPSG:4326说明是经纬度像元大小约0.00027778度换算成米大约30.8米如果看到EPSG:4547这类投影坐标像元单位直接就是米。在tif同目录下还有福州市DEM.tfw这是世界描述文件。GeoTIFF本身通常已内嵌地理参考tfw可以看作备份当tif头部信息损坏或者要把tif导入CAD、Photoshop时tfw就派上用场。它共有六个数字依次是x方向像元大小、旋转项、y方向像元大小、左上角x坐标、左上角y坐标第2、3个旋转项一般都为0。福州市DEM.tif.ovr是金字塔文件用于缩放视图时读取低分辨率概览层加快显示福州市DEM.tif.aux.xml保存栅格统计信息和色彩映射。这两个文件不是必须的但建议保留频繁删除ovr会让大数据量tif在每次缩放时重新计算统计值明显卡顿。2.2 边界文件组福州市范围.shp和它的三件套Shapefile是矢量格式但它依赖多个文件同时存在。只有.shp是无法完整打开的至少需要shp、shx、dbf三件套。压缩包里的文件对应关系如下文件类型作用缺失时的影响福州市范围.shp核心矢量存储几何坐标面/多边形无法读取矢量福州市范围.shx索引记录几何位置偏移量加速读取部分软件拒绝加载福州市范围.dbf属性表保存地名、行政区代码等字段属性表为空或无属性福州市范围.prj投影描述保存坐标系WKT文本叠加时容易出现偏移福州市范围.sbn/.sbx空间索引在ArcGIS中加速空间查询可重建通常不影响打开福州市范围.shp.xml元数据记录数据来源和更新时间一般不影响使用福州市范围.sbn和.sbx是ArcGIS生成的二进制空间索引。QGIS即使没有它们也能打开ArcGIS里做属性查询时略慢但不会被判定为缺失文件。福州市DEM.tif.xml是tif的ISO元数据比如采集日期和坐标系描述对普通用户价值不大但归档时仍建议保留。2.3 为什么坐标系统一决定数据能不能对齐DEM是栅格shp是矢量两个图层能不能严丝合缝取决于坐标系是否一致。tif内部和shp的.prj必须指向同一个地理坐标系或投影坐标系否则在GIS里会看到DEM跑偏到海里。拿到数据后先做一次体检是解开一切问题的前提。常见做法是先把shp的.prj用文本编辑器打开看里面是GEOGCS[GCS_WGS_1984还是China Geodetic Coordinate System 2000再对照tif的投影信息决定是否需要重投影。2.4 截取边界坐标到txt方便脚本传参很多命令行工具不直接读shp只认矩形范围。这时需要从福州市范围.shp里取外接矩形边界。最快的方式是用ogrinfoogrinfo -so -al 福州市范围.shp | grep Extent # -so 表示只读概要-al 表示读取全部图层grep 只留下范围行输出是四个数字分别对应西、南、东、北。把数值按minX minY maxX maxY的顺序填到gdalwarp -te后面即可。这个方法在批量处理多个区县时非常省事不用每次都去QGIS属性里查看“图层范围”。如果只想把shp的属性表转成带经纬度的txt也可以用ogr2ogr -f CSV导出dbf内容但取范围用ogrinfo更快。3. 解压与体检用GDAL确认坐标系、范围和分辨率拿到zip包的第一反应不是双击解压到桌面而是先建一个英文目录把压缩包放进去用命令行释放尽可能减少中文路径带来的“找不到文件”问题。尤其是Windows下使用Python的rasterio或gdal库时中文路径偶尔会触发UnicodeDecodeError提前规避能省很多事。3.1 在Linux/macOS下解压不乱码mkdir -p ~/gis/fuzhou_dem cd ~/gis/fuzhou_dem unzip -O UTF-8 福建省福州市DEM数字高程数据30m.zip # -O 参数用于把压缩包内文件名按UTF-8重新编码如果解压后文件名出现乱码试试unzip -O GBKWindows下用系统右键“全部解压缩”即可。解压完成后建议把整个目录移动到类似C:\gis\fuzhou_dem这种不带中文的路径后续调用Python脚本会少很多问题。3.2 打开QGIS看一眼范围启动QGIS后直接把福州市DEM.tif拖进图层区再把福州市范围.shp拖入。如果DEM显示黑色或白色在图层样式的“单波段灰度”里把拉伸方法改成“最小/最大”渲染就会恢复正常。若shp和DEM完全重合说明坐标基准一致可以进入下一步地形分析。如果只是边界shp压住DEM而你想看单个区县的边界在属性表里找到行政区字段用表达式选中目标feature再右键“导出—保存所选要素”生成一个新shp。这和“只保留外边界线”是同一个逻辑用矢量—地理处理—消除多边形边界把内部边界消除而不是直接删字段。3.3 用GDAL做一次量化的数据体检第2章的gdalinfo已经够用再加两个参数做更完整的检查gdalinfo -stats 福州市DEM.tif | grep -E Upper Left|Lower Right|Pixel Size|EPSG|Statistics # -stats 会计算最小值、最大值和标准差grep 过滤出关键行需要重点看四组信息Upper Left和Lower Right决定范围确认是否完整覆盖福州且比行政边界略大因为原始下载是按矩形范围裁的。Pixel Size确认是否为30米级别如果出现0.00027或30.8说明坐标单位是度后续算坡度要加比例因子。EPSG代码决定后续重采样参数。Statistics里Min和Max如果出现类似-3.4e38的异常值说明tif内部有无效像素建议先执行gdal_translate -a_nodata -9999统一无效值避免坡度、等高线计算时把黑边一起算进去。3.4 坐标对不上时用gdalwarp重新投影如果DEM和shp之间存在明显偏移大概率是一个WGS84经纬度一个是CGCS2000高斯投影。这时先读取shp的.prj根据文本里的坐标系名称确认EPSG代码然后重投影gdalwarp -t_srs EPSG:4547 -r cubic -overwrite 福州市DEM.tif 福州市DEM_4547.tif # -r cubic 是三次卷积内插适合连续表面若只是显示用 -r bilinear 更快选择EPSG时福州地区常用CGCS2000 / 3-degree Gauss-Kruger zone 38对应的EPSG是4547。如果只是在网页端做了展示直接输出成EPSG:4326也能与WGS84坐标的shp对齐。注意重投影会重采样会轻微改变高程统计值尽量只做一次不要在多个坐标系之间反复切换。4. 从DEM到实用产品坡度、山体阴影和等高线输出DEM本身不是成品坡度图、山体阴影和等高线才是能交付的地形产品。GDAL自带一个gdaldem命令行工具一行命令就能生成多种地形分析图层比ArcGIS里逐一步骤操作更快也更适合写进自动化脚本。4.1 为什么要处理比例因子和边缘像元坡度算法的原理是围绕中心像元做一个3×3窗口用水平距离和垂直高差计算最大变化率。DEM如果是经纬度坐标x/y方向单位是度垂直方向单位是米直接算出来的坡度会完全错误。gdaldem针对这种情况提供-s参数表示水平与垂直单位比例gdaldem slope 福州市DEM.tif 福州市_slope.tif -s 111120 -p -compute_edges # -s 111120 表示1度约等于111120米-p 输出坡度百分比-compute_edges 处理边缘像元福州在北纬25度至26度之间也可以用111320*cos(lat)估算更精确的比例因子但对30m坡度的结果影响可以忽略。如果DEM本身就是投影坐标系x/y单位就是米不需要-s参数。4.2 用Python批处理生成坡度、山体阴影和晕染更灵活的方式是调用subprocess控制gdaldem配合范围shp批量生成多个产品import subprocess from pathlib import Path dem_path Path(福州市DEM.tif) out_dir Path(products) out_dir.mkdir(exist_okTrue) commands [ [gdaldem, hillshade, str(dem_path), str(out_dir / hillshade.tif), -z, 2, -az, 315, -alt, 45], [gdaldem, slope, str(dem_path), str(out_dir / slope.tif), -s, 111120, -p], [gdaldem, color-relief, str(dem_path), str(out_dir / color.txt), str(out_dir / relief.tif)], ] for cmd in commands: subprocess.run(cmd, checkTrue)-z是垂直拉伸倍数平原地区可以适当调大到2或3让微地形更明显-az是太阳方位角山体阴影常用315度模拟西北方向光-alt是太阳高度角45度最接近人眼阅读习惯。checkTrue保证中间某一步失败时立即停止避免拿到残缺的产品目录。4.3 输出50米间距等高线并叠加shp等高线是根据DEM网格线性内插出的等值线最常用的是gdal_contourgdal_contour -a elev 福州市DEM.tif 福州市_contour.shp -i 50 -off 0 # -a elev 表示把高程值写入elev属性-i 50 表示每50米一条-off 0 表示从0米起算如果文件较大建议先gdal_translate -outsize 10% 10%生成低分辨率预览等参数调好再全分辨率输出。生成的福州市_contour.shp要和福州市范围.shp一起叠加检查看等高线是否被行政边界压实。如果发现等高线在某个区域密集异常大概率是tif里有无效值需要回到第3.3节处理nodata。4.4 成果导出为png便于评审交付处理完的坡度和山体阴影可以用matplotlib输出成图片塞进PPT或评审文档import matplotlib.pyplot as plt import rasterio with rasterio.open(products/hillshade.tif) as hs: hillshade hs.read(1, maskedTrue) plt.imshow(hillshade, cmapgray) plt.axis(off) plt.savefig(products/hillshade.png, dpi150, bbox_inchestight)maskedTrue会忽略nodata区域避免图片边缘出现黑框。注意如果输出图片尺寸太大可以用plt.figure(figsize(8, 8))控制画布但实际打印尺寸主要受dpi影响。5. 裁剪DEM到边界shp时避免黑边的实用参数最后一步经常变成“翻车现场”把福州市范围.shp作为裁剪边界输出的DEM四周出现一圈黑色无值区域。这不是数据问题而是gdalwarp参数没写全。5.1 最小可行命令gdalwarp -cutline 福州市范围.shp -crop_to_cutline -dstnodata -9999 \ 福州市DEM.tif 福州市DEM_clip.tif # -dstnodata 必须设置否则边缘可能被计算成0或nan关键是-crop_to_cutline。如果不加这个参数输出范围仍然是原DEM的外接矩形只是把边界外区域设为nodata看上去还是一整块透明黑。加上之后输出栅格的宽高会紧贴shp的多边形边界边缘干净很多。5.2 进一步压缩体积和统一坐标如果只需要规则矩形研究区用-te指定西北和东南角坐标配合-tr强制像元大小gdalwarp -te 118.8 25.6 120.2 26.4 -tr 0.00027777778 0.00027777778 \ -r bilinear 福州市DEM.tif 福州_crop.tif # -te 依次是 minX minY maxX maxY-r bilinear 用双线性内插重采样如果你手里还有福建省内其他城市的边界shp把福州市范围.shp替换成对应shp流程完全一致。裁剪完先不要删除原始zip等验证dem和shp边界线完全咬合后再清理不迟。本文还有配套的精品资源点击获取
RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进