ARTICLE DETAIL

资讯详情

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

云南省30米DEM地形分析全流程:从投影转换到坡度起伏度提取

云南省30米DEM地形分析全流程:从投影转换到坡度起伏度提取 简介这份云南省30米分辨率DEM数据面向地理信息、测绘、环境研究与城市规划等领域的从业者和学习者提供覆盖全省的高程地形基础数据可用于洪水风险评估、地质灾害分析、地形地貌研究及交通线路设计等场景。压缩包共10个文件约1MB以GeoTIFF栅格影像为核心辅以shp、shx、dbf、prj等矢量边界文件以及tfw、xml、sbn、sbx等坐标与索引辅助文件采用WGS84坐标系便于在QGIS、ArcGIS等软件中直接加载与空间分析。目前已有1057人学习下载。数据基于ASTER GDEM V3版本精度为30米读者可据此完成云南省地形可视化、坡度坡向提取、流域分析等操作并借助配套边界文件快速裁剪研究区是开展区域地理研究、课程实验与项目建模的实用底图资料。1. 云南省DEM30米分辨率从数据认知到地形分析的落地路径拿到一份省级30米DEM很多人第一反应是“不就是个高程栅格吗”结果一打开ArcGIS就卡在投影不对、范围对不上、nodata值离谱这些破事上。云南省DEM30米分辨率覆盖全省39.4万平方公里从滇西北海拔6740米的卡瓦格博峰到河口县76米的河谷高差超过6600米这种极端地形对数据精度和后续分析流程的要求比平原地区苛刻得多。这份资源适合做水文分析、坡度坡向提取、地形起伏度计算、工程选址、遥感影像正射校正的从业者也适合需要省级尺度地形底图做空间建模的研究人员。它解决的核心问题是给你一套可以直接进入GIS工作流的高程数据省掉从公开源拼接、裁剪、重投影的重复劳动。但能不能用好取决于你对投影、分辨率、nodata和地形分析参数的理解深度。2. 30米DEM的技术底座投影、基准与分辨率选型2.1 为什么是30米而不是12.5米或90米30米分辨率对应的是SRTM航天飞机雷达地形测绘任务和ASTER GDEM的经典格网间距1弧秒约等于30米。这个尺度在省级尺度上是一个平衡点90米太粗做小流域水文分析时河道会断线12.5米更精细但云南省面积大全量数据量会膨胀到几十GB普通工作站处理起来内存吃紧。30米在云南省范围内单波段16位整型压缩后大约几百MB到1GB出头笔记本16GB内存能扛住分块处理。常见做法是省级宏观分析用30米重点区域再用更高分辨率数据补充。选型时还要注意一个坑不同来源的30米DEM垂直精度差异很大。SRTM v3的绝对高程精度约±16米ASTER GDEM v3约±17米但在地形陡峭区误差会放大。云南省滇西北高山峡谷区DEM高程误差在坡度大于30度的区域可能达到20米以上做坡度计算时这个误差会被进一步放大。如果你的应用对高程绝对值敏感比如大坝选址30米DEM只能做初筛不能做最终设计依据。2.2 投影与基准WGS84地理坐标还是UTM投影拿到数据第一件事是确认坐标系。云南省DEM常见分发格式是WGS84地理坐标系EPSG:4326单位是度像元大小0.000277778度。这种格式适合存储和分发但不适合做面积、距离、坡度计算因为经纬度格网在高纬度地区不是正方形。做地形分析前我一般会先投影到适合云南的投影坐标系。云南省跨UTM 47N和48N两个带中央经线分别是99°E和105°E。更推荐用Albers等面积投影自定义参数如下# GDAL 投影转换WGS84 地理坐标 - Albers 等面积投影云南适用 gdalwarp -t_srs projaea lat_125 lat_229 lat_027 lon_0102 x_00 y_00 datumWGS84 unitsm no_defs \ -tr 30 30 \ -r bilinear \ -of GTiff \ -co COMPRESSLZW \ -co TILEDYES \ yunnan_dem_wgs84.tif \ yunnan_dem_albers.tif这段命令做了三件事把地理坐标转成Albers等面积投影重采样到30米×30米像元用LZW压缩并分块存储。-tr 30 30指定输出分辨率-r bilinear是双线性插值适合连续型高程数据。-co TILEDYES让后续按块读取更快不然GDAL每次都要读整行。注意如果原始数据已经是投影坐标不要再转一次否则重采样会引入额外误差。2.3 nodata与无效值处理云南省DEM的nodata值常见是-9999或-32768。如果直接拿去做坡度计算这些值会被当成真实高程参与运算结果就是边界出现巨大异常值。处理方式是在投影转换时同时设nodatagdalwarp -t_srs projaea lat_125 lat_229 lat_027 lon_0102 x_00 y_00 datumWGS84 unitsm no_defs \ -tr 30 30 \ -r bilinear \ -srcnodata -9999 \ -dstnodata -9999 \ -of GTiff \ -co COMPRESSLZW \ yunnan_dem_wgs84.tif \ yunnan_dem_albers.tif-srcnodata告诉GDAL原始数据里哪个值是无效的-dstnodata指定输出无效值。如果原始nodata是-32768就改成-32768。处理完之后用gdalinfo -stats看一眼统计值确认最小值不是-9999否则说明nodata没设对。3. 从DEM到地形因子坡度、坡向与起伏度提取实操3.1 坡度与坡向算法选择与单位陷阱坡度计算有两种常用算法Horn算法和Zevenbergen-Thorne算法。GDAL默认用HornArcGIS默认用Zevenbergen-Thorne。两者在平缓地区差异不大但在陡峭地形区Zevenbergen-Thorne对噪声更敏感。云南省高山峡谷区我一般用Horn结果更平滑。用GDAL算坡度# 坡度计算输出单位为度 gdaldem slope \ -alg Horn \ -compute_edges \ -of GTiff \ -co COMPRESSLZW \ yunnan_dem_albers.tif \ yunnan_slope.tif # 坡向计算输出0-360度正北为0 gdaldem aspect \ -compute_edges \ -of GTiff \ -co COMPRESSLZW \ yunnan_dem_albers.tif \ yunnan_aspect.tif-alg Horn指定算法-compute_edges让边缘像元也能算出值不然边界一圈是nodata。坡度输出默认是度如果你要弧度加-p。坡向输出是0到360度0是正北90是正东。注意坡向在平坦区域坡度接近0没有意义算出来是随机值后续分析时要按坡度阈值掩膜掉。3.2 地形起伏度窗口大小的选择逻辑地形起伏度是指定窗口内最大高程与最小高程之差反映地形破碎程度。窗口大小直接决定结果窗口太小结果噪声大窗口太大起伏度被平滑失去局部特征。云南省地形复杂我一般用3×3窗口做微观起伏用11×11窗口做中观起伏用21×21窗口做宏观起伏。用Python和GDAL计算起伏度from osgeo import gdal import numpy as np from scipy.ndimage import maximum_filter, minimum_filter # 打开DEM ds gdal.Open(yunnan_dem_albers.tif) band ds.GetRasterBand(1) dem band.ReadAsArray().astype(np.float32) nodata band.GetNoDataValue() # 把nodata设为nan避免参与计算 dem[dem nodata] np.nan # 定义窗口大小 window_size 11 # 计算最大值和最小值滤波 max_dem maximum_filter(dem, sizewindow_size) min_dem minimum_filter(dem, sizewindow_size) # 起伏度 最大值 - 最小值 relief max_dem - min_dem # 恢复nodata relief[np.isnan(relief)] nodata # 写出结果 driver gdal.GetDriverByName(GTiff) out_ds driver.Create(yunnan_relief_11x11.tif, ds.RasterXSize, ds.RasterYSize, 1, gdal.GDT_Float32) out_ds.SetGeoTransform(ds.GetGeoTransform()) out_ds.SetProjection(ds.GetProjection()) out_band out_ds.GetRasterBand(1) out_band.WriteArray(relief) out_band.SetNoDataValue(nodata) out_band.FlushCache() out_ds None这段代码用scipy.ndimage的maximum_filter和minimum_filter做窗口统计比手写循环快几个数量级。window_size11对应11×11像元在30米分辨率下约等于330米×330米的窗口。如果你要算21×21改成21就行。注意maximum_filter和minimum_filter对nan的处理是传播nan所以窗口内只要有一个nan输出就是nan这正好符合nodata不参与计算的需求。3.3 水文分析填洼与流向计算DEM做水文分析前必须填洼否则水流方向会被局部凹陷打断。GDAL没有内置填洼常用的是richdem或WhiteboxTools。用richdemimport richdem as rd # 读取DEM dem rd.LoadGDAL(yunnan_dem_albers.tif) # 填洼 dem_filled rd.FillDepressions(dem, epsilonTrue, in_placeFalse) # 计算流向D8算法 flow_dir rd.FlowDirections(dem_filled, methodD8) # 保存结果 rd.SaveGDAL(yunnan_dem_filled.tif, dem_filled) rd.SaveGDAL(yunnan_flow_dir.tif, flow_dir)epsilonTrue表示填洼时保留微小坡度避免大面积平坦区导致流向不确定。methodD8是八方向流向算法每个像元流向八个邻居中坡度最陡的那个。填洼后的DEM再算流向河网才连续。注意richdem对内存要求较高省级DEM建议分块处理或者用WhiteboxTools的FillDepressions它支持分块。4. 避坑与排查30米DEM处理中的五个血泪教训4.1 现象坡度计算结果出现大面积0值或异常高值原因nodata值没设对或者投影转换时重采样方法选错。如果原始nodata是-32768你没设GDAL会把它当真实高程算出来的坡度要么是0因为-32768周围都是-32768要么是巨大值因为-32768和真实高程之间高差几千米。解决用gdalinfo -stats看原始数据的统计值确认nodata。然后在gdalwarp和gdaldem里都显式指定-srcnodata和-dstnodata。如果已经算错了重新从投影转换那一步开始不要试图在坡度结果上修补。4.2 现象投影转换后范围偏移和矢量边界对不上原因原始数据的坐标系定义缺失或错误。有些DEM分发时没有嵌入投影信息GDAL默认按WGS84地理坐标处理但实际可能是其他基准。解决先用gdalinfo看有没有Coordinate System is这一行。如果没有用gdal_edit.py -a_srs EPSG:4326补上。如果补上后还对不上用gdalsrsinfo对比矢量边界的坐标系确认基准是否一致。云南省常见的是WGS84和CGCS2000两者在30米尺度上差异很小但严格来说不能混用。4.3 现象填洼后DEM出现大面积平坦区流向计算失败原因填洼算法把大片区域填成同一高程D8算法在平坦区无法确定流向。解决用richdem的epsilonTrue选项填洼时保留微小坡度。如果已经填平了用rd.FlowDirections的methodD8配合rd.ResolveFlats处理平坦区。或者换用WhiteboxTools的FillDepressions它默认带坡度保持。4.4 现象起伏度计算结果在边界出现异常值原因窗口滤波在边界处窗口不完整maximum_filter和minimum_filter默认用reflect模式填充边界导致边界值失真。解决在计算起伏度前先把DEM边缘裁剪掉窗口半径的宽度或者用modeconstant并指定cvalnp.nan让边界输出nan。我一般用后者然后在后处理时把边界nan掩膜掉。4.5 现象省级DEM处理时内存溢出原因30米云南省DEM全量读入内存约几GB加上中间结果16GB内存容易爆。解决用GDAL的分块读取每次处理一个block。或者用rasterio的block_windowsimport rasterio from rasterio.windows import Window with rasterio.open(yunnan_dem_albers.tif) as src: for ji, window in src.block_windows(1): data src.read(1, windowwindow) # 处理data # 写出到目标文件对应windowblock_windows按GeoTIFF内部的分块大小逐块读取内存占用可控。注意分块处理时窗口滤波需要额外读取边界像元不然块与块之间会出现接缝。常见做法是每个块向外扩窗口半径处理完再裁掉。5. 进阶技巧用DEM做地形阴影与三维可视化验证5.1 地形阴影光照方位角与高度角的参数选择地形阴影hillshade是验证DEM质量最直观的方式。GDAL的gdaldem hillshade默认方位角315度、高度角45度这是制图惯例。但云南省地形复杂默认参数在南北向河谷会显得过暗。我一般用方位角315度、高度角30度让阴影更柔和细节更丰富。gdaldem hillshade \ -az 315 \ -alt 30 \ -z 2 \ -compute_edges \ -of GTiff \ -co COMPRESSLZW \ yunnan_dem_albers.tif \ yunnan_hillshade.tif-z 2是垂直 exaggeration 因子把高程放大2倍让地形起伏更明显。-az 315是光照方位角315度是西北方向符合北半球制图习惯。-alt 30是光照高度角30度比45度更能突出微地形。注意-z不要设太大超过3会让阴影失真看起来像褶皱。5.2 三维可视化用QGIS和Blender快速出图QGIS里加载DEM和hillshade把hillshade放在DEM下面DEM设半透明就能做出立体地形图。如果要更高质量的三维渲染导出DEM为OBJ格式用Blender加光照和材质。# DEM 转 OBJ gdal_translate -of OBJ \ -co FACESYES \ -co FORMATOBJ \ yunnan_dem_albers.tif \ yunnan_dem.objFACESYES生成三角面片FORMATOBJ指定输出格式。导出的OBJ可以直接拖进Blender加一个太阳光调整角度渲染出来就是三维地形。注意省级DEM转OBJ文件会很大建议先降采样或者裁剪重点区域。5.3 验证DEM质量的三个快速检查第一看hillshade有没有明显接缝或条带有的话说明数据拼接有问题。第二算坡度看坡度大于60度的区域占比如果超过5%可能是DEM噪声。第三和已知高程点对比云南省内选几个GNSS控制点看DEM高程和实测高程的差值30米DEM在平坦区差值应在10米以内陡峭区20米以内。从那以后我每次拿到新DEM都强制走一遍“gdalinfo看元数据 → 投影转换 → nodata检查 → hillshade目视 → 坡度统计”这个流程不跳过任何一步。希望帮到你。本文还有配套的精品资源点击获取
返回列表