简介：福建省福州市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 90m | 90m | WGS84地理坐标 | 细节不够 |
| SRTM 1弧秒 | 30m | WGS84地理坐标 | 部分版本有空洞 |
| ALOS AW3D30 | 30m | WGS84/UTM | 下载分块多 |
| ALOS 12.5m | 12.5m | UTM | 数据量大 |
| ASTER GDEM | 30m | WGS84 | 影像伪影多 |
这张表不负责判断你手上这份是哪种，但它告诉你：同样是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 COMPRESS=DEFLATE -co TILED=YES \ dem_30m.tif fuzhou_dem_clip.tif参数解释：-cutline指定裁剪矢量，shp可直接使用；-crop_to_cutline把输出栅格范围收紧到矢量边界的实际范围，而不是只做掩膜；-dstnodata -9999把边界外的像素填充为-9999，防止后续统计把0值当真实高程；-co COMPRESS=DEFLATE启用压缩，栅格体积通常会小一半，-co TILED=YES把内部存储改为瓦片，后续读取局部区域更快。
如果只想按矩形范围快速切一块，不加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 GEOMETRY=AS_XY-t_srs EPSG:4326把输出坐标转成经纬度；-lco GEOMETRY=AS_XY在CSV中额外生成X、Y两列；图层里的MultiPolygon会用重复的点表达，需要在外层处理。如果GDAL版本较老，GEOMETRY=AS_XY可能不生效，用Python读更可靠。
from osgeo import ogr ds = ogr.Open("fuzhou_boundary.shp") layer = ds.GetLayer(0) with open("fuzhou_boundary.txt", "w", encoding="utf-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是WGS84，shp是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 | -compass | 0度方向错误 |
| 山体阴影 | 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], crop=True) data = out_img[0] valid = data[data != src.nodata] print(row["name"], valid.min(), valid.max(), valid.mean())注意mask()的第二个参数需要几何对象列表，传入DataFrame整列会有隐藏坑。crop=True保证输出范围对齐shp边界。src.nodata取的是文件真实填充值，如果文件里写的是-32768而代码里写-9999，统计结果会混入无效像元。需要土方量时，再对valid乘以30×30即可。
5. 用无头脚本校验DEM与shp的完整性和坐标一致性
5.1 用Python脚本检查栅格与shp是否套合
收到zip后，在交付别人之前，我会跑一段快速校验脚本，确认三件事：shp和DEM是否覆盖同一范围，分辨率是否标注错误，nodata是否被当作真实高程。脚本只依赖osgeo：
from 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设为透明即可。
本文还有配套的精品资源，点击获取
