ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

梧州DEM数据处理全流程:从坐标系检查到坡度河网提取与批量出图

梧州DEM数据处理全流程:从坐标系检查到坡度河网提取与批量出图 简介这份资源面向地理信息、城乡规划、环境与测绘等方向的学习者和从业者提供广西梧州市30米分辨率的DEM数字高程数据并附带市域边界矢量文件可用于高程提取、坡度坡向计算、排水网络模拟及灾害风险评估等分析场景。压缩包共12个文件约49.43MB以tif栅格为主辅以prj投影、tfw地理坐标、ovr快速预览、xml元数据以及shp、shx、dbf组成的边界矢量文件另含sbn、sbx索引文件方便在ArcGIS、QGIS等软件中直接加载与叠加分析。目前已有367人学习下载。数据覆盖梧州市行政范围并外扩至周边区域便于获得更完整的地理背景适合作为GIS课程练习、项目底图或区域研究的基础地形资料帮助读者快速上手空间数据处理与地形分析流程。1. 梧州 DEM 数据到手后先搞清楚这套 tif shp 能干什么做广西地形分析的人多半都经历过满网找 DEM 的绝望。这次拆的是一份梧州市的 DEM 数字高程 tif 数据附带市范围的 shp 文件。拿到手第一件事不是急着打开 ArcMap而是先判断它能不能用。tif 是栅格高程每个像元存一个海拔值shp 是矢量边界圈定梧州市的行政范围。两者配合能干的活很具体坡度坡向提取、流域划分、淹没模拟、选址分析、三维地形底图。适合做 GIS 教学、区域规划、水文预研的从业者。但别指望它直接给你 DSM 那种带建筑的高度——DEM 是裸地形这是两码事。下面按「是什么、怎么用、坑在哪」一路拆下去。2. 梧州 DEM 与 shp 的坐标系、分辨率与选型逻辑2.1 为什么先看坐标系而不是先看分辨率很多人拿到 tif 第一反应是看分辨率其实坐标系错了分辨率再高也是废的。梧州地处东经 110°18′ 到 111°40′北纬 22°37′ 到 24°10′属于典型的华南丘陵地带。国内常见的 DEM 坐标系有两种地理坐标系 WGS84经纬度单位度和投影坐标系如 CGCS2000 3 度带或 UTM 49N单位米。如果你拿到的 tif 是经纬度坐标直接算坡度会出大问题——因为经纬度不是等距的纬度方向一度约 111 公里经度方向在梧州纬度约 102 公里ArcMap 算坡度时按像元尺寸当平面距离算结果会偏。判断方法很简单在 ArcMap 里右键图层看属性或者用 gdalinfo 命令行gdalinfo wuzhou_dem.tif输出里看 Coordinate System 那一行。如果是GEOGCS[WGS 84]单位是 degree如果是PROJCS[CGCS2000_3_Degree_GK_Zone_37]之类单位是 metre。梧州落在 CGCS2000 3 度带第 37 带中央经线 111°E常见做法是先把经纬度 DEM 投影到 CGCS2000 3 度带再做地形分析。投影用 ArcToolbox 的 Project Raster或者 gdalwarpgdalwarp -t_srs EPSG:4547 -r bilinear wuzhou_dem.tif wuzhou_dem_proj.tifEPSG:4547 是 CGCS2000 3 度带 CM 111E 的编码。-r bilinear是重采样方式高程连续数据用双线性比最近邻更平滑。这一步做完像元尺寸才是真正的米坡度单位才是度或百分比。2.2 分辨率决定了你能做多大尺度的分析这份梧州 DEM 的分辨率常见的有 30 米SRTM、ASTER和 12.5 米ALOS两种。30 米分辨率意味着一个像元代表地面 30×30 米梧州总面积约 1.26 万平方公里算下来大概 1400 万个像元文件大小在几十到一百多 MB 之间。12.5 米的话像元数翻五倍多文件更大。分辨率直接决定分析尺度。做全市域坡度分级、流域划分30 米够用做小流域洪水淹没、单个园区选址30 米就太粗了河网会断、坡长失真。我一般会先算一下研究区面积和最小地形单元尺寸如果最小关心的地形特征小于 5 个像元就该换更高分辨率的数据。这份数据如果只有 30 米别硬拿去做村级规划。2.3 shp 边界文件的两个用途和常见误用市范围 shp 不只是用来出图好看的。它有两个硬用途一是裁剪把 DEM 裁到梧州行政边界内减少数据量、避免邻区干扰二是分区统计按县区统计平均高程、坡度分布。但这里有个高频翻车点shp 和 tif 坐标系不一致。如果 shp 是 WGS84 经纬度tif 被你投影成了 CGCS2000 米制直接裁剪会报错或裁出空白。正确做法是先统一坐标系再裁剪。另外shp 的边界如果是行政区划线裁剪时边缘像元会被切掉一半做统计时边缘县区的均值会有偏差这个后面避坑章节细说。3. 用 ArcMap 裁剪 DEM掩码提取和按面裁剪到底选哪个3.1 两种裁剪方式的本质区别热搜里有人问「arcmap 中依靠面图层裁剪 dem 栅格 tif 文件和依靠面图层掩码提取有啥区别」这个问题问到点子上了。ArcMap 里有两个工具长得像但行为不同Clip数据管理 → 栅格 → 栅格处理 → 裁剪按面要素的几何边界切割栅格输出范围是面的最小外接矩形面外的像元被设为 NoData。注意它输出的是矩形范围不是精确贴合边界。Extract by Mask空间分析 → 提取分析 → 按掩码提取同样把面外像元设为 NoData但它是逐像元判断输出范围也是矩形但 NoData 区域严格按面边界。实际差别在哪Clip 在处理复杂边界时边缘像元可能保留不完整Extract by Mask 更精确但需要 Spatial Analyst 许可。我一般做正式分析用 Extract by Mask做快速预览用 Clip。# ArcPy 方式按掩码提取适合批量处理 import arcpy from arcpy.sa import ExtractByMask arcpy.env.workspace rD:\wuzhou arcpy.CheckOutExtension(Spatial) dem rwuzhou_dem_proj.tif mask rwuzhou_boundary.shp out ExtractByMask(dem, mask) out.save(rwuzhou_dem_clipped.tif)这段代码先签出 Spatial 扩展许可然后调用 ExtractByMask。dem是投影后的 DEMmask是市范围 shp。输出保存为裁剪后的 tif。注意arcpy.env.workspace设成你的实际路径路径里别带中文ArcPy 对中文路径的支持时好时坏这是血泪经验。3.2 裁剪后的 NoData 处理和统计陷阱裁剪完你会发现边界外是 NoData边界内是正常高程。这时候做统计要小心ArcMap 的栅格统计默认忽略 NoData但如果你用 Python 的 rasterio 或 numpy 读NoData 可能被当成 0 或 -9999直接算均值就废了。import rasterio import numpy as np with rasterio.open(wuzhou_dem_clipped.tif) as src: data src.read(1) nodata src.nodata print(NoData value:, nodata) valid data[data ! nodata] print(有效像元数:, valid.size) print(平均高程:, valid.mean()) print(最高点:, valid.max(), 最低点:, valid.min())src.nodata读出 NoData 值然后用布尔索引过滤。如果 nodata 是 None说明文件没定义 NoData这时候要手动判断——常见做法是看数据里有没有 -9999 或 -32768 这种异常值。梧州最低海拔在几十米最高在 1200 米左右具体看数据超出这个范围的极值基本就是 NoData 没处理好。3.3 按县区批量统计高程与坡度梧州市下辖万秀、长洲、龙圩三个城区和苍梧、藤县、蒙山、岑溪等县市。如果你有县区级 shp可以用 Zonal Statistics as Table 批量统计。import arcpy from arcpy.sa import ZonalStatisticsAsTable arcpy.CheckOutExtension(Spatial) dem rwuzhou_dem_clipped.tif zones rwuzhou_counties.shp out_table rwuzhou_stats.dbf ZonalStatisticsAsTable(zones, NAME, dem, out_table, DATA, ALL)NAME是 shp 里的县区名字段ALL表示统计所有指标均值、最大、最小、标准差等。输出是 dbf 表可以导进 Excel 看。这一步做完你就能知道哪个县平均海拔最高、哪个县地形起伏最大。蒙山县在梧州北部山地多平均高程应该明显高于市区。4. 从 DEM 提取坡度、坡向与河网参数怎么设才不翻车4.1 坡度计算Z 因子和单位选择坡度是 DEM 最常用的派生数据。ArcMap 的 Slope 工具里有个 Z factor 参数很多人直接默认 1结果算出来的坡度偏小。Z factor 是高程单位与平面单位的比值。如果你的 DEM 是米制投影坐标高程也是米Z factor 1 没问题。但如果 DEM 是经纬度坐标平面单位是度高程是米就必须设 Z factor否则坡度完全不对。梧州纬度约 23°经纬度 DEM 的 Z factor 大约等于 1/111320每度约 111.32 公里但更稳妥的做法是先投影再算坡度别在经纬度上硬算。from arcpy.sa import Slope from arcpy import CheckOutExtension CheckOutExtension(Spatial) dem_proj rwuzhou_dem_proj.tif slope_deg Slope(dem_proj, DEGREE, 1) slope_deg.save(rwuzhou_slope_deg.tif)DEGREE输出角度制1是 Z factor。如果输出百分比坡度改成PERCENT_RISE。梧州丘陵区坡度多在 5° 到 25° 之间超过 25° 的陡坡集中在山区这些区域做建设选址要避开。4.2 坡向提取与阴坡阳坡分析坡向用 Aspect 工具输出 0 到 360 度0 是正北90 是正东。梧州属亚热带阳坡南向、东南向蒸发强、土壤干阴坡北向、西北向湿度大、植被好。做农业或生态分析时坡向常和坡度叠加使用。from arcpy.sa import Aspect aspect Aspect(rwuzhou_dem_proj.tif) aspect.save(rwuzhou_aspect.tif)坡向数据有个坑平坦区域坡度接近 0的坡向是随机的因为算法靠邻域像元高差判断方向平地没有明显方向。所以坡向分析前最好先用坡度做个掩码把坡度小于 1° 的区域排除。4.3 河网提取填洼、流向、流量累积三步走从 DEM 提取河网是水文分析的经典流程也是翻车重灾区。标准步骤是填洼Fill→ 流向Flow Direction→ 流量累积Flow Accumulation→ 栅格河网Raster Calculator 阈值→ 矢量化Stream to Feature。from arcpy.sa import Fill, FlowDirection, FlowAccumulation, Con from arcpy import CheckOutExtension CheckOutExtension(Spatial) dem rwuzhou_dem_proj.tif filled Fill(dem) filled.save(rwuzhou_filled.tif) flow_dir FlowDirection(filled, NORMAL) flow_dir.save(rwuzhou_flowdir.tif) flow_acc FlowAccumulation(flow_dir) flow_acc.save(rwuzhou_flowacc.tif) # 阈值设 500意味着集水面积超过 500 个像元的才算河网 stream Con(flow_acc 500, 1) stream.save(rwuzhou_stream.tif)Fill填平洼地避免水流断头。FlowDirection用NORMAL表示按最大坡度下降。FlowAccumulation算每个像元上游汇水像元数。阈值 500 是经验值30 米分辨率下大约对应 0.45 平方公里集水面积。阈值越小河网越密但假河道也越多。梧州这种湿润区阈值可以设 300 到 800 之间具体看你要多细的河网。5. 避坑与排查梧州 DEM 处理中最容易翻车的五件事5.1 裁剪后统计均值偏高或偏低现象裁剪后的 DEM 算平均高程和裁剪前差了好几十米。原因shp 边界外的 NoData 被当成了 0 参与计算或者 shp 本身有重叠面、空洞。解决先检查 shp 的几何有效性用 Repair Geometry 修一遍统计时确认 NoData 被正确识别Python 里用data ! nodata过滤ArcMap 里确认环境设置中 NoData 处理方式。5.2 坡度图出现网格状条纹现象坡度图上有规律的方格纹理像马赛克。原因DEM 本身有噪声或者投影时用了最近邻重采样导致高程值出现阶梯。解决投影重采样用双线性或三次卷积如果 DEM 噪声大先做一次低通滤波Focal Statistics窗口 3×3均值再算坡度。但滤波会平滑真实地形慎用。5.3 河网提取出来是断线或乱线现象提取的河网不连续或者平地上出现大量平行线。原因填洼不彻底或者流向算法在平坦区失效。解决填洼时检查是否有大范围平地必要时用Fill后再做一次FlowDirection验证平坦区可以用FORCE选项强制流向但会引入人为偏差。更稳妥的做法是换用 TauDEM 或 WhiteboxTools 的河网提取算法对平地处理更好。5.4 shp 和 tif 坐标系不一致导致裁剪空白现象裁剪结果全是 NoData或者只裁出一小块。原因shp 是地理坐标系tif 是投影坐标系ArcMap 虽然能显示在一起但裁剪时按各自坐标算对不上。解决用 Project 工具把 shp 投影到和 tif 一样的坐标系或者反过来。别依赖 ArcMap 的实时投影裁剪工具不认。5.5 大文件处理时 ArcMap 崩溃现象处理 12.5 米分辨率 DEM 时 ArcMap 无响应或报内存错误。原因32 位 ArcMap 内存上限约 2GB大栅格运算容易爆。解决换 ArcGIS Pro64 位或者用 Python rasterio/GDAL 在命令行处理。如果必须用 ArcMap先把 DEM 裁到研究区范围再算别拿全市数据直接跑流量累积。6. 进阶技巧用 Python 批量出图与高程分级配色前面几步走完你手里应该有裁剪后的 DEM、坡度、坡向、河网。最后一章说个实际工作中提效的技巧批量出图。做区域规划时经常要按县区各出一张地形图手动在 ArcMap 里调符号、加图例、导出一个县五分钟七个县就是半小时。用 Python 脚本可以压到两分钟。核心思路是用 arcpy.mp 操作 ArcGIS Pro 的工程文件或者用 matplotlib rasterio 直接画。后者不依赖 ArcGIS 许可更适合批量。import rasterio import matplotlib.pyplot as plt import numpy as np from matplotlib.colors import LinearSegmentedColormap # 自定义高程配色低海拔绿中海拔黄高海拔棕 colors [#2d7f3e, #8fbc5a, #e8d96a, #c9a24b, #8b5a2b, #ffffff] cmap LinearSegmentedColormap.from_list(elevation, colors) with rasterio.open(wuzhou_dem_clipped.tif) as src: dem src.read(1) nodata src.nodata dem_masked np.ma.masked_equal(dem, nodata) fig, ax plt.subplots(figsize(10, 8), dpi150) im ax.imshow(dem_masked, cmapcmap) ax.set_title(Wuzhou DEM Elevation, fontsize14) ax.axis(off) cbar plt.colorbar(im, axax, shrink0.7) cbar.set_label(Elevation (m)) plt.tight_layout() plt.savefig(wuzhou_dem_map.png, bbox_inchestight) plt.close()这段代码用 rasterio 读 tifnumpy 做 NoData 掩码matplotlib 画图。LinearSegmentedColormap.from_list定义了一个从绿到白的渐变色带比默认的 jet 配色更符合地形直觉。bbox_inchestight去掉多余白边。输出是 PNG可以直接插进报告。如果要按县区批量出图把 shp 读进来循环每个县做掩码裁剪再画。geopandas 读 shprasterio.mask 做裁剪流程和上面一样只是多一层循环。import geopandas as gpd from rasterio.mask import mask counties gpd.read_file(wuzhou_counties.shp) with rasterio.open(wuzhou_dem_clipped.tif) as src: for idx, row in counties.iterrows(): geom [row.geometry] out_image, out_transform mask(src, geom, cropTrue) out_meta src.meta.copy() out_meta.update({ height: out_image.shape[1], width: out_image.shape[2], transform: out_transform }) name row[NAME] with rasterio.open(fdem_{name}.tif, w) as dst: dst.write(out_image)mask函数按几何裁剪cropTrue裁掉外围。out_meta更新尺寸和变换矩阵。每个县输出一个 tif后续可以单独出图或统计。注意row[NAME]要换成你 shp 里实际的字段名别照抄。从那以后我每次拿到新 DEM都强制走一遍「gdalinfo 看坐标系 → 投影到米制 → 裁剪 → 坡度验证」的流程不再凭感觉直接开算。希望帮到你。本文还有配套的精品资源点击获取
返回列表