简介这份广东省30米分辨率DEM高程数据面向GIS从业者、地理科研人员及规划决策者用于地形分析、坡度坡向计算、排水方向模拟与自然灾害风险评估等场景。数据以Shapefile矢量格式组织同时附带GeoTIFF栅格影像可在ArcGIS、QGIS等软件中直接加载支撑地形图制作、基础设施选址与土地利用规划。资源包共26个文件涵盖shp、dbf、shx、prj等Shapefile核心组件以及tif、tfw、ovr、aux等栅格与辅助文件另有adf、xml元数据压缩包约53.97MB结构完整便于直接调用。目前已有1642人学习下载适合需要广东省高精度高程底图的中高级GIS用户快速开展空间分析与专题制图。1. 广东省30米DEM高程数据从拿到数据到能跑地形分析中间隔着什么做华南地区地形分析的人多半遇到过这个场景项目范围划到广东省需要一份能直接进GIS的高程栅格精度要求不算苛刻但也不能太糙30米分辨率刚好卡在“够用”和“体积可控”之间。这份广东省30米DEM高程数据就是冲着这个需求来的——覆盖全省陆域GeoTIFF格式配套shp边界文件拿到手就能往ArcGIS、QGIS或者Python的rasterio里丢。它适合做流域提取、坡度坡向计算、淹没模拟前期底图、风电光伏选址初筛这类活。但别急着双击打开DEM这东西看着简单坐标系、无效值、投影方式任何一个没对齐后面算出来的坡度能离谱到让你怀疑人生。下面按我实际拆包和跑数据的顺序把这份资源从结构到用法到坑一层层说清楚。2. 拆开数据包30米DEM的文件构成与坐标系确认2.1 栅格与矢量文件的配套关系拿到一份区域DEM第一件事不是打开看效果图而是搞清楚里面到底有几个文件、各自什么角色。广东省30米DEM通常以两种形式分发一种是单幅全省拼接的GeoTIFF另一种是按标准分幅切好的多幅瓦片。前者省事后者灵活。配套的shp文件一般是省界或市界矢量用来做裁剪和掩膜。我一般会先用命令行把目录结构过一遍确认文件类型和大小分布# 查看数据包内文件类型与体积分布 find ./guangdong_dem -type f | head -30 du -sh ./guangdong_dem/*如果看到一堆.tif加一个.shp/.dbf/.shx/.prj的组合说明是分幅加边界的标准配置。.prj文件千万别忽略它决定了坐标系定义后面所有分析都依赖这个。2.2 坐标系与投影WGS84还是CGCS2000广东省的DEM数据常见坐标系有两种WGS84地理坐标系EPSG:4326和CGCS2000EPSG:4490部分产品会用Web MercatorEPSG:3857做预处理。30米分辨率下如果原始数据是经纬度存储像元在纬度方向上的实际地面距离会随纬度变化广东省跨北纬20°到25°左右这个差异不能忽略。用gdalinfo快速确认# 查看DEM的坐标系、范围、像元大小和无效值定义 gdalinfo guangdong_dem.tif | grep -E Coordinate|Origin|Pixel|NoData|Band重点看三行Coordinate System确认EPSG代码Pixel Size确认是度还是米NoData Value确认无效值是多少。如果Pixel Size显示的是0.000277...这类小数说明是经纬度存储后续做坡度计算前必须重投影到投影坐标系如UTM 49N/50N或Albers等积投影否则坡度单位是“度/度”没有物理意义。提示广东省东西跨度大UTM分带横跨49N和50N做全省分析建议用Albers等积投影自定义参数避免分带接边处变形。2.3 用Python快速读取并检查数据完整性命令行看完下一步用rasterio做一次程序化体检确认行列数、波段数、统计值和无效值占比import rasterio import numpy as np # 打开DEM文件读取元数据和第一波段 with rasterio.open(guangdong_dem.tif) as src: print(CRS:, src.crs) print(尺寸:, src.width, x, src.height) print(波段数:, src.count) print(像元大小:, src.res) print(无效值:, src.nodata) # 读取第一波段数据 dem src.read(1) # 统计有效值范围 if src.nodata is not None: valid dem[dem ! src.nodata] else: valid dem.flatten() print(最小值:, np.min(valid)) print(最大值:, np.max(valid)) print(均值:, np.mean(valid)) print(无效像元占比: {:.2f}%.format( (dem.size - valid.size) / dem.size * 100 ))这段代码的逻辑很直接先读元数据确认坐标系和无效值定义再把无效值剔除后统计高程范围。广东省陆地高程从沿海接近0米到北部山区约1900米如果统计出来的最大值超过2000或最小值是负数且绝对值很大大概率是无效值没被正确识别或者数据里混入了海洋区域的负值填充。参数说明src.nodata返回的是文件头里定义的无效值常见的是-9999或-32768。如果返回None说明文件没定义无效值需要根据数据分布手动判断。src.res返回像元大小元组地理坐标系下单位是度投影坐标系下单位是米。3. 从DEM到地形因子坡度、坡向与流域提取的实操链路3.1 重投影为什么这一步不能省前面确认了坐标系如果原始数据是EPSG:4326直接算坡度会得到一个无量纲的比值因为水平距离单位是度、高程单位是米量纲不统一。正确做法是先重投影到投影坐标系。用gdalwarp做重投影同时指定重采样方法# 将WGS84地理坐标系重投影为Albers等积投影 # 自定义Albers参数适配广东省范围 gdalwarp -t_srs projaea lat_120 lat_225 lat_022.5 lon_0113 datumWGS84 unitsm \ -tr 30 30 \ -r bilinear \ -srcnodata -9999 -dstnodata -9999 \ -co COMPRESSLZW -co TILEDYES \ guangdong_dem.tif guangdong_dem_aea.tif参数逐个说-t_srs指定目标投影这里用Albers等积投影标准纬线设20°和25°覆盖广东全省中央经线113°大致在省域中部。-tr 30 30强制输出30米像元保持分辨率一致。-r bilinear用双线性插值DEM连续表面适合这种重采样不要用最近邻否则高程会出现阶梯状伪影。-srcnodata和-dstnodata保持无效值一致。-co是创建选项LZW压缩减小体积TILEDYES让后续分块读取更快。重投影完成后用gdalinfo再确认一次Pixel Size应该显示30, 30单位是米。3.2 坡度与坡向计算参数选择决定结果可信度重投影后的DEM就可以算地形因子了。GDAL自带gdaldem工具支持坡度、坡向、山体阴影等# 计算坡度输出单位为度 gdaldem slope guangdong_dem_aea.tif slope_deg.tif -p -s 1.0 -compute_edges # 计算坡向输出0-360度 gdaldem aspect guangdong_dem_aea.tif aspect.tif -compute_edges # 计算山体阴影用于可视化 gdaldem hillshade guangdong_dem_aea.tif hillshade.tif -z 2 -az 315 -alt 45-p表示输出坡度用度为单位不加这个参数默认输出弧度。-s 1.0是高程缩放因子如果DEM高程单位是米、水平单位也是米保持1.0即可。-compute_edges让边缘像元也参与计算否则边界一圈会是无效值。坡向的-compute_edges同理。山体阴影的-z 2是垂直夸张系数2倍能让地形起伏更明显-az 315是光源方位角西北方向-alt 45是光源高度角。算完坡度后建议做一次直方图统计确认分布合理import rasterio import numpy as np import matplotlib.pyplot as plt with rasterio.open(slope_deg.tif) as src: slope src.read(1, maskedTrue) # 绘制坡度分布直方图 plt.hist(slope.compressed(), bins50, edgecolorblack) plt.xlabel(坡度 (度)) plt.ylabel(像元数) plt.title(广东省坡度分布) plt.savefig(slope_hist.png, dpi150)广东省以丘陵和平原为主坡度直方图应该呈现明显的右偏分布——大量像元集中在0到10度尾部延伸到30度以上。如果直方图出现异常尖峰或截断检查是否有无效值混入或重投影时边缘处理不当。3.3 流域提取填洼与流向计算的关键参数用DEM做流域提取是常见需求流程是填洼、算流向、算汇流累积、提取河网。用WhiteboxTools或TauDEM都可以这里用Python的pysheds演示核心步骤import numpy as np from pysheds.grid import Grid # 加载重投影后的DEM grid Grid.from_raster(guangdong_dem_aea.tif) dem grid.read_raster(guangdong_dem_aea.tif) # 填洼处理消除内部凹陷点 flooded grid.fill_depressions(dem) # 计算流向D8算法 fdir grid.flowdir(flooded) # 计算汇流累积量 acc grid.accumulation(fdir) # 按阈值提取河网阈值根据研究区面积调整 # 广东省面积约18万平方公里30米像元下阈值设5000比较合适 streams acc 5000 print(河网像元数:, np.sum(streams))填洼这步不能省原始DEM里总有各种原因造成的凹陷点不填的话流向计算会在这些点中断。fill_depressions用的是优先洪水填充算法计算量不小全省数据建议分块处理或者用WhiteboxTools的FillDepressionsC实现速度快很多。汇流累积阈值决定了河网密度5000这个值对应约4.5平方公里的集水面积是中等分辨率DEM提取河网的常用起点具体项目里要根据目标河流等级调整。注意填洼会改变原始高程如果后续要做淹没模拟应该用原始DEM做水位分析填洼DEM只用于流向和汇流计算。4. 避坑与排查30米DEM处理中最容易翻车的五个地方4.1 无效值被当成真实高程参与计算现象坡度图出现大片异常高值区域或者流域提取结果里出现不合理的长直线。原因DEM文件头里定义了NoData为-9999但某些处理工具默认不识别把-9999当成真实高程参与运算导致局部坡度接近90度。解决每次读取数据后先做无效值掩膜。用rasterio的maskedTrue参数或者手动dem[dem src.nodata] np.nan。GDAL命令行工具用-srcnodata显式指定。处理前用gdalinfo确认NoData值处理后用统计工具验证无效像元占比是否合理。4.2 经纬度数据直接算坡度导致量纲错误现象坡度值全部集中在0到1之间或者出现大量小于0.1的“平地”。原因DEM是EPSG:4326地理坐标系像元大小是度约0.000277度水平距离单位和高程单位不一致GDAL计算时按数值直接算比值得到的是“米/度”这种无意义单位。解决计算任何地形因子之前先重投影到投影坐标系。用gdalwarp转成Albers或UTM确认gdalinfo输出的Pixel Size单位是米。这一步没有捷径所有依赖水平距离的计算都必须先做。4.3 分幅数据拼接后出现接边缝隙现象全省拼接DEM在分幅边界处出现细线状异常坡度图上表现为一条突变的条带。原因相邻分幅的DEM在高程赋值上存在微小差异或者重采样时边缘像元处理方式不一致拼接后接边处不连续。解决优先使用已经拼接好的全省单幅数据。如果只有分幅数据用gdal_merge.py或gdalbuildvrt拼接时加-vrtnodata参数统一无效值拼接后再做一次3x3中值滤波平滑接边处。更彻底的办法是用gdalwarp的-multi和-wo参数做带重叠区的重投影拼接。4.4 填洼算法在大区域上内存溢出现象跑全省填洼时程序崩溃报MemoryError或者进度条卡在某个百分比不动。原因优先洪水填充算法需要维护一个优先队列全省30米DEM约2亿个像元Python实现的内存开销可能超过16GB。解决分块处理用gdalwarp的-te参数按经纬度范围切片每块单独填洼后再拼接。或者换用WhiteboxTools的FillDepressions底层是C实现内存效率高很多。另一个选择是RichDEM支持并行填洼。4.5 汇流累积阈值设错导致河网密度失真现象提取的河网要么密密麻麻像毛细血管要么稀疏到只剩几条主干。原因汇流累积阈值没有根据研究区面积和DEM分辨率换算。不同分辨率的DEM相同阈值对应的集水面积完全不同。解决先算清楚目标集水面积再反推阈值。公式是阈值 目标集水面积 / (像元大小²)。比如想要提取集水面积大于1平方公里的河流30米像元下阈值 1000000 / 900 ≈ 1111。实际项目中一般会试几个阈值叠加影像对比后确定最优值。5. 进阶技巧用DEM做淹没模拟与地形可视化5.1 基于DEM的简易淹没分析拿到重投影后的DEM做淹没模拟其实不复杂。核心逻辑是设定一个水位高程把所有低于该高程的像元标记为淹没区import rasterio import numpy as np def flood_simulation(dem_path, water_level, output_path): 简易淹没模拟标记低于水位的像元 dem_path: DEM文件路径 water_level: 水位高程米 output_path: 输出淹没范围栅格路径 with rasterio.open(dem_path) as src: dem src.read(1, maskedTrue) profile src.profile.copy() # 创建淹没掩膜低于水位且非无效值 flooded (dem water_level) (~dem.mask) # 输出为8位栅格1表示淹没区0表示非淹没区 profile.update(dtyperasterio.uint8, nodata255, compresslzw) with rasterio.open(output_path, w, **profile) as dst: dst.write(flooded.astype(rasterio.uint8), 1) # 统计淹没面积 with rasterio.open(output_path) as src: arr src.read(1) pixel_area src.res[0] * src.res[1] # 像元面积平方米 flooded_pixels np.sum(arr 1) print(f水位 {water_level} 米时淹没面积: {flooded_pixels * pixel_area / 1e6:.2f} 平方公里) # 模拟水位10米时的淹没范围 flood_simulation(guangdong_dem_aea.tif, 10, flood_10m.tif)这段代码的关键在于dem.mask的使用——rasterio读取时如果指定maskedTrue会自动把NoData像元掩膜掉~dem.mask确保无效值不会被误判为淹没区。输出用uint8类型体积小方便后续做叠加分析。像元面积从src.res获取重投影后单位是米乘积就是平方米。这个简易模型没有考虑水文连通性——一个低洼内陆盆地即使低于水位也不会被淹没因为它和海洋不连通。如果需要考虑连通性要用scipy.ndimage.label做连通域分析只保留与边界连通的淹没区。5.2 地形可视化山体阴影叠加高程着色单纯的高程灰度图看不出地形细节把山体阴影和高程着色叠加是标准做法import rasterio import numpy as np import matplotlib.pyplot as plt from matplotlib.colors import LightSource # 读取DEM和山体阴影 with rasterio.open(guangdong_dem_aea.tif) as src: dem src.read(1, maskedTrue) extent [src.bounds.left, src.bounds.right, src.bounds.bottom, src.bounds.top] with rasterio.open(hillshade.tif) as src: hillshade src.read(1, maskedTrue) # 创建地形着色 ls LightSource(azdeg315, altdeg45) rgb ls.shade(dem, cmapplt.cm.terrain, vert_exag2, blend_modesoft) # 叠加山体阴影增强立体感 fig, ax plt.subplots(figsize(16, 12)) ax.imshow(rgb, extentextent) ax.imshow(hillshade, extentextent, cmapgray, alpha0.3) ax.set_title(广东省30米DEM地形可视化) ax.set_xlabel(东向坐标 (米)) ax.set_ylabel(北向坐标 (米)) plt.tight_layout() plt.savefig(guangdong_terrain.png, dpi200)LightSource.shade的vert_exag2做2倍垂直夸张blend_modesoft让着色和阴影融合更自然。叠加时山体阴影用alpha0.3降低透明度避免盖住高程颜色信息。这样出来的图能同时看到高程分布和地形起伏细节。5.3 一个我踩过的坑早些年做某个区域的地形分析DEM拿到手直接算坡度结果出来一片“平地”最大值不到2度。排查了半天以为是数据质量问题后来才发现原始数据是EPSG:4326像元大小0.000277度GDAL按数值算比值坡度自然小得离谱。重投影到UTM后同一区域坡度最大值跳到35度完全对得上实地情况。从那以后我每次拿到DEM第一件事就是gdalinfo看坐标系和像元大小确认单位是米之后才往下走。这个习惯帮我省了至少三次返工。希望这份广东省30米DEM数据和上面的处理链路能让你少走一段弯路。本文还有配套的精品资源点击获取 
长河君 · 2026/10/13 15:24:28 · 阅读 1.2k 广东省30米DEM高程数据：从数据包到地形分析全流程
简介这份广东省30米分辨率DEM高程数据面向GIS从业者、地理科研人员及规划决策者用于地形分析、坡度坡向计算、排水方向模拟与自然灾害风险评估等场景。数据以Shapefile矢量格式组织同时附带GeoTIFF栅格影像可在ArcGIS、QGIS等软件中直接加载支撑地形图制作、基础设施选址与土地利用规划。资源包共26个文件涵盖shp、dbf、shx、prj等Shapefile核心组件以及tif、tfw、ovr、aux等栅格与辅助文件另有adf、xml元数据压缩包约53.97MB结构完整便于直接调用。目前已有1642人学习下载适合需要广东省高精度高程底图的中高级GIS用户快速开展空间分析与专题制图。1. 广东省30米DEM高程数据从拿到数据到能跑地形分析中间隔着什么做华南地区地形分析的人多半遇到过这个场景项目范围划到广东省需要一份能直接进GIS的高程栅格精度要求不算苛刻但也不能太糙30米分辨率刚好卡在“够用”和“体积可控”之间。这份广东省30米DEM高程数据就是冲着这个需求来的——覆盖全省陆域GeoTIFF格式配套shp边界文件拿到手就能往ArcGIS、QGIS或者Python的rasterio里丢。它适合做流域提取、坡度坡向计算、淹没模拟前期底图、风电光伏选址初筛这类活。但别急着双击打开DEM这东西看着简单坐标系、无效值、投影方式任何一个没对齐后面算出来的坡度能离谱到让你怀疑人生。下面按我实际拆包和跑数据的顺序把这份资源从结构到用法到坑一层层说清楚。2. 拆开数据包30米DEM的文件构成与坐标系确认2.1 栅格与矢量文件的配套关系拿到一份区域DEM第一件事不是打开看效果图而是搞清楚里面到底有几个文件、各自什么角色。广东省30米DEM通常以两种形式分发一种是单幅全省拼接的GeoTIFF另一种是按标准分幅切好的多幅瓦片。前者省事后者灵活。配套的shp文件一般是省界或市界矢量用来做裁剪和掩膜。我一般会先用命令行把目录结构过一遍确认文件类型和大小分布# 查看数据包内文件类型与体积分布 find ./guangdong_dem -type f | head -30 du -sh ./guangdong_dem/*如果看到一堆.tif加一个.shp/.dbf/.shx/.prj的组合说明是分幅加边界的标准配置。.prj文件千万别忽略它决定了坐标系定义后面所有分析都依赖这个。2.2 坐标系与投影WGS84还是CGCS2000广东省的DEM数据常见坐标系有两种WGS84地理坐标系EPSG:4326和CGCS2000EPSG:4490部分产品会用Web MercatorEPSG:3857做预处理。30米分辨率下如果原始数据是经纬度存储像元在纬度方向上的实际地面距离会随纬度变化广东省跨北纬20°到25°左右这个差异不能忽略。用gdalinfo快速确认# 查看DEM的坐标系、范围、像元大小和无效值定义 gdalinfo guangdong_dem.tif | grep -E Coordinate|Origin|Pixel|NoData|Band重点看三行Coordinate System确认EPSG代码Pixel Size确认是度还是米NoData Value确认无效值是多少。如果Pixel Size显示的是0.000277...这类小数说明是经纬度存储后续做坡度计算前必须重投影到投影坐标系如UTM 49N/50N或Albers等积投影否则坡度单位是“度/度”没有物理意义。提示广东省东西跨度大UTM分带横跨49N和50N做全省分析建议用Albers等积投影自定义参数避免分带接边处变形。2.3 用Python快速读取并检查数据完整性命令行看完下一步用rasterio做一次程序化体检确认行列数、波段数、统计值和无效值占比import rasterio import numpy as np # 打开DEM文件读取元数据和第一波段 with rasterio.open(guangdong_dem.tif) as src: print(CRS:, src.crs) print(尺寸:, src.width, x, src.height) print(波段数:, src.count) print(像元大小:, src.res) print(无效值:, src.nodata) # 读取第一波段数据 dem src.read(1) # 统计有效值范围 if src.nodata is not None: valid dem[dem ! src.nodata] else: valid dem.flatten() print(最小值:, np.min(valid)) print(最大值:, np.max(valid)) print(均值:, np.mean(valid)) print(无效像元占比: {:.2f}%.format( (dem.size - valid.size) / dem.size * 100 ))这段代码的逻辑很直接先读元数据确认坐标系和无效值定义再把无效值剔除后统计高程范围。广东省陆地高程从沿海接近0米到北部山区约1900米如果统计出来的最大值超过2000或最小值是负数且绝对值很大大概率是无效值没被正确识别或者数据里混入了海洋区域的负值填充。参数说明src.nodata返回的是文件头里定义的无效值常见的是-9999或-32768。如果返回None说明文件没定义无效值需要根据数据分布手动判断。src.res返回像元大小元组地理坐标系下单位是度投影坐标系下单位是米。3. 从DEM到地形因子坡度、坡向与流域提取的实操链路3.1 重投影为什么这一步不能省前面确认了坐标系如果原始数据是EPSG:4326直接算坡度会得到一个无量纲的比值因为水平距离单位是度、高程单位是米量纲不统一。正确做法是先重投影到投影坐标系。用gdalwarp做重投影同时指定重采样方法# 将WGS84地理坐标系重投影为Albers等积投影 # 自定义Albers参数适配广东省范围 gdalwarp -t_srs projaea lat_120 lat_225 lat_022.5 lon_0113 datumWGS84 unitsm \ -tr 30 30 \ -r bilinear \ -srcnodata -9999 -dstnodata -9999 \ -co COMPRESSLZW -co TILEDYES \ guangdong_dem.tif guangdong_dem_aea.tif参数逐个说-t_srs指定目标投影这里用Albers等积投影标准纬线设20°和25°覆盖广东全省中央经线113°大致在省域中部。-tr 30 30强制输出30米像元保持分辨率一致。-r bilinear用双线性插值DEM连续表面适合这种重采样不要用最近邻否则高程会出现阶梯状伪影。-srcnodata和-dstnodata保持无效值一致。-co是创建选项LZW压缩减小体积TILEDYES让后续分块读取更快。重投影完成后用gdalinfo再确认一次Pixel Size应该显示30, 30单位是米。3.2 坡度与坡向计算参数选择决定结果可信度重投影后的DEM就可以算地形因子了。GDAL自带gdaldem工具支持坡度、坡向、山体阴影等# 计算坡度输出单位为度 gdaldem slope guangdong_dem_aea.tif slope_deg.tif -p -s 1.0 -compute_edges # 计算坡向输出0-360度 gdaldem aspect guangdong_dem_aea.tif aspect.tif -compute_edges # 计算山体阴影用于可视化 gdaldem hillshade guangdong_dem_aea.tif hillshade.tif -z 2 -az 315 -alt 45-p表示输出坡度用度为单位不加这个参数默认输出弧度。-s 1.0是高程缩放因子如果DEM高程单位是米、水平单位也是米保持1.0即可。-compute_edges让边缘像元也参与计算否则边界一圈会是无效值。坡向的-compute_edges同理。山体阴影的-z 2是垂直夸张系数2倍能让地形起伏更明显-az 315是光源方位角西北方向-alt 45是光源高度角。算完坡度后建议做一次直方图统计确认分布合理import rasterio import numpy as np import matplotlib.pyplot as plt with rasterio.open(slope_deg.tif) as src: slope src.read(1, maskedTrue) # 绘制坡度分布直方图 plt.hist(slope.compressed(), bins50, edgecolorblack) plt.xlabel(坡度 (度)) plt.ylabel(像元数) plt.title(广东省坡度分布) plt.savefig(slope_hist.png, dpi150)广东省以丘陵和平原为主坡度直方图应该呈现明显的右偏分布——大量像元集中在0到10度尾部延伸到30度以上。如果直方图出现异常尖峰或截断检查是否有无效值混入或重投影时边缘处理不当。3.3 流域提取填洼与流向计算的关键参数用DEM做流域提取是常见需求流程是填洼、算流向、算汇流累积、提取河网。用WhiteboxTools或TauDEM都可以这里用Python的pysheds演示核心步骤import numpy as np from pysheds.grid import Grid # 加载重投影后的DEM grid Grid.from_raster(guangdong_dem_aea.tif) dem grid.read_raster(guangdong_dem_aea.tif) # 填洼处理消除内部凹陷点 flooded grid.fill_depressions(dem) # 计算流向D8算法 fdir grid.flowdir(flooded) # 计算汇流累积量 acc grid.accumulation(fdir) # 按阈值提取河网阈值根据研究区面积调整 # 广东省面积约18万平方公里30米像元下阈值设5000比较合适 streams acc 5000 print(河网像元数:, np.sum(streams))填洼这步不能省原始DEM里总有各种原因造成的凹陷点不填的话流向计算会在这些点中断。fill_depressions用的是优先洪水填充算法计算量不小全省数据建议分块处理或者用WhiteboxTools的FillDepressionsC实现速度快很多。汇流累积阈值决定了河网密度5000这个值对应约4.5平方公里的集水面积是中等分辨率DEM提取河网的常用起点具体项目里要根据目标河流等级调整。注意填洼会改变原始高程如果后续要做淹没模拟应该用原始DEM做水位分析填洼DEM只用于流向和汇流计算。4. 避坑与排查30米DEM处理中最容易翻车的五个地方4.1 无效值被当成真实高程参与计算现象坡度图出现大片异常高值区域或者流域提取结果里出现不合理的长直线。原因DEM文件头里定义了NoData为-9999但某些处理工具默认不识别把-9999当成真实高程参与运算导致局部坡度接近90度。解决每次读取数据后先做无效值掩膜。用rasterio的maskedTrue参数或者手动dem[dem src.nodata] np.nan。GDAL命令行工具用-srcnodata显式指定。处理前用gdalinfo确认NoData值处理后用统计工具验证无效像元占比是否合理。4.2 经纬度数据直接算坡度导致量纲错误现象坡度值全部集中在0到1之间或者出现大量小于0.1的“平地”。原因DEM是EPSG:4326地理坐标系像元大小是度约0.000277度水平距离单位和高程单位不一致GDAL计算时按数值直接算比值得到的是“米/度”这种无意义单位。解决计算任何地形因子之前先重投影到投影坐标系。用gdalwarp转成Albers或UTM确认gdalinfo输出的Pixel Size单位是米。这一步没有捷径所有依赖水平距离的计算都必须先做。4.3 分幅数据拼接后出现接边缝隙现象全省拼接DEM在分幅边界处出现细线状异常坡度图上表现为一条突变的条带。原因相邻分幅的DEM在高程赋值上存在微小差异或者重采样时边缘像元处理方式不一致拼接后接边处不连续。解决优先使用已经拼接好的全省单幅数据。如果只有分幅数据用gdal_merge.py或gdalbuildvrt拼接时加-vrtnodata参数统一无效值拼接后再做一次3x3中值滤波平滑接边处。更彻底的办法是用gdalwarp的-multi和-wo参数做带重叠区的重投影拼接。4.4 填洼算法在大区域上内存溢出现象跑全省填洼时程序崩溃报MemoryError或者进度条卡在某个百分比不动。原因优先洪水填充算法需要维护一个优先队列全省30米DEM约2亿个像元Python实现的内存开销可能超过16GB。解决分块处理用gdalwarp的-te参数按经纬度范围切片每块单独填洼后再拼接。或者换用WhiteboxTools的FillDepressions底层是C实现内存效率高很多。另一个选择是RichDEM支持并行填洼。4.5 汇流累积阈值设错导致河网密度失真现象提取的河网要么密密麻麻像毛细血管要么稀疏到只剩几条主干。原因汇流累积阈值没有根据研究区面积和DEM分辨率换算。不同分辨率的DEM相同阈值对应的集水面积完全不同。解决先算清楚目标集水面积再反推阈值。公式是阈值 目标集水面积 / (像元大小²)。比如想要提取集水面积大于1平方公里的河流30米像元下阈值 1000000 / 900 ≈ 1111。实际项目中一般会试几个阈值叠加影像对比后确定最优值。5. 进阶技巧用DEM做淹没模拟与地形可视化5.1 基于DEM的简易淹没分析拿到重投影后的DEM做淹没模拟其实不复杂。核心逻辑是设定一个水位高程把所有低于该高程的像元标记为淹没区import rasterio import numpy as np def flood_simulation(dem_path, water_level, output_path): 简易淹没模拟标记低于水位的像元 dem_path: DEM文件路径 water_level: 水位高程米 output_path: 输出淹没范围栅格路径 with rasterio.open(dem_path) as src: dem src.read(1, maskedTrue) profile src.profile.copy() # 创建淹没掩膜低于水位且非无效值 flooded (dem water_level) (~dem.mask) # 输出为8位栅格1表示淹没区0表示非淹没区 profile.update(dtyperasterio.uint8, nodata255, compresslzw) with rasterio.open(output_path, w, **profile) as dst: dst.write(flooded.astype(rasterio.uint8), 1) # 统计淹没面积 with rasterio.open(output_path) as src: arr src.read(1) pixel_area src.res[0] * src.res[1] # 像元面积平方米 flooded_pixels np.sum(arr 1) print(f水位 {water_level} 米时淹没面积: {flooded_pixels * pixel_area / 1e6:.2f} 平方公里) # 模拟水位10米时的淹没范围 flood_simulation(guangdong_dem_aea.tif, 10, flood_10m.tif)这段代码的关键在于dem.mask的使用——rasterio读取时如果指定maskedTrue会自动把NoData像元掩膜掉~dem.mask确保无效值不会被误判为淹没区。输出用uint8类型体积小方便后续做叠加分析。像元面积从src.res获取重投影后单位是米乘积就是平方米。这个简易模型没有考虑水文连通性——一个低洼内陆盆地即使低于水位也不会被淹没因为它和海洋不连通。如果需要考虑连通性要用scipy.ndimage.label做连通域分析只保留与边界连通的淹没区。5.2 地形可视化山体阴影叠加高程着色单纯的高程灰度图看不出地形细节把山体阴影和高程着色叠加是标准做法import rasterio import numpy as np import matplotlib.pyplot as plt from matplotlib.colors import LightSource # 读取DEM和山体阴影 with rasterio.open(guangdong_dem_aea.tif) as src: dem src.read(1, maskedTrue) extent [src.bounds.left, src.bounds.right, src.bounds.bottom, src.bounds.top] with rasterio.open(hillshade.tif) as src: hillshade src.read(1, maskedTrue) # 创建地形着色 ls LightSource(azdeg315, altdeg45) rgb ls.shade(dem, cmapplt.cm.terrain, vert_exag2, blend_modesoft) # 叠加山体阴影增强立体感 fig, ax plt.subplots(figsize(16, 12)) ax.imshow(rgb, extentextent) ax.imshow(hillshade, extentextent, cmapgray, alpha0.3) ax.set_title(广东省30米DEM地形可视化) ax.set_xlabel(东向坐标 (米)) ax.set_ylabel(北向坐标 (米)) plt.tight_layout() plt.savefig(guangdong_terrain.png, dpi200)LightSource.shade的vert_exag2做2倍垂直夸张blend_modesoft让着色和阴影融合更自然。叠加时山体阴影用alpha0.3降低透明度避免盖住高程颜色信息。这样出来的图能同时看到高程分布和地形起伏细节。5.3 一个我踩过的坑早些年做某个区域的地形分析DEM拿到手直接算坡度结果出来一片“平地”最大值不到2度。排查了半天以为是数据质量问题后来才发现原始数据是EPSG:4326像元大小0.000277度GDAL按数值算比值坡度自然小得离谱。重投影到UTM后同一区域坡度最大值跳到35度完全对得上实地情况。从那以后我每次拿到DEM第一件事就是gdalinfo看坐标系和像元大小确认单位是米之后才往下走。这个习惯帮我省了至少三次返工。希望这份广东省30米DEM数据和上面的处理链路能让你少走一段弯路。本文还有配套的精品资源点击获取