简介：这份福建省漳州市30米分辨率DEM数字高程数据包，面向GIS学习者、城乡规划与地质灾害评估人员，可用于地形分析、坡度坡向提取、洪水模拟等场景。压缩包共12个文件，大小约35.1MB，核心为漳州市DEM.tif高程栅格，配套行政范围Shapefile（含shp、dbf、prj、sbn/sbx/shx索引）及tfw坐标参考、xml元数据等，可确保在ArcGIS、QGIS中正确加载与配准。已有478人学习下载。通过该数据可获取漳州市完整地形模型，结合边界文件快速裁剪出研究区域，进而计算坡度、坡向、山脊山谷线，为城市选址、道路选线、生态保护区划分提供基础数据支撑。
1. 拿到“福建省漳州市DEM数字高程数据30m（含区域范围shp文件）.zip”后要做什么？
这一个数据包看起来只是某次项目交付中的常见产物，但它实际覆盖了地理数据生产里一条非常完整的链路：以覆盖漳州市范围的30米分辨率数字高程模型（DEM）栅格为主数据，再附带一个用于定位和切边界的行政边界矢量文件（shp）。解压后，你通常会同时看到TIF/IMG格式的高程栅格，以及一组由.dbf、.shp、.shx、.prj组成的矢量文件。对于做规划、国土、交通、通信覆盖或户外选址的工程师，这份数据能直接用于坡度分析、可视域计算、流域提取和三维地形底图；对刚接触GIS或遥感的人来说，它也是难得的“有真边界、有真地形”的学习素材。下面我按自己处理这类数据的常规顺序，从数据体检、坐标归一、按shp裁剪、批量提取到成果验证一步步展开。
2. 弄清楚30米DEM和shp文件的家底：格式、命名与坐标参考
2.1 30米分辨率DEM的定位：SRTM、ASTER与ALOS的差异
标题里强调“30m”，意味着栅格每个像元对应地面约30米见方。这个尺度在国土空间分析和区域规划里是黄金比例：比90米数据能看清山脊线和沟谷变化，比12.5米数据（常见来源是ALOS PALSAR）体量小一个数量级，普通笔记本就能流畅处理。通常30米DEM源是NASA的SRTM（覆盖全球北纬60°到南纬56°）和ASTER GDEM（覆盖更广但噪声略高）。国内很多数据集也基于SRTM做了填补、重投影和按行政区裁剪。和12.5米ALOS数据相比，30米在平缓平原上细节差异不大，但在漳州这种西北多山、东南沿海丘陵的地形里，做路径规划或基站选址时12.5米能多看出一些微地形，不过处理时间和存储差不多要翻倍。下面这个表是我平时对比数据源时常用到的参考：
| 数据源 | 分辨率 | 常见格式 | 特点 | 常用场景 |
|---|---|---|---|---|
| SRTM | 30m | GeoTIFF | 全球覆盖，山地区域表现稳定 | 区域规划、水文分析、制图底图 |
| ASTER GDEM | 30m | GeoTIFF | 覆盖纬度更高，但局部有伪地形 | 大范围初筛，需后处理 |
| ALOS PALSAR | 12.5m | GeoTIFF | 细节丰富，东南亚地区覆盖好 | 小流域、精细坡度计算 |
| 标题中的漳州DEM | 30m | GeoTIFF或IMG | 已按漳州行政边界裁剪，附带shp | 直接进入业务分析 |
如果你拿到的不是标准SRTM分幅文件名，而是类似“zhangzhou_dem.tif”这种自定义名，第一步建议先查看元数据，不要直接扔进ArcGIS或QGIS里出图。因为栅格有效范围、像素深度、NoData值都会影响后续分析。
2.2 zip包内文件组成：栅格文件与shp要素类的健康度检查
一个常见的“含区域范围shp”交付包，压缩包里至少会有两类东西。栅格DEM是单波段高程图像，可能是GeoTIFF（.tif）、IMG或GRID；shp边界则是一整套不可拆分的文件集合，缺了任何一个，ArcGIS或GDAL都可能识别失败。典型的shp家族包括.shp（几何）、.shx（索引）、.dbf（属性）、.prj（投影信息），如果是UTF-8编码可能还有.cpg。拿到后先不要急于单独复制一个.shp文件走人，应该整目录解压，否则会遇到“无法打开要素类”的报错。
在终端里用GDAL自带工具做体检是最快的。Linux或macOS下安装了gdal后，可以执行：
unzip 福建省漳州市DEM数字高程数据30m（含区域范围shp文件）.zip -d zhangzhou_dem gdalinfo zhangzhou_dem/*.tif | head -30 ogrinfo -so -al zhangzhou_dem/*.shp如果压缩包里有多个tif或shp，上面带通配符的命令可能无法正确匹配，需要逐个指定文件名。gdalinfo输出里的Size is和Origin告诉你栅格行列数、起始坐标；Coordinate System is告诉你投影，比如WGS 84或CGCS2000 / 3-degree Gauss-Kruger zone 38。这些信息直接决定后面裁剪和叠加时的处理策略。ogrinfo -so -al会输出shp的要素数量、图层坐标范围、属性字段名，例如看到CNTY_CODE和NAME字段，就能知道边界是区县级还是镇级。
2.3 检查坐标系：不可跳过的对齐步骤
很多拿到数据的人直接在QGIS里把tif和shp拖进去，看到两者位置对不上，第一反应是数据坏了。其实大部分原因是栅格是WGS84经纬度，shp是CGCS2000投影坐标；又或者一个是EPSG:4490，另一个是EPSG:4547之类。判断方式很简单，用ogrinfo看shp的prj信息，用gdalinfo看tif的投影。如果都是WGS84，位置偏差只是可视化上的拉伸；如果范围差很多，就必须先统一坐标系。
统一坐标系要区分“数据本身用哪种投影”和“计算时用哪种投影”。算坡度可以用经纬度直接算，但算面积、距离，以及把DEM和shp做叠加裁剪，我建议统一到Albers等积圆锥投影或CGCS2000高斯投影。福建省内常用CGCS2000 / Gauss-Kruger CM 117E（EPSG:4547）或Web墨卡托（EPSG:3857）作为中间输出。注意Web墨卡托会造成面积失真，漳州纬度在24度左右，南北变形不大，但做严谨的坡度面积统计时请用EPSG:4547。
3. 用GDAL系列命令把DEM和shp对齐、裁剪并派生坡度坡向
3.1 栅格投影转换：gdalwarp的参数不必每次都从零记
当你发现DEM和shp投影不一致，直接覆盖写一份投影一致的临时文件，比每次调用都实时转换更稳妥。最常见做法是用gdalwarp把栅格重投影到与shp相同的坐标系。先通过ogrinfo -al -so拿到shp的EPSG代码，比如是4547，然后执行：
gdalwarp -t_srs EPSG:4547 -r bilinear -of GTiff zhangzhou_dem.tif zhangzhou_dem_4547.tif-t_srs定义输出投影，-r bilinear是重采样算法。注意：重采样算法不是随便选的。DEM是连续高程表面，双线性（bilinear）或三次卷积（cubic）更适合保持地形平滑；如果选最邻近（nearest），会出现台阶状纹理。-of GTiff指定输出为GeoTIFF。命令执行完后再次gdalinfo确认Pixel Size从0.0003度左右变成30米级别。
另外，重投影后最好再执行一次gdalinfo，对比输出文件和shp的范围边界。有些时候由于DEM在覆盖范围边缘存在拉伸或缺失，warp后会生成一小块黑边。这时用-dstnodata显式指定无效值，并且在后续所有统计命令里统一使用该值，能让结果干净不少。我在处理漳州这类丘陵地形时，通常把无效值设为-9999，而原始SRTM的无效值可能是-32768，两者混用会直接污染统计结果。
3.2 按shp边界裁剪：两条命令和一条黄金参数
裁剪是这类数据包最核心的操作。既然压缩包里已经带了漳州市的范围shp，那就不需要自己从全国矢量数据里抠了。裁剪前先确认shp和DEM是否同一投影；如果刚才已经做了-t_srs步骤，现在可以直接执行：
gdalwarp -cutline zhangzhou.shp -crop_to_cutline -of GTiff zhangzhou_dem_4547.tif zhz_dem_clip.tif-cutline指定边界矢量，-crop_to_cutline是黄金参数：它让输出栅格的范围收缩到shp的最小外接矩形，同时把边界外的像元设置为无效值。不加这个参数，输出范围会沿用DEM原始范围，只是把外部像元遮蔽，文件体积一点没小。对整包数据来说，裁剪后可能从几百万像元缩到几十万，后续计算会快很多。另外需要注意，-cutline搭配-crop_to_cutline时，GDAL只支持ESRI Shapefile、GeoJSON等若干矢量驱动；如果shp路径里有中文，偶尔会报错，建议先把路径换成英文。
裁剪后检查一下无效值。用gdalinfo -stats zhz_dem_clip.tif查看STATISTICS_VALID_PERCENT，正常应该在98%以上。因为行政边界是锯齿状，完全贴边的像元不一定都是有效值，少量边界像元成为NoData属于正常。但如果有效百分比低于90%，说明shp范围与DEM重叠区域太少，要回去检查投影或边界文件是否拿错。下面这张表列出常用参数：
| 参数 | 作用 | 建议值 |
|---|---|---|
| -cutline | 指定裁剪矢量 | 优先用GeoJSON避免中文路径问题 |
| -crop_to_cutline | 以矢量范围输出栅格 | 必加 |
| -dstnodata | 设置无效值 | 统一为-9999 |
| -r | 重采样算法 | 连续表面用bilinear或cubic |
| -of | 输出格式 | GTiff |
3.3 一次算出坡度、坡向和山体阴影：不装ArcGIS也能出图
拿到裁剪后的DEM，做地形分析最快的方式是用GDAL自带的gdaldem，它支持hillshade、slope、aspect、color-relief等模式。以漳州西部山区为例，要给规划报告配三张图，直接执行三条命令：
gdaldem slope zhz_dem_clip.tif zhz_slope.tif -p -s 1 -of GTiff gdaldem aspect zhz_dem_clip.tif zhz_aspect.tif -zero_for_flat -of GTiff gdaldem hillshade zhz_dem_clip.tif zhz_hillshade.tif -z 1.0 -az 315 -alt 45 -of GTiff-p让坡度输出为百分制而非度制，做地质灾害评估或坡度分级时更直观；-s是垂直比例因子，平面坐标加米制高程通常设置为1；如果DEM是经纬度坐标，这个值要设为约111320（每度长度约111公里），否则坡度会被严重放大。-az和-alt分别是山体阴影的太阳方位角和高度角，默认315度和45度在平原地区没问题，山区建议先看地形走向再调。如果发现输出tif全是黑色，多半是输入DEM的NoData值没有被正确识别，重投影时用-dstnodata -9999显式指定一下即可。
这里有个常见误区：gdaldem的-s参数和三维显示里的垂直夸张因子不是一回事。三维场景里为了视觉效果经常把高程拉伸到2倍或5倍，但坡度、坡向计算要求在水平和垂直方向属于同一度量制。用经纬度坐标直接算坡度时，如果不乘以111320，在福建山区算出来的坡度能差出接近一倍，所以在裁剪前完成投影转换很有必要。
4. 结合shp做批量区域提取、高程导出和等高线叠加
4.1 用Python循环处理多个区县shp：拆分、裁剪、统计一步完成
漳州下辖多个区县，当压缩包里给的shp是全市范围时，可以用自己的区县边界做更细的统计。常见做法是：循环读shp中的每个要素，逐个调用gdalwarp，再计算高程均值和分位数。用Python的subprocess最省事，不必为每个功能单独调C++接口：
import subprocess import os from osgeo import gdal, ogr ds = ogr.Open('zhangzhou.shp') layer = ds.GetLayer() for feat in layer: name = feat.GetField('NAME') geom = feat.GetGeometryRef() tmp_geojson = f'{name}.geojson' # 把单个要素导出为GeoJSON临时文件 geojson_ds = ogr.GetDriverByName('GeoJSON').CreateDataSource(tmp_geojson) geojson_layer = geojson_ds.CreateLayer('boundary', geom_type=ogr.wkbPolygon) geojson_layer.CreateFeature(feat.Clone()) geojson_ds = None cmd = f'gdalwarp -cutline {tmp_geojson} -crop_to_cutline -of GTiff zhangzhou_dem_4547.tif {name}_dem.tif' subprocess.run(cmd, shell=True) info = gdal.Info(f'{name}_dem.tif', stats=True) # 这里可以解析valid_percent等字段继续处理 os.remove(tmp_geojson)这个脚本把每个行政区边界先转成GeoJSON，再交给gdalwarp，可以避开中文要素名的编码问题。注意feat.Clone()返回的要素可能和原始图层坐标系统不一致，最好先调用geom.TransformTo(layer.GetSpatialRef())；如果shp坐标和DEM一致，这里不需要额外操作。subprocess.run里的shell=True在正式的生产脚本里建议换成参数列表形式，避免路径空格或特殊字符导致命令失效。
4.2 将DEM高程点导出为txt/csv：对接勘察设备与Excel统计
现场踏勘或通信覆盖仿真时，需要把栅格高程变成离散点。这个问题和“shp转txt”是同一类需求：把DEM按固定间隔采样，生成带坐标和高程的文本文件。最简单的路径是先用gdal_translate输出XYZ文本：
gdal_translate -of XYZ zhz_dem_clip.tif zhz_elevation.xyz head -10 zhz_elevation.xyz这样生成的xyz文件每行是“经度 纬度 高程”，但它是把所有像元全覆盖输出。漳州全境30米分辨率会得到约2000万行，直接给Excel会卡死。建议先用gdalwarp重采样到100米或200米，再用上面的命令，或者直接用Python的rasterio读取后按步长抽样：
import rasterio import pandas as pd rows = [] with rasterio.open('zhz_dem_clip.tif') as src: data = src.read(1) for i in range(0, data.shape[0], 5): for j in range(0, data.shape[1], 5): if data[i, j] > -9999: x, y = src.xy(i, j) rows.append((x, y, data[i, j])) df = pd.DataFrame(rows, columns=['lon', 'lat', 'elevation']) df.to_csv('zhz_sample.csv', index=False)步长5表示每隔5个像元取一个点，相当于150米间隔。适合快速生成300米间隔的勘察点。如果要做更专业的外业布点，建议再加一层随机抖动，避免点位落在规则格网上导致空间自相关。src.xy(i, j)返回的是栅格像元中心，不是像元左上角；要和GPS轨迹对齐时，应允许10米到20米的平面误差。不同步长对应的成果规模大致如下：
| 采样步长 | 间隔约 | 输出点数（漳州范围） | 适合用途 |
|---|---|---|---|
| 1 | 30m | 2000万 | 全量高程库 |
| 3 | 90m | 220万 | 详细地形建模 |
| 5 | 150m | 80万 | Excel可打开 |
| 10 | 300m | 20万 | 快速概览、路线初勘 |
4.3 等高线提取与shp叠加：出地图的常见顺序
很多人的目标是得到一幅带地形层和边界层的地图。做法是从裁剪后的DEM用gdal_contour提取等高线，再和shp边界叠在QGIS或ArcGIS里制图排版：
gdal_contour -a ELEV -i 50 -nln contour zhz_dem_clip.tif zhz_contour.shp-a ELEV指定把高程值写入线要素的字段名，-i 50表示每隔50米生成一条等高线。漳州从海边的0米到西部山区的约1000米，50米间隔大约能出20条主线，放在图例里不拥挤。如果要让等高线更平滑，可以先用gdalwarp -r bilinear降低分辨率到60米或90米，但我不建议对DEM做平滑后再提取等高线，因为会破坏真实地形的微起伏。更专业的做法是使用GRASS的r.contour生成带标注的矢量线，不过生产环境里gdal_contour已经足够。
等高线shp生成后，用ogr2ogr整理编码，或者在QGIS里设置标注列为ELEV。叠加边界时，因为等高线是从裁剪后的DEM提取，而DEM边界是shp的平滑版，所以等高线和边界之间会出现几个像元左右宽的空白带，这不是错误。如果坚持让等高线严格截止到边界，可以用ogr2ogr -clipsrc zhangzhou.shp再裁剪一次。
5. 进阶用法和验证：从12.5米换数据源到shp边界线的细节处理
5.1 用ALOS 12.5米数据验证30米DEM的可信区间
关于“12.5米dem下载”的检索热度很高，原因是30米DEM在山谷地形里会把窄谷和细小沟壑压成一片平地。我在拿到这份漳州数据后，通常选漳州西北部的南靖县或华安县一小块范围，下载ALOS 12.5米数据做交叉验证。注意这不是要替代30米，而是看同一位置的高程差范围。做法：把ALOS裁剪到同样shp范围，重采样到30米，然后用栅格计算器计算差值（两组高程相减），查看平均值和标准差。标准差小于10米说明30米数据在区域内基本可靠；若出现几处带状负差，则可能是原DEM的河网区域被过度平滑，后续做淹没分析时需要倍加小心。
5.2 用像素统计快速验收DEM是否经历过“坏点填充”
判断DEM要不要做二次修复，不需要打开软件，Python结合numpy就能给出结论。用rasterio读入数据后，统计高程直方图，检测是否有异常条带，比如大量像元完全相同、传感器坏线留下的矩形空洞，同时观察NoData分布：
import numpy as np import rasterio with rasterio.open('zhz_dem_clip.tif') as src: arr = src.read(1).astype(np.float32) nodata = src.nodata if nodata is not None: arr[arr == nodata] = np.nan valid = arr[~np.isnan(arr)] print('有效像元数:', len(valid)) print('高程范围:', np.nanmin(valid), np.nanmax(valid)) print('异常高值占比:', np.sum(valid > np.nanpercentile(valid, 99.7)) / len(valid))如果异常高值占比超过5%，建议先对DEM做中值滤波再使用。不过要注意，这里的异常值检测只是统计学意义上的初步筛选，不代表地理学合理性。例如漳州沿海有海拔为负的滩涂，出现-10米以内的负高程是合理的；但内陆山区出现-50米，大概率是空洞被内插成了错误值。
5.3 shp只保留外边界线：从完整面要素提取单条边界
最后一个常用技巧：很多人拿到的是整个漳州面的shp，但想做剖面图或写报告时只需要最外边界一条线。QGIS里可以用“矢量几何工具-边界”，命令行则用GDAL的SQLite方言：
ogr2ogr -dialect sqlite -sql "SELECT ST_ExteriorRing(geometry) AS geometry FROM zhangzhou" zhangzhou_boundary.shp zhangzhou.shp注意，如果原始shp包含多个不相连的面要素（比如东山岛），直接使用ST_ExteriorRing会输出多条外环线。正确做法是先对几何做ST_Union合并，再取外环。因为漳州主体是一个连通行政区，可以按上面这条命令处理；如果要保留岛屿，就换用ST_Boundary。生成后的边界线和等高线叠加时如果出现明显断头，说明DEM裁边处像元偏碎，建议先做一次3x3的焦点统计，让边界线处的像元更平滑。
这套从gdalinfo体检到ogr2ogr导出的流程，几乎覆盖了“福建漳州DEM+shp”数据包的所有常规操作。下一次拿到别的省市同类型30米高程数据时，你只要替换路径和EPSG代码，就能无缝复用。
本文还有配套的精品资源，点击获取
