简介这份福建省南平市30米分辨率DEM数字高程数据包面向GIS学习者、城乡规划与地理分析从业者适用于地形起伏分析、环境规划、城市设计、灾害评估等场景。包内含南平市全域及周边范围的高程栅格数据并附行政范围Shapefile边界文件便于直接裁剪、制图与叠加分析。压缩包共12个文件主要以tif高程数据、shp矢量边界及配套dbf属性、prj投影、sbn/sbx/shx空间索引等GIS通用格式组成另有tfw/ovr/xml等辅助元数据文件整体约97.12MB结构完整可直接在ArcGIS、QGIS等平台加载。已有298人学习适合需要快速获取标准化区域DEM数据、练习地形分析或开展项目研究的用户使用。1. 拿到“福建省南平市DEM数字高程数据30m含区域范围shp文件.zip”后第一件事不是解压第一次拿到福建省南平市DEM数字高程数据30m含区域范围shp文件.zip这个压缩包的人通常会把它当成一张“带边界的高程图”解压、用GIS打开、看一眼颜色然后就直接开始做分析。我的建议是先别急着解压。这份数据真正值钱的是里面几十万个浮点像元每个像元记录的是南平市某个约30m×30m的格网内的平均地面高程单位一般是米有了它你才能继续做坡度、坡向、山体阴影和汇水区划分而shp文件负责告诉所有工具“分析范围到哪里为止”。这篇笔记按我处理这类数据包的先后顺序来写怎么确认数据可信、怎么裁剪投影、哪里最容易翻车以及最后怎么验证结果。整个过程不复杂但每一步选错参数后面都要花几倍时间返工。2. 先弄清30m DEM和shp的约定再考虑怎么处理很多问题不是因为DEM本身坏掉了而是因为你在还不知道“坐标系是什么、NoData是多少”的情况下就让工具开始干活。一个常见场景打开DEM后发现用gdalinfo看到的分辨率是两个小数说明数据还是经纬度坐标系而shp可能也是经纬度也可能是投影坐标系。如果你没注意到这个差异后面裁剪出来的结果甚至可能旋转了90度。所以处理这套数据的第一步不是改数据而是把它的“身份信息”读出来。2.1 30m分辨率不等于精度30米30m分辨率说的是每个像元在地面上覆盖的面积约为900平方米而不是说它的高程误差是30米。在高程起伏较大的区域一个30m像元实际上是把一个坡段平均成了一块平面所以山脊线、沟谷的细节会被明显平滑。这也是为什么很多人拿30m数据做小流域分析时总觉得汇水线位置和野外对不上——不是数据错了是像元尺度限制了细节。另一个要注意的点是DEM的像元值是可以为负的比如某些洼地或特殊地形的高程记录千万不要把负值一律当成异常值。这份压缩包里的DEM常见的是单个GeoTIFF文件也可能是按图幅分成多个tif。shp则是一组侧车文件组成的至少包含.shp、.shx、.dbf三个文件如果能看到.prj就说明带了坐标信息。解压后可以先看一眼文件清单文件类型常见扩展名作用DEM栅格.tif / .img / .grd存储30m高程值单波段浮点边界矢量.shp区域范围面要素shp索引/属性.shx / .dbf矢量空间索引与属性表投影描述.prj记录矢量坐标系辅助世界文件.tfw部分tif附带便于其他软件识别如果解压后只有一个.shp而没有.shx和.dbf那你可能需要先用平台自带的修复工具重建索引否则后面读取时容易报错。我的习惯是先复制原始文件到工作目录再做任何修改不要把源数据留在zip里反复解压。2.2 用Python打开DEM先回答四个问题我一般会用rasterio和geopandas把DEM和shp同时读进来看四件事坐标系、像元大小、NoData、有效范围。下面这个代码可以当作起步模板import rasterio import geopandas as gpd dem_path nanping_dem_30m.tif shp_path nanping_boundary.shp with rasterio.open(dem_path) as src: print(坐标系:, src.crs) print(像元大小(度或米):, src.res) print(影像范围:, src.bounds) print(NoData:, src.nodata) print(数据类型:, src.dtypes[0]) gdf gpd.read_file(shp_path) print(shp坐标系:, gdf.crs) print(shp要素数量:, len(gdf)) print(shp总范围:, gdf.total_bounds)这段代码是这套流程的“开机自检”。src.crs告诉你数据用的是什么坐标参考src.res如果返回两个一样的数说明是正方形像元如果两个数差很多可能是投影转换后的结果没做对齐后面需要重采样。src.nodata是重中之重很多DEM用 -9999 表示无效值但有的数据源用 -3.402823e38 这类浮点极值读出来才能对症处理。gdf.crs如果显示为 None说明shp缺少.prj那就要先解决投影信息而不是直接拿去裁剪。2.3 当shp边界和DEM边缘对不齐的根源经常有人发现明明都是从“南平市”来的数据shp边界却超出了DEM的覆盖范围或者反过来。原因通常有三个坐标系不一致、重投影后像元范围变化、shp边界版本和DEM生产年代不同。坐标系不一致时最稳妥的办法不是手动输入一个EPSG而是直接读取shp的.prj文件把DEM重投影到它的坐标系。如果shp没有.prj就根据坐标数值量级判断上百度的数就是经纬度六位数的数就是投影坐标米制。重投影可以用gdalwarp但要注意重投影本身会改变像元大小和DEM数值的平滑度所以能避免就避免。如果只是范围边缘少量错位可以在裁剪时给shp做一个小缓冲让边界不完全贴合后面分析再按严格边界裁剪避免边缘像元被切掉一半。3. 预处理从ZIP到能用于坡度、水文分析的DEM拿到解压后的数据我推荐的顺序是先拼接分幅再统一投影然后裁剪最后填洼。这个顺序的好处是每一步都在处理“看得见”的问题分幅拼接解决碎片化投影统一解决坐标系混乱裁剪解决范围填洼解决水文分析的前提条件。不要一开始就填洼因为填洼会改变高程值边界外如果还有NoData填洼结果会把边界附近的地形弄出奇怪的帽子。3.1 解压、拼接、裁剪的最小命令在Linux或者macOS下解压带中文文件名压缩包时部分终端会出现乱码。最简单的方式是用系统自带的图形界面解压或者用Python来做避免编码问题。解压后如果发现多个tif文件先用gdal_merge.py拼成一张unzip 福建省南平市DEM数字高程数据30m含区域范围shp文件.zip -d nanping_dem cd nanping_dem gdal_merge.py -o merged_dem.tif -n -9999 -a_nodata -9999 -ot Float32 tile_1.tif tile_2.tif tile_3.tif-n -9999表示输入数据里值为 -9999 的像素视为无效像素-a_nodata -9999则把输出文件的NoData重新设置为 -9999-ot Float32保证高程值不降级成整数。如果分幅文件之间重叠或者有接边横纹这一步不会消除接边误差所以拼接后最好生成一张山体阴影图检查接边处有没有断层。之后裁剪到shp边界gdalwarp -cutline nanping_boundary.shp -crop_to_cutline \ -dstnodata -9999 -r bilinear -tr 30 30 \ merged_dem.tif nanping_dem_cropped.tif参数含义-cutline告诉gdal用哪个矢量文件裁剪-crop_to_cutline让输出范围贴合边界外框-dstnodata强制输出无效值为 -9999-r bilinear指定重采样算法-tr 30 30把输出像元大小固定为30米。如果你本身就是30m数据也建议加上-tr因为投影转换过程中原始像元网格往往会被旋转偏移明确指定大小可以让后续处理更可控。3.2 重采样参数不要所有场景都用nearest很多人会在gdalwarp里习惯性地用-r near这个参数用于离散类别数据没错但用在DEM上会让高程值出现明显的块状阶梯。对连续地表数据常见做法是双线性或三次卷积。我一般按下面这个表选择需求推荐重采样算法原因坡度、山体阴影、水文分析bilinear / cubic保持表面连续减少锯齿从30m降到90maverage聚合时取均值最接近真实起伏只要确定像元对应位置bilinear速度快且稳定分类数据、专题图near保持类别编号不被打乱如果你把30m重采样到60m我个人更倾向-r average而不是cubic。因为cubic会在陡峭地形中产生轻微过冲导致山顶高程被抬高虽然视觉上更平滑但水文分析里这种系统性偏差很难排查。重采样完之后再次用rasterio检查src.res和src.nodata确保没被悄悄改掉。3.3 填洼做水文分析前最重要的一步填洼几乎是所有DEM水文分析的前置操作。因为DEM或者天然地形中存在凹陷水流到那里就停了无法继续往下游汇流。常见做法是使用专门的地形分析工具集填洼工具会根据周围的像元把洼地“抬平”让水流可以顺利通过。下面这行命令是可复制的版本whitebox_tools -r FillDepressions -i nanping_dem_cropped.tif \ -o nanping_dem_fill.tif --flat_increment 0.01-i是输入DEM-o是输出文件--flat_increment表示平坦区域加入的微小高程增量。这个增量很重要如果设为0所有平坦像元高程完全相同后续D8流向计算会在平地上出现大量随机方向。0.01米是一个很小的值不会干扰整体地形趋势但足以让流向算法有稳定的方向。填洼会改变原始高程所以永远要在复制出来的文件上做原始DEM保留用于坡度计算。我这里踩过一次坑直接在原始DEM上填洼后做坡度结果把山谷边坡全部抹平了所有水文分析结果都偏缓。4. 避坑处理南平市30m DEM时常见的5个翻车现场以下五类问题是我在多个DEM数据处理场景中反复遇到的现象很典型原因各不相同但都可排查。4.1 面积统计和shp差很多甚至偏差几百平方公里现象在GIS里用shp算面积再统计DEM有效像元数量乘以900平方米两者对不上。原因shp本身是经纬度坐标系直接算出来的“面积”单位是度不是平方米或者DEM被重投影后像元面积变化但统计时仍用900平方米。解决先把shp和DEM都转到同一投影坐标系再算面积。对于南平市这样的范围投影带选择要参考数据自带的中央经线信息不要随便选一个通用横轴墨卡托。转换后用下面的方式核对import geopandas as gpd gdf gpd.read_file(nanping_boundary.shp) shp_area gdf.to_crs(gdf.estimate_utm_crs()).area.sum() print(shp面积(平方米):, shp_area)estimate_utm_crs是geopandas提供的一个估算方法会基于shp自身范围选一个合适UTM带。这样算出的面积才是可比较的。4.2 NoData值混进高程统计坡度图一片黑现象高程统计最小值是 -9999最大高程几百米平均高程却被拉低到负值坡度图上边界一圈像被烧焦了一样。原因裁剪或拼接后没有设置NoData工具把 -9999 当作真实高程参与邻域计算。解决每次输出都显式指定dstnodata读取时也要用mask。建议用代码固定处理import rasterio import numpy as np with rasterio.open(nanping_dem_cropped.tif) as src: dem src.read(1, maskedTrue) print(有效高程最小值:, np.nanmin(dem)) print(有效高程最大值:, np.nanmax(dem))read(1, maskedTrue)会自动把NoData像元标记为mask后续统计用np.nanmin就不会被 -9999 污染。如果你已经生成出了一个带异常值的图层最快补救是用gdal_edit.py -a_nodata -9999给文件补上NoData声明而不是手动去改像元值。4.3 裁剪边界锯齿明显但不是数据坏了现象按shp边界裁剪后矢量边界上的阶梯状锯齿非常显眼放大看山脊线也变成了不规则的方块边缘。原因栅格裁剪本质上是用矩形像元去逼近矢量边界任何30m数据都会有这种锯齿。解决如果你只是拿DEM做区域统计锯齿不影响结果不用管如果要做边缘平滑的可视化产品可以先对shp做一个小的缓冲再裁剪让边界留出几米余量然后用叠加方式隐藏锯齿。注意不要为了平滑边缘对DEM做低通滤波那样会把真实的高程起伏也抹掉。4.4 shp缺少.prj文件叠加时整体偏出几千米现象gdf.crs输出为None把shp叠加到DEM上时边界跑到了数据范围外或者位置明显旋转。原因压缩包里的shp经过多次拷贝投影信息文件丢失。解决先根据坐标数值判断。如果shp坐标是像117.5, 26.8这样的小数它大概率是经纬度坐标系可以直接把它转成与DEM一致的投影坐标系。如果坐标是像26551623, 3361584这样的七位数那就是投影坐标。在没有.prj的情况下最稳妥是找相邻区域的相同来源shp读取其.prj内容补到这个shp里。补错了比没有更危险补完一定要重新读取坐标并和DEM叠加验证。4.5 平缓地区坡度图出现“方格子”纹理现象在丘陵缓坡区域坡度图不是平滑的渐变色而是一格一格的锯齿状纹理像像素画。原因DEM可能被前面某步用nearest重采样过或者输入数据是整型存储高程值被取整导致同一个平面区域内只有几个离散的高程值。解决回到原始浮点DEM重新走流程如果没有原始数据可以用gdalwarp -r bilinear再做一次轻量平滑但要意识到这已经改变了数据。还有一个常见来源填洼后没有在原DEM上做坡度而是直接在填洼结果上算坡度那样也会把平地和平缓坡变得过于规整。5. 让30m DEM发挥价值坡度、山体阴影与水文分析预处理干净后DEM才能开始产生实际分析价值。通常我的输出链路是原始DEM用于坡度和山体阴影填洼DEM用于水文分析。两条链分开走原始数据始终不动。5.1 坡度坡向提取先投影再计算计算坡度时工具需要知道像元横向的实际距离。如果数据是经纬度坐标系像元横向单位是度而高程单位是米两者不能直接除。所以第一步要把DEM重投影到米制坐标系再做坡度计算gdaldem slope nanping_dem_projected.tif slope_deg.tif -p -s 1-p表示输出坡度单位是度-s 1表示水平比例因子为1。在投影坐标系中水平单位已经是米因此-s 1是正确写法。如果有人在经纬度数据上直接执行通常会给-s 111320但这个值在南平市的纬度上会有明显误差因为每个纬度的弧长不同。所以我的建议永远是把结论放在更前面先投影再算坡度。5.2 山体阴影一张图暴露90%的DEM问题坡度图只能看出地形陡缓看不出DEM是否出现条带、拼接缝或倒置的负值。山体阴影是更快的地形体检方式gdaldem hillshade nanping_dem_projected.tif hillshade.tif \ -az 315 -alt 45 -z 1.2-az 315表示太阳方位角315度西北方向-alt 45表示太阳高度角45度-z 1.2是垂直夸大系数。处理30m DEM时-z用1到1.5之间比较合适太小看不出地形起伏太大会让人觉得像假山。山体阴影图导出后放大到边界和山脊线附近看如果有明显的直线断层通常就是分幅拼接时的接边误差如果边界一圈出现黑色光晕就是NoData没处理好。5.3 水文分析链路填洼、流向与累积流量这是DEM最高频的应用方向。做完填洼后先计算D8流向栅格再计算累积流量whitebox_tools -r D8Pointer -i nanping_dem_fill.tif \ -o pointer.tif --esri_pntr whitebox_tools -r D8FlowAccumulation -i pointer.tif \ -o flowacc.tif --pntr --esri_pntr第一条命令的输出是每个像元指向其下游像元的方向编码--esri_pntr表示采用ESRI风格的流向编码1表示东2表示东南4表示南8表示西南16表示西32表示西北64表示北128表示东北。第二条命令读取流向栅格累积每个像元上游的像元数量--pntr告诉工具输入的是流向栅格类型。累积流量值越大代表该像元越接近谷线通常设置一个阈值就能提取出河网。这个流程里的坑是如果填洼时--flat_increment设置过大比如0.5米那整片平坦区域都会被抬成微凸的地形流向会沿着假想的“脊线”乱跑提取出的河网会偏离实际位置。我一般先用0.01跑一遍如果累积流量图在平坦区域出现很多条细小平行线再加大增量重跑。6. 最后留一个习惯用shp边界掩膜结果对比确认每一步都没有改坏数据我处理DEM到后面养成了一个习惯每生成一个中间文件就用shp边界把结果栅格化做一次有效像元统计对比。这样做的好处是能快速判断裁剪、重投影和填洼有没有把范围或者数值弄坏。下面这套脚本就是我现在常用的“体检工具”import rasterio import numpy as np import geopandas as gpd from rasterio.features import rasterize def mask_from_shp(shp_path, ref_tif, out_mask): with rasterio.open(ref_tif) as src: meta src.meta.copy() gdf gpd.read_file(shp_path) shapes [g for g in gdf.geometry if g is not None] mask rasterize(shapes, out_shape(src.height, src.width), transformsrc.transform, fill0, default_value1) meta.update(dtypeuint8, count1, nodata0) with rasterio.open(out_mask, w, **meta) as dst: dst.write(mask, 1) def report_dem(path): with rasterio.open(path) as src: dem src.read(1, maskedTrue) valid dem[~dem.mask] print(path) print(有效像元数:, len(valid)) print(有效占比:, (~dem.mask).sum() / dem.size * 100) print(高程范围:, float(np.nanmin(valid)), 到, float(np.nanmax(valid))) print(均值/标准差:, float(np.nanmean(valid)), float(np.nanstd(valid)))mask_from_shp的作用是把shp边界转成和DEM完全同一分辨率的0/1掩膜哪些像元应该在边界内、哪些应该在边界外一目了然。report_dem则算出这个文件里实际有效像元的比例和统计值。如果你裁剪后发现有效占比比掩膜算出的面积少了10%以上那说明裁剪参数或坐标系出了问题。如果高程最小值从0米变成-9999说明NoData又混了进来。数值上的偏差比任何目视检查都更能说明问题。用这套方法我能把每一步的影响控制在可解释范围内。前段时间处理一个类似高程包我就是在填洼后顺手做了一次报告发现最低高程从原来的85米变成了84.5米心里有底了后面水文分析才敢继续用。做这类数据宁可多花几分钟跑一次统计也不要等分析完成后再对着一个异常结果反复猜。希望帮到你。本文还有配套的精品资源点击获取