简介这份福建省漳州市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米能多看出一些微地形不过处理时间和存储差不多要翻倍。下面这个表是我平时对比数据源时常用到的参考数据源分辨率常见格式特点常用场景SRTM30mGeoTIFF全球覆盖山地区域表现稳定区域规划、水文分析、制图底图ASTER GDEM30mGeoTIFF覆盖纬度更高但局部有伪地形大范围初筛需后处理ALOS PALSAR12.5mGeoTIFF细节丰富东南亚地区覆盖好小流域、精细坡度计算标题中的漳州DEM30mGeoTIFF或IMG已按漳州行政边界裁剪附带shp直接进入业务分析如果你拿到的不是标准SRTM分幅文件名而是类似“zhangzhou_dem.tif”这种自定义名第一步建议先查看元数据不要直接扔进ArcGIS或QGIS里出图。因为栅格有效范围、像素深度、NoData值都会影响后续分析。2.2 zip包内文件组成栅格文件与shp要素类的健康度检查一个常见的“含区域范围shp”交付包压缩包里至少会有两类东西。栅格DEM是单波段高程图像可能是GeoTIFF.tif、IMG或GRIDshp边界则是一整套不可拆分的文件集合缺了任何一个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 117EEPSG: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输出格式GTiff3.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_typeogr.wkbPolygon) geojson_layer.CreateFeature(feat.Clone()) geojson_ds None cmd fgdalwarp -cutline {tmp_geojson} -crop_to_cutline -of GTiff zhangzhou_dem_4547.tif {name}_dem.tif subprocess.run(cmd, shellTrue) info gdal.Info(f{name}_dem.tif, statsTrue) # 这里可以解析valid_percent等字段继续处理 os.remove(tmp_geojson)这个脚本把每个行政区边界先转成GeoJSON再交给gdalwarp可以避开中文要素名的编码问题。注意feat.Clone()返回的要素可能和原始图层坐标系统不一致最好先调用geom.TransformTo(layer.GetSpatialRef())如果shp坐标和DEM一致这里不需要额外操作。subprocess.run里的shellTrue在正式的生产脚本里建议换成参数列表形式避免路径空格或特殊字符导致命令失效。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, indexFalse)步长5表示每隔5个像元取一个点相当于150米间隔。适合快速生成300米间隔的勘察点。如果要做更专业的外业布点建议再加一层随机抖动避免点位落在规则格网上导致空间自相关。src.xy(i, j)返回的是栅格像元中心不是像元左上角要和GPS轨迹对齐时应允许10米到20米的平面误差。不同步长对应的成果规模大致如下采样步长间隔约输出点数漳州范围适合用途130m2000万全量高程库390m220万详细地形建模5150m80万Excel可打开10300m20万快速概览、路线初勘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导出的流程几乎覆盖了“福建漳州DEMshp”数据包的所有常规操作。下一次拿到别的省市同类型30米高程数据时你只要替换路径和EPSG代码就能无缝复用。本文还有配套的精品资源点击获取