福州市30m DEM数据处理实战：从解包裁剪到坡度等高线生成

简介福建省福州市30米分辨率的DEM数字高程数据包面向GIS从业者、城乡规划人员与地理信息专业学生满足地形起伏分析、坡度坡向计算、洪水淹没模拟、交通选线等应用需求30米分辨率在宏观尺度上兼顾了精度与数据量。压缩包共12个文件核心是福州市DEM.tif以30米×30米栅格表达地表海拔覆盖福州市全境并包含周边部分区域同时附有福州市范围.shp矢量边界配套prj坐标参考文件、dbf属性表、sbn/sbx空间索引以及xml元数据用户在ArcGIS、QGIS等软件中加载后即可直接使用无需自行配准、定义投影或裁剪范围。整个压缩包约35.83MB体积紧凑便于分发和本地存储。目前已有1300人学习下载可用于城乡规划、灾害评估和GIS教学实训包内DEM影像还能生成等高线、制作三维地形场景、提取坡度坡向因子配合shp边界可快速统计不同行政区内的高程分布特征适合科研项目、工程实践与课堂教学数据支撑。1. 福建省福州市DEM数字高程数据zip的解包与数据组织“福建省福州市DEM数字高程数据30m含区域范围shp文件.zip”这个压缩包名字已经把内容拆得很清楚一份标称30米分辨率的DEM栅格外加一份福州市区域范围矢量边界。凡是做地形分析、洪涝淹没模拟、道路选线、甚至光伏选址的人手里大概率都存过同款格式的数据包。真正拿到手以后问题通常不在“有没有数据”而在“数据能不能被直接用”DEM坐标基准是什么shp和栅格是否叠加在同一位置nodata填充值是多少如果这些不确定后面所有坡度、面积统计都是白算。下面按“先看元数据再裁剪对齐然后做应用最后校验交付”的顺序展开。用到的是GDAL、QGIS和Python依赖小、脚本可带走。无论你是第一次打开这个zip还是给甲方做数据预处理都可以照下面的参数先跑一遍再去调整细节。2. DEM数字高程数据30m的分辨率、数据源与坐标系判断2.1 30m像元能看见什么看不见什么DEM的30m指的是每个像元覆盖地面上约30m×30m范围内的高程平均值。福州市区丘陵居多闽江沿岸和鼓山一带高差大。30m网格能说清整个山体的走势能算出坡度是5度还是15度但看不清一条10m宽的冲沟也压不住单栋建筑的高程影响。相比之下热词里常被问到的ALOS 12.5m DEM就能多分辨出一些细小地形。12.5m像元面积是30m像元的约六分之一理论上细节多数据量也接近6倍处理时间不只线性上涨。在实际项目里选哪个分辨率通常不是“越高越好”而是和底图分辨率匹配。如果你只是叠加2020道路shp做宏观选线30m足够如果要细化到一个小型地块的汇水边界最好换12.5m或机载LiDAR。30m的优势是覆盖范围规整、生产标准成熟很多区域范围shp和流域数据都是基于30m DEM切出来的交接成本低。2.2 用gdalinfo和ogrinfo查清坐标系和分辨率拿到数据包里的dem文件第一件事不是拿去算而是用命令行看元数据。以GDAL 3.x环境为例gdalinfo dem_30m.tif输出里重点看三行Size规定栅格宽高像素数Origin规定左上角坐标Pixel Size规定每个像元在地图单位下的宽高Coordinate System is之后的文字是空间参考。如果Pixel Size是0.00027777778说明栅格是地理坐标系这个值约等于1秒也就是赤道附近30m。如果Pixel Size是30说明已经是投影坐标单位是米。再看shp矢量边界ogrinfo -so fuzhou_boundary.shp fuzhou_boundary-so是“summary only”的意思只输出图层概要而不是逐要素打印。重点看Extent空间范围和Layer SRS空间参考。如果shp的范围和dem栅格范围对不上或者一个在经纬度一个在投影米后续步骤必须统一。这里有个容易忽略的细节ogrinfo的输出里shp的Extent单位必须与SRS一致不要只看到数字就以为是经纬度。2.3 没有readme时怎么推断DEM来源压缩包如果只给一个tif和shp没有readme我会用三个指标判断来源。第一看nodata常见SRTM衍生数据经常用-32768或-9999。第二看覆盖范围如果tif的经纬度范围正好覆盖福州市辖区通常说明生产方用市级行政边界做了一次裁剪。第三看边缘30m栅格在海岸带经常出现负值这是因为部分地区的高程基准和海洋深度数据拼接时没有完全对齐。常见公共DEM数据源对比数据源标称分辨率常见坐标系主要弱项SRTM 90m90mWGS84地理坐标细节不够SRTM 1弧秒30mWGS84地理坐标部分版本有空洞ALOS AW3D3030mWGS84/UTM下载分块多ALOS 12.5m12.5mUTM数据量大ASTER GDEM30mWGS84影像伪影多这张表不负责判断你手上这份是哪种但它告诉你同样是30m不同生产机构的精度和空洞分布并不一样。如果之后看到坡度结果里有长条状的低值多半是原始数据源本身的问题不是你的处理流程不对。3. 用GDAL把福州市DEM与shp裁剪到同一坐标系和范围3.1 解压后先核对shp的四个伴生文件shp文件不是单文件。解压后能看到fuzhou_boundary.shp、fuzhou_boundary.shx、fuzhou_boundary.dbf、fuzhou_boundary.prj等。.shx是索引.dbf是属性表.prj是坐标系定义。如果压缩包里只有.shp没有.prj后面在QGIS里会默认用WGS84估算一旦带入真实工程就会偏位。收到数据先跑一句unzip -l 福建省福州市DEM数字高程数据30m含区域范围shp文件.zip列出zip内容清单核对关键文件。如果发现缺.prj或者文件名为乱码可以在QGIS里手动指定坐标系但更稳妥的办法是让提供方补一份因为没有坐标系的空间边界不具备可交换性。shp配套表可以参考扩展名作用缺失影响.shp几何要素主文件不可用.shx几何索引无法在GIS中正常绘制.dbf属性表要素无法带属性.prj坐标系定义无法准确落图.cpg属性编码中文属性乱码3.2 裁剪dem数据gdalwarp按shp精确裁剪最常见的需求是“把福州之外的高程都去掉”这是对dem文件的常规裁剪操作。用gdalwarp按shp边界裁剪gdalwarp -cutline fuzhou_boundary.shp -crop_to_cutline \ -dstnodata -9999 -co COMPRESSDEFLATE -co TILEDYES \ dem_30m.tif fuzhou_dem_clip.tif参数解释-cutline指定裁剪矢量shp可直接使用-crop_to_cutline把输出栅格范围收紧到矢量边界的实际范围而不是只做掩膜-dstnodata -9999把边界外的像素填充为-9999防止后续统计把0值当真实高程-co COMPRESSDEFLATE启用压缩栅格体积通常会小一半-co TILEDYES把内部存储改为瓦片后续读取局部区域更快。如果只想按矩形范围快速切一块不加shp也行gdal_translate -projwin 118.8 26.4 120.2 25.2 \ -a_srs EPSG:4326 dem_30m.tif west_fuzhou.tif-projwin的参数顺序是左上角x、左上角y、右下角x、右下角y。新手容易在25.2和26.4之间颠倒上下来回改gdal_translate不会报错只会输出一张空栅格。注意gdalwarp是按边界形状精确裁剪gdal_translate只能切外接矩形两者用途不同。提示裁剪后用QGIS加载fuzhou_dem_clip.tif如果边界外显示黑色并不是数据坏了而是背景被填了-9999。把渲染的nodata选项勾上即可。3.3 shp转txt导出福州边界坐标用于交接和校验“shp转txt”的本质是把矢量几何转成可读文本便于交给不装GIS的开发或写报告。我一般直接转CSV并把坐标系改成WGS84经纬度ogr2ogr -f CSV fuzhou_boundary_wgs84.csv fuzhou_boundary.shp \ -t_srs EPSG:4326 -lco GEOMETRYAS_XY-t_srs EPSG:4326把输出坐标转成经纬度-lco GEOMETRYAS_XY在CSV中额外生成X、Y两列图层里的MultiPolygon会用重复的点表达需要在外层处理。如果GDAL版本较老GEOMETRYAS_XY可能不生效用Python读更可靠。from osgeo import ogr ds ogr.Open(fuzhou_boundary.shp) layer ds.GetLayer(0) with open(fuzhou_boundary.txt, w, encodingutf-8) as fp: for feature in layer: geom feature.geometry() if geom is not None: fp.write(geom.ExportToWkt() \n)这段代码把每个要素的几何导出成WKT字符串一行一个多边形。WKT比CSV更适合表达复杂边界因为它保留了环的顺序和子多边形边界。实际交接时如果对方只要坐标点再在WKT基础上做字符串拆分即可。3.4 统一坐标系shp与DEM投影不一致时优先做哪个如果gdalinfo显示dem是WGS84shp是CGCS2000高斯投影gdalwarp裁剪时会在内部临时转换但结果容易卡在重投影环节。更可控的做法是先统一到一种坐标再裁剪。判断基准很简单看shp的.prj如果单位是米而栅格的Pixel Size是度就先把shp转成地理坐标再裁剪。gdalwarp也支持-t_srs但建议把裁剪和重投影分两步执行哪个环节出错一目了然。转换shp用一条命令即可ogr2ogr -t_srs EPSG:4326 fuzhou_boundary_wgs84.shp fuzhou_boundary.shp这样生成的shp文件可以用QGIS打开也能作为后续gdalwarp的cutline。福州市常用投影带是120度中央经线的东西带直接选WGS84经纬度能避免带号出错的麻烦。4. 把福州市30m DEM转成坡度、等高线和区域统计结果4.1 坡度坡向与山体阴影三个gdaldem参数表DEM最常见的二次产品是坡度、坡向和山体阴影。GDAL的gdaldem命令一次只能算一种优点是不用装桌面软件。对福州市这种起伏较大的区域30m栅格算坡度得到的角度能反映沟谷密度和断裂走向。常用参数gdaldem slope dem_30m.tif fuzhou_slope.tif -p -s 111120 gdaldem aspect dem_30m.tif fuzhou_aspect.tif -compass gdaldem hillshade dem_30m.tif fuzhou_hillshade.tif -az 315 -alt 45第一条命令里-p要求输出坡度百分比如果去掉-p则输出0到90的度数值不少人的坡度直方图和自己预期对不上往往就是百分比和度数的差别。-s 111120是垂直放大倍数只在地理坐标系下生效含义是“1度约等于111120米”投影坐标系下应该删除这个参数。第二条-compass把坡向输出为北向0度顺时针到360度如果不加0度指东和常用的方位表示不一致。第三条的-az 315设置光源来自西北-alt 45设置太阳高度角45度适合观察福州西北向山体形态。参数调整可以参考输出命令必调参数常见误用坡度gdaldem slope-p / -s投影坐标系下仍留-s坡向gdaldem aspect-compass0度方向错误山体阴影gdaldem hillshade-az -alt高度角过高没有阴影4.2 DSM和DEM的关系什么时候需要从DSM生成DEM很多城市级高程包其实是DSM表面高程包含了建筑和树冠。DSM直接算坡度屋顶和道路连接处会形成断层。热词“dsm生成dem”描述的正是这种前处理把DSM中的地物盖层去掉还原裸露地表高程。概念上DEM是地形表面DSM是地物表面两者之差是冠层高度。处理时常用移动窗口低通滤波或者用LiDAR点云做地面分类。对30m格网来说直接在DSM上做低通滤波会把真正的山脊线一起抹平所以我会先测一次滤波前后的坡度差异如果差异集中在城区再决定是否做DSM转DEM。4.3 用gdal_contour生成shp文件等高线直接交付CAD等高线是dem文件的矢量产品也是“生成shp文件”最常见的场景。命令gdal_contour -a elev -i 20 fuzhou_dem_clip.tif fuzhou_contour_20m.shp-a elev表示在输出shp属性表里加一个叫elev的字段记录高程值-i 20表示每20米画一条等高线。福州市区闽江两岸高差小20米一条会显得稀疏到了鼓山一带又过密所以绘制大区域时常用分段间距先按30米生成再在关键区域单独加密到5米。gdal_contour输出的shp线条数量在城区可能达到几万条建议生成后用ogr2ogr -simplify做线条简化再交付。4.4 用shp对DEM做分区统计各区县平均高程这一步回答“福州市某一片区平均高程是多少”这类问题。用rasterio的mask函数先按shp边界裁剪再对像元值做统计。import geopandas as gpd import rasterio from rasterio.mask import mask import numpy as np # 如果你的shp是区县边界循环计算每个区域 city gpd.read_file(fuzhou_counties.shp) with rasterio.open(fuzhou_dem_clip.tif) as src: for _, row in city.iterrows(): geom row.geometry out_img, out_transform mask(src, [geom], cropTrue) data out_img[0] valid data[data ! src.nodata] print(row[name], valid.min(), valid.max(), valid.mean())注意mask()的第二个参数需要几何对象列表传入DataFrame整列会有隐藏坑。cropTrue保证输出范围对齐shp边界。src.nodata取的是文件真实填充值如果文件里写的是-32768而代码里写-9999统计结果会混入无效像元。需要土方量时再对valid乘以30×30即可。5. 用无头脚本校验DEM与shp的完整性和坐标一致性5.1 用Python脚本检查栅格与shp是否套合收到zip后在交付别人之前我会跑一段快速校验脚本确认三件事shp和DEM是否覆盖同一范围分辨率是否标注错误nodata是否被当作真实高程。脚本只依赖osgeofrom osgeo import gdal, ogr dem gdal.Open(fuzhou_dem_clip.tif) gt dem.GetGeoTransform() cols, rows dem.RasterXSize, dem.RasterYSize width cols * abs(gt[1]) height rows * abs(gt[5]) shp ogr.Open(fuzhou_boundary.shp) extent shp.GetLayer(0).GetExtent() extent_width extent[1] - extent[0] extent_height extent[3] - extent[2] x_ok abs(width - extent_width) / extent_width 0.01 y_ok abs(height - extent_height) / extent_height 0.01 print(栅格范围与shp匹配, x_ok and y_ok)这段代码从GetGeoTransform()取出像元宽高乘以栅格行列数得到DEM的总跨度再与shp的Extent()比较。相对误差小于1%就认为匹配避免直接在经纬度和米之间比较。如果裁剪前有大量余量宽高会比shp范围大出几百倍打印结果立刻能发现问题。5.2 大范围任务时用渔网分割shp降内存福州一个市的范围单块30m栅格通常只有几百万像元直接处理没问题。但如果把整个福建省的DEM一起灌进来或者要做逐格网统计内存会先爆。热词“渔网分割shp”就是这个场景的解法在QGIS的Vector creation – Create grid工具里生成矩形渔网再用Clip把shp按格网切片每片单独算坡度最后用gdalbuildvrt拼回一个VRT。批量循环里一次只打开一个切片文件峰值内存可以降到原来的几十分之一。5.3 输出文件的最终验证数据交接的最后一公里经常是“文件打开后显示黑色”原因多半是nodata没有做透明处理。确认方法gdalinfo -stats fuzhou_dem_clip.tif-stats会计算最大值最小值。如果最小值和nodata一致说明栅格沿边界仍有留空区在QGIS中使用“Singleband pseudocolor”渲染时选择“Min/Max”并把nodata设为透明即可。本文还有配套的精品资源点击获取