简介这份广东省30米分辨率DEM高程数据面向GIS从业者、地理科研人员及规划决策者用于地形分析、坡度坡向计算、排水方向模拟与自然灾害风险评估等场景。数据以Shapefile矢量格式组织同时附带GeoTIFF栅格文件可在ArcGIS、QGIS等软件中直接导入与人口、交通、气候等数据叠加开展空间分析。资源包共26个文件涵盖shp、dbf、shx、prj等Shapefile核心组件以及tif、adf、ovr、xml、tfw等栅格与元数据文件压缩包约53.97MB目录结构完整便于按需检索。目前已有1642人学习下载适合需要广东省高精度地形底图的中高级GIS用户可快速支撑地形图制作、基础设施选址与土地利用规划等工作。1. 广东省30米高清DEM高程数据从下载到能用的地形分析全流程做华南地区任何一个跟地形沾边的项目绕不开的第一份底图就是广东省30米高清DEM高程数据。我最早接触它是因为一个洪水淹没模拟的活儿——拿到手才发现网上讲DEM的文章要么只讲概念要么直接甩个下载链接就完事真正从数据源选型、坐标系统一到坡度坡向计算、再到汇流分析的完整链路几乎没人讲透。30米分辨率意味着每个像素代表地面30×30米的区域对省级尺度的地形分析来说这个精度刚好卡在「够用」和「算得动」的平衡点上。这篇文章面向的是需要拿DEM做实际分析的GIS工程师、水文方向的研究生、以及做遥感应用的开发者我会把从数据获取到典型地形因子的提取流程拆开讲包括每一步的参数怎么设、坑在哪、结果怎么验证。2. 广东省30米DEM的数据源选型与获取路径2.1 为什么是30米而不是12.5米或90米选分辨率这件事本质是在精度、数据量和计算成本之间做三角权衡。广东省陆地面积约17.98万平方公里如果用12.5米分辨率的DEM单省栅格数量大约11.5亿个像素一个浮点型GeoTIFF轻松超过4GB做全省汇流分析时内存直接爆炸。90米分辨率虽然数据量小但在珠三角冲积平原这种地势起伏极小的区域90米网格会严重平滑微地形导致水流方向计算出现大面积平行伪影。30米分辨率的来源主要是ASTER GDEM和Copernicus DEM这两套全球数据集。ASTER GDEM V3的空间分辨率是1弧秒约30米覆盖范围南北纬83度之间2009年发布第一版2019年更新到V3修复了大量V2中的异常值。Copernicus DEM的GLO-30产品同样是30米格网基于TanDEM-X雷达干涉测量垂直精度优于ASTER。两套数据在广东省的表现差异主要体现在对比项ASTER GDEM V3Copernicus GLO-30数据来源光学立体像对雷达干涉测量垂直精度±17米RMSE±4米LE90广东省覆盖完整完整水体处理部分水面有异常值水体区域经过编辑数据格式GeoTIFFGeoTIFF / COG获取方式公开分发公开分发从实际使用体验看如果你做的是水文分析优先选Copernicus GLO-30因为它的水体编辑做得更好能减少大量填洼前的预处理工作。如果做的是宏观地形可视化ASTER GDEM V3也够用。2.2 下载与裁剪从全球瓦片到广东省范围无论选哪套数据下载下来的都是按经纬度分幅的瓦片。广东省的经纬度范围大致是东经109°45′到117°20′北纬20°12′到25°31′。按1°×1°分幅需要覆盖的瓦片包括N20E109到N25E117这个范围。下载完成后第一步是镶嵌和裁剪。用GDAL命令行可以完成整个流程# 第一步将下载的多个瓦片合并为一个VRT虚拟栅格 gdalbuildvrt gd_dem.vrt *.tif # 第二步按广东省行政边界裁剪假设已有shapefile gdalwarp -cutline guangdong_boundary.shp \ -crop_to_cutline \ -dstnodata -9999 \ -t_srs EPSG:4326 \ -tr 0.0002777778 0.0002777778 \ -r bilinear \ -of GTiff \ -co COMPRESSLZW \ -co TILEDYES \ gd_dem.vrt gd_dem_30m.tif这里几个参数值得展开说。-tr 0.0002777778 0.0002777778 是把输出分辨率设为1弧秒换算过来就是约30米这个值必须显式指定否则gdalwarp会按源数据分辨率自动推算可能引入微小偏差。-r bilinear 指定重采样方法为双线性插值对高程数据来说比最近邻更平滑但如果你后续要做精确的坡度计算建议改成 -r cubicspline三次样条在保持地形细节上更好。-dstnodata -9999 是给裁剪边界外的区域设一个无效值后续分析时要记得把这个值排除。注意裁剪后的文件用EPSG:4326地理坐标系单位是度。做坡度、坡向、汇流这些分析之前必须投影到等面积或等距投影坐标系否则计算出来的坡度值是错的。广东省常用的投影是CGCS2000 3度带高斯-克吕格投影中央经线114°E对应EPSG:4547。2.3 坐标转换与投影最容易翻车的一步我见过太多人拿着EPSG:4326的DEM直接算坡度结果整个广东省的坡度都偏小因为纬度越高一个经度对应的实际距离越短。在赤道附近1度经度约111公里到了北纬23度只有约102公里差了将近9%。投影转换的命令# 将地理坐标系DEM投影到CGCS2000 3度带 gdalwarp -s_srs EPSG:4326 \ -t_srs EPSG:4547 \ -tr 30 30 \ -r cubicspline \ -of GTiff \ -co COMPRESSLZW \ -co TILEDYES \ -co BIGTIFFIF_SAFER \ gd_dem_30m.tif gd_dem_30m_proj.tif-tr 30 30 把输出像元大小设为30米×30米这是在投影坐标系下的真实地面距离。-co BIGTIFFIF_SAFER 是个保险措施投影后的文件可能超过4GB普通TIFF格式存不下BIGTIFFIF_SAFER会在文件超过4GB时自动切换为BigTIFF格式。转换完成后用gdalinfo检查一下gdalinfo gd_dem_30m_proj.tif | head -30重点看几个字段Pixel Size应该是(30, 30)Coordinate System应该是CGCS2000 / 3-degree Gauss-Kruger zone 38EPSG:4547对应的描述NoData Value应该是-9999。如果Pixel Size显示的是(0.000277..., 0.000277...)说明投影没成功还在地理坐标系下。3. 用GDAL和Rasterio提取坡度坡向与地形因子3.1 坡度计算算法选择与边界处理坡度是DEM最基础也最常用的衍生因子。计算坡度的核心是拟合每个像素与其8个邻域像素构成的三次曲面然后求该曲面在x和y方向的一阶偏导。GDAL的gdaldem工具内置了Horn算法这也是ArcGIS默认使用的算法。# 计算坡度输出单位为度 gdaldem slope gd_dem_30m_proj.tif gd_slope.tif \ -of GTiff \ -compute_edges \ -co COMPRESSLZW \ -co TILEDYES # 计算坡向输出0-360度正北为0顺时针增加 gdaldem aspect gd_dem_30m_proj.tif gd_aspect.tif \ -of GTiff \ -zero_for_flat \ -compute_edges \ -co COMPRESSLZW-compute_edges 参数让GDAL计算边缘像素的坡度值不加这个参数的话图像最外圈一圈像素的坡度会是NoData。-zero_for_flat 让平坦区域坡度小于某个阈值的坡向输出为0而不是-9999在做坡向统计时更友好。坡度结果验证有个简单办法广东省的地形大致是北高南低粤北南岭山区坡度普遍在15°到35°之间珠三角平原坡度多在2°以下雷州半岛台地坡度在3°到8°。如果你的坡度图在珠三角区域出现大片大于10°的值大概率是投影没转对或者DEM里有异常高值没处理。3.2 坡向与地形起伏度从单因子到组合因子坡向计算本身不复杂但坡向的后续应用有个坑坡向是圆形数据0°和360°是同一个方向不能直接做算术平均。如果你要统计某个区域的平均坡向必须用向量合成法——把每个像素的坡向转成单位向量分别求x分量和y分量的均值再反算角度。import numpy as np import rasterio # 读取坡向栅格 with rasterio.open(gd_aspect.tif) as src: aspect src.read(1) nodata src.nodata profile src.profile # 排除NoData valid_mask aspect ! nodata aspect_valid aspect[valid_mask] # 将角度转为弧度 aspect_rad np.deg2rad(aspect_valid) # 向量合成法计算平均坡向 mean_x np.mean(np.cos(aspect_rad)) mean_y np.mean(np.sin(aspect_rad)) mean_aspect np.rad2deg(np.arctan2(mean_y, mean_x)) if mean_aspect 0: mean_aspect 360 print(f区域平均坡向: {mean_aspect:.2f}度)这段代码的关键在np.arctan2的使用它能正确处理所有象限的角度比np.arctan靠谱得多。另外注意坡向的0°在GIS里通常指正北但数学上0弧度指正东所以如果你的坡向数据定义不同需要先做偏移校正。地形起伏度是另一个常用组合因子定义为单位面积内最高海拔与最低海拔之差。在30米分辨率下常用的统计窗口是3×3、5×5和11×11像元。窗口越大起伏度反映的地形尺度越宏观。from scipy.ndimage import maximum_filter, minimum_filter with rasterio.open(gd_dem_30m_proj.tif) as src: dem src.read(1).astype(np.float32) nodata src.nodata dem[dem nodata] np.nan # 3x3窗口起伏度 window_size 3 max_dem maximum_filter(dem, sizewindow_size) min_dem minimum_filter(dem, sizewindow_size) relief_3x3 max_dem - min_dem # 保存结果 profile src.profile profile.update(dtypefloat32, nodatanp.nan) with rasterio.open(gd_relief_3x3.tif, w, **profile) as dst: dst.write(relief_3x3, 1)maximum_filter和minimum_filter来自scipy.ndimage对NaN的处理需要注意——这两个函数会把NaN当作极值参与比较导致结果异常。稳妥的做法是先把NaN替换为一个远小于真实高程的值比如-9999算完再还原。3.3 汇流分析前的填洼处理做水文分析绕不开填洼。DEM中的洼地分两种一种是真实地形中的洼地比如喀斯特溶蚀洼地另一种是数据误差造成的伪洼地。填洼算法会把所有洼地填平这在消除伪洼地的同时也会抹掉真实洼地所以填洼前后一定要做对比验证。# 使用RichDEM进行填洼需要先安装richdem python -c import richdem as rd dem rd.LoadGDAL(gd_dem_30m_proj.tif) dem_filled rd.FillDepressions(dem, epsilonTrue, in_placeFalse) rd.SaveGDAL(gd_dem_filled.tif, dem_filled) epsilonTrue 使用的是Wang Liu算法填洼后在平坦区域加入微小坡度默认1e-5保证水流方向计算时不会出现平行流。如果epsilonFalse填出来的区域是绝对平坦的后续D8算法会在大片平坦区产生平行伪影。填洼前后的差值图能告诉你哪些区域被修改了。如果差值超过5米的区域面积占比超过10%说明原始DEM的伪洼地太多可能需要考虑换一套数据源或者先做一次中值滤波再填洼。4. 广东省30米DEM实操避坑与常见问题排查4.1 坑一NoData值被当成真实高程参与计算现象坡度图在沿海区域出现大片异常高值或者汇流分析时水流从海洋往陆地流。原因DEM裁剪后边界外和海洋区域的NoData值被设为-9999但很多分析工具默认不识别NoData标记直接把-9999当成真实高程参与坡度计算导致-9999与相邻的50米高程之间产生巨大落差坡度直接爆表。解决在做任何分析之前先用gdalwarp或rasterio把NoData区域统一处理。最稳妥的做法是在Python里显式设置import numpy as np import rasterio with rasterio.open(gd_dem_30m_proj.tif) as src: dem src.read(1).astype(np.float32) nodata src.nodata # 显式将NoData替换为NaN dem[dem nodata] np.nan # 后续所有计算基于NaN传播用NaN而不是-9999的好处是NumPy的大多数统计函数会自动忽略NaN不需要每次手动排除。4.2 坑二投影转换后分辨率不是整30米现象gdalinfo查看投影后的DEMPixel Size显示29.9998或30.0003。原因gdalwarp在投影转换时如果-tr参数没有显式指定会按源数据的地理坐标分辨率自动推算投影后的分辨率这个推算过程涉及椭球体参数和投影参数浮点运算必然引入微小误差。解决始终显式指定-tr 30 30。如果已经转换了用gdalwarp再做一次重采样gdalwarp -tr 30 30 -r near gd_dem_30m_proj.tif gd_dem_30m_proj_fix.tif-r near 用最近邻重采样不会改变高程值只调整网格对齐。4.3 坑三ASTER GDEM在水体区域的高程异常现象珠江口、近岸海域的DEM高程显示为正值比如20米、50米而不是0或负值。原因ASTER GDEM基于光学立体像对水面对光学传感器来说几乎是均质表面匹配算法容易失败产生随机高程值。虽然V3版本做了大量修复但广东省沿海的养殖塘、河口浅水区仍然存在异常。解决用最新的水体矢量数据比如从遥感影像提取的年度水体范围对DEM做掩膜处理水体区域高程强制设为0。如果做洪水淹没模拟这一步是必须的。import geopandas as gpd import rasterio from rasterio.mask import mask with rasterio.open(gd_dem_30m_proj.tif) as src: water gpd.read_file(water_mask.shp) # 将水体区域掩膜为0 out_image, out_transform mask(src, water.geometry, cropFalse, filledTrue, nodata0) profile src.profile with rasterio.open(gd_dem_water_fixed.tif, w, **profile) as dst: dst.write(out_image)4.4 坑四填洼后平坦区域面积过大现象填洼后的DEM在平原区域出现大片完全相同的高程值汇流分析结果出现长距离平行线。原因填洼算法把所有洼地填到溢出点高程如果原始DEM在平原区有大量微小洼地高程差只有几厘米填洼后会形成大面积平坦区。D8算法在平坦区无法确定唯一流向只能按固定规则分配产生平行伪影。解决两个办法。一是填洼时加epsilon参数在平坦区引入微小坡度二是改用D-infinity算法它能处理平坦区的多向流分配。RichDEM和WhiteboxTools都支持D-infinity。4.5 坑五大文件处理时的内存溢出现象用Python读取全省30米DEM时程序被系统杀掉或报MemoryError。原因广东省30米DEM投影后约2亿个像素float32格式下约800MB看起来不大但很多操作比如坡度计算、滤波会产生多个中间数组内存占用轻松超过8GB。解决分块处理。rasterio支持窗口读取每次只读一块with rasterio.open(gd_dem_30m_proj.tif) as src: for ji, window in src.block_windows(1): block src.read(1, windowwindow) # 对block做处理 # 处理完写回block_windows会自动按栅格的内部块大小切分通常每个块是256×256或512×512像素内存占用可控。5. 从DEM到汇流网络一个可复现的完整分析链路前面把数据获取、投影转换、因子提取和避坑都讲完了这一章把整条链路串起来用一个具体的目标收尾从广东省30米DEM出发提取全省汇流网络并计算每个子流域的面积。整个流程分五步填洼、计算流向、计算汇流累积量、提取河网、划分子流域。我用WhiteboxTools来做因为它对大数据量的处理比纯Python实现快很多而且命令行调用方便集成到脚本里。# 第一步填洼 whitebox_tools -rFillDepressions -v \ --demgd_dem_30m_proj.tif \ --outputgd_dem_filled.tif # 第二步计算D8流向 whitebox_tools -rD8Pointer -v \ --demgd_dem_filled.tif \ --outputgd_flowdir.tif # 第三步计算汇流累积量 whitebox_tools -rD8FlowAccumulation -v \ --inputgd_flowdir.tif \ --outputgd_flowacc.tif \ --pntr \ --out_typecells # 第四步提取河网阈值设为5000个汇流像元 whitebox_tools -rExtractStreams -v \ --flow_accumgd_flowacc.tif \ --outputgd_streams.tif \ --threshold5000 # 第五步划分子流域 whitebox_tools -rBasins -v \ --d8_pntrgd_flowdir.tif \ --outputgd_basins.tif \ --esri_pntr阈值5000个汇流像元是什么概念每个像元30×30米5000个像元约4.5平方公里。这个阈值下提取的河网密度大致对应广东省中小河流的级别。如果你只关心主要河流把阈值提到50000约45平方公里如果要提取细小沟谷降到1000约0.9平方公里。提取完河网后用Python统计每个子流域的面积和平均坡度import numpy as np import rasterio from rasterio.features import shapes import geopandas as gpd from shapely.geometry import shape # 读取子流域栅格 with rasterio.open(gd_basins.tif) as src: basins src.read(1) transform src.transform # 读取坡度栅格 with rasterio.open(gd_slope.tif) as src_slope: slope src_slope.read(1) # 矢量化子流域 basin_shapes shapes(basins, transformtransform) records [] for geom, value in basin_shapes: if value 0: # 跳过背景值 continue mask basins value mean_slope np.mean(slope[mask]) area_km2 np.sum(mask) * 30 * 30 / 1e6 records.append({ basin_id: int(value), area_km2: round(area_km2, 2), mean_slope_deg: round(float(mean_slope), 2), geometry: shape(geom) }) gdf gpd.GeoDataFrame(records, crsEPSG:4547) gdf.to_file(gd_basins_stats.shp, encodingutf-8) print(f共提取子流域 {len(gdf)} 个) print(f最大流域面积: {gdf[area_km2].max():.2f} km²) print(f最小流域面积: {gdf[area_km2].min():.2f} km²)这段代码里有个细节值得注意np.mean(slope[mask]) 计算的是子流域内的平均坡度但坡度数据里可能有NoData值比如水体区域如果NoData是-9999均值会被严重拉偏。稳妥的做法是在读取坡度后先做一次NoData替换slope slope.astype(np.float32) slope[slope 0] np.nan # 坡度不可能为负 mean_slope np.nanmean(slope[mask])用np.nanmean替代np.mean自动忽略NaN。最后一步验证把提取的河网叠加到DEM山体阴影图上目视检查河网是否沿着山谷走。如果发现河网横穿山脊说明流向计算有问题大概率是填洼不彻底或者投影坐标系搞错了。另一个验证方法是检查河网的源头高程——广东省主要河流的源头大多在粤北山区海拔500米以上如果提取的河网源头出现在珠三角平原那肯定是哪里出了错。我自己做这类分析的习惯是每完成一步就导出一张中间结果图看一眼不要等全流程跑完再检查。填洼后的DEM和原始DEM做个差值图流向图的NoData区域是不是和DEM的NoData区域一致汇流累积量的最大值是不是出现在珠江口——这些中间检查花不了几分钟但能帮你省下大量返工时间。希望帮到你。本文还有配套的精品资源点击获取