简介一份适用于地理信息、测绘、环境与城市规划等领域的湖北省30米分辨率DEM数据包源自ASTER GDEM V3影像以GeoTIFF格式存储并附带WGS84坐标系统可直接在ArcGIS、QGIS等软件中加载。栅格影像提供地形高程、坡度、坡向等关键信息可用于地表稳定性评估、洪涝灾害分析、基础设施选址、生态保护及地质灾害风险预测等应用配合附带的湖北省行政边界Shapefile可快速完成省域/市县级裁剪、叠加与分区统计便于区域专题制图和建模。压缩包内共11个文件包括.tif栅格DEM、.shp省界矢量、.prj投影参数、.tfw坐标参考及.dbf属性表等整体大小约233.94MB。已有760人学习/下载适合GIS初学者及需要开展湖北地形分析、科研或项目应用的工程技术人员使用。1. 一份30米DEM文件为什么值得存五年拿到湖北的ASTER GDEM V3数据时很多人第一反应是拖进ArcGIS看一眼山形然后关掉等真正做水文分析或坡度分级时才发现像素尺寸、投影和空值都没确认结果算出来的汇流累积量完全不可用。这份30米分辨率的湖北省DEM原始文件命名里已经透露了它的底细ASTGTMV003是2019年8月发布的全球第三版ASTER高程产品比SRTM 1弧秒更稳定比ALOS AW3D30更省事。对做地质灾害评估、流域划分、基建选址的人来说它够用但要盯住几个坑才能用。2. GeoTIFF与Shapefile湖北省压缩包里的文件分工2.1 .tif与.tfw像素怎么对到地球坐标压缩包里的主文件是HuBei_DEM_30m_ASTGTMV003.tif。GeoTIFF在普通TIFF的基础上内嵌了地理空间元数据GDAL读取后就能识别坐标系和像素尺寸。真正决定“影像的每个像素落在哪里”的是紧随其后的.tfw世界文件。tfw是一个六参数仿射变换文本内容大致是0.0002777778 0 0 -0.0002777778 112.0001388889 33.0001388889六个数依次表示像素宽度、行旋转、列旋转、像素高度负值表示影像从左上角向下排列、左上角X坐标、左上角Y坐标。这里有个容易误读的地方像素宽度和高度是0.0002777778度约等于30米说明该文件采用WGS84经纬度坐标组织而不是UTM投影。很多人在ArcGIS里直接量距离发现“单位是度”问题就出在这里。GeoTIFF内嵌元数据足够时tfw可以省略但ASTER分发的文件习惯保留它因为部分老GIS软件优先读tfw而非内嵌信息。.vat.dbf是把栅格值映射为属性说明的数据库表对连续高程栅格来说通常没有实际内容可以忽略。2.2 那一堆.dbf/.shx/.sbx是什么压缩包里的湖北省.shp是行政边界的矢量数据但Shapefile不是单文件格式而是至少由.shp、.shx、.dbf三件套组成。.shp记录几何坐标.shx是几何索引.dbf存属性字段.sbn和.sbx是ArcGIS生成的空间索引.prj记录坐标系.shp.xml是元数据。这套边界文件的价值在于DEM是按规则网格组织的栅格边界是矢量多边形两者在坐标系一致的前提下才能正确裁剪。处理前先用一个文件管理器确认六个基础文件齐全缺了.sbn不影响读取但缺了.dbf会导致属性字段全部丢失。我一般会在拿到数据后先跑一次ogrinfo确认边界坐标系与DEM一致再决定是否需要重投影。2.3 ASTER GDEM V3为什么它是DEM数据下载的首选V3版本号称对V2做了一次系统性修复核心变化是补掉了大量水体区域高程异常和云层导致的空洞并在低纬度区域替换了新的立体像对数据。对于中国中部省份来说V3的可用性明显优于V2尤其是长江干流两侧的平缓地带V2有时候会出现异常突起的“气泡”V3要干净得多。缺点也有ASTER GDEM的本质是光学立体像对反演的地表高度它描述的是“地表以上首次反射面”不是裸地高程。森林覆盖区和城区会偏高这个特性在第5章详细展开。这也是为什么严谨的处理流程里拿到30米DEM后必须先做质量检查再做裁剪和分析而不是直接出图。3. 裁剪DEM数据的GDAL实操从整片到湖北省界3.1 先看清元数据gdalinfo读坐标系与nodata不管是用QGIS还是纯命令行第一步永远是检查源数据的描述信息。在终端里进入解压目录执行gdalinfo HuBei_DEM_30m_ASTGTMV003.tif输出里重点看四类内容。Coordinate System字段标明坐标系可能显示GEOGCS[WGS 84]也可能显示PROJCS的UTM投影Size是像素行列数结合Pixel Size能估算覆盖范围NoData Value如果没写说明该tif直接用0度或负值表示空值统计时会污染结果最后一行的Corner Coordinates给出四个角点经纬度。ASTER GDEM V3官方分发的tif不统一设置NoData空值区域可能以-9999或0填充这在后续坡度计算中会造成明显的“沟壑”伪影。因此我的处理习惯是第一步就用gdal_translate统一赋空值gdal_translate -a_nodata -9999 -ot Int16 \ HuBei_DEM_30m_ASTGTMV003.tif \ hubei_dem_nodata.tif-a_nodata把-9999写为NoData-ot Int16将高程存储从Float32压到16位有符号整型因为中国东部省份的高程一般不超过9000米Int16足够文件体积直接减半。注意这一步只改元数据不动像素值如果原空值是0而真实地形里也有海拔0米的点请谨慎使用对湖北而言问题不大。3.2 gdalwarp裁剪cutline与crop_to_cutline拿到湖北省.shp后用gdalwarp裁剪这是“裁剪dem数据”的通用做法。不要用ArcGIS的“提取掩膜”因为命令行可复现且不会生成庞大的临时文件gdalwarp -cutline 湖北省.shp \ -crop_to_cutline \ -dstnodata -9999 \ -of GTiff \ hubei_dem_nodata.tif \ hubei_dem_clip.tif参数含义-cutline指定矢量裁剪边界-crop_to_cutline让输出栅格的范围严格贴合多边形外接矩形而不是保留原始数据的完整矩形范围。“湖北省.shp”这里的路径如果有中文在Windows上建议先复制到纯英文目录再执行否则部分GDAL版本会抛出dataset access failed。裁剪完成后必须检查边界处的黑边。如果-dstnodata与源数据nodata不一致边界外的填充黑边会在后续坡度计算中被当作真实高程产生一圈异常陡坡。我习惯在裁剪前先跑每条边的直方图gdalinfo -hist的NoData统计可以快速判断空值占比是否异常。3.3 重投影计算坡度前必须先做投影变换这是很多教程跳过的一步湖北省DEM原始数据如果是WGS84经纬度坐标像素宽度是0.0002778度而纬度方向1度约111公里经度方向1度在湖北纬度约96公里x、y两个方向的单位不等长。直接在这份栅格上算坡度得到的角度是错的因为三角函数的前提是x和y量纲一致。常见做法是转成投影坐标系。湖北位于东经109度至116度跨UTM 49N和50N两个分带直接用单带投影会把省份切成两段。我通常采用自定义中央经线的横轴墨卡托投影或用全国统一的Albers等积圆锥投影。GDAL里临时定义投影gdalwarp -t_srs projtmerc lat_00 lon_0111 k1 x_0500000 y_00 ellpsWGS84 unitsm \ -r cubic \ -dstnodata -9999 \ hubei_dem_clip.tif \ hubei_dem_utm.tiflon_0111把中央经线设在湖北中心附近unitsm保证输出单位为米-r cubic用三次卷积重采样。ASTER的原始30米在重采样后空间分辨率略受影响但对坡度这类派生参数三次卷积比双线性更平滑。如果后续要与其它数据叠加记得同时把矢量边界用ogr2ogr转成同一投影。4. 湖北DEM的地形参数计算与河网提取实战4.1 用gdal.DEMProcessing算坡度与坡向投影确认无误后坡度计算交给GDAL自带的DEMProcessing比ArcGIS的Slope工具更快且结果完全一致。Python调用方式如下from osgeo import gdal gdal.DEMProcessing( hubei_slope.tif, hubei_dem_utm.tif, slope, formatGTiff, slopeFormatdegree, computeEdgesTrue )slopeFormatdegree输出0到90度的坡角如果地面分析需要百分比坡度改成percent。computeEdgesTrue用来补算影像边缘的坡道否则最外圈一圈像素的坡度是NoData镶嵌到更大的图幅里会留下细缝。坡向计算同样一条命令把slope换成aspect即可。坡向输出0到360度其中-1或0表示平地。GDAL的aspect以正北为0度顺时针增加和气象上习惯的方位角一致ArcGIS默认输出也是这个约定。4.2 高程分级与面积统计分析湖北地形时高程分级比裸DEM更直观。我习惯用rasterio做分段统计避免在GIS里反复目估import numpy as np import rasterio with rasterio.open(hubei_dem_utm.tif) as src: dem src.read(1) meta src.profile # 按湖北实际地形定义分级区间 bins [0, 100, 300, 800, 1500, 3000] labels [平原, 丘陵, 低山, 中山, 高山] reclass np.digitize(dem, bins) - 1 reclass np.where(dem 0, 0, reclass) with rasterio.open( hubei_elev_class.tif, w, driverGTiff, heightdem.shape[0], widthdem.shape[1], count1, dtypeuint8, crssrc.crs, transformsrc.transform ) as dst: dst.write(reclass.astype(uint8), 1)高程区间划分没有唯一标准。湖北的平原高程大多在100米以下鄂西山地最高峰神农顶约3105米中间跨越2000米以上的梯度。bins数组的边界设计直接影响分类面积统计给规划部门出图时我一般用[0, 50, 100, 300, 800, 1500]六档能更清楚反映江汉平原的平坦特征。分段结果用uint8存储栅格文件体积极小。分类完成后统计各级面积可以用numpy的bincount再乘以单像素面积。30米分辨率下每像素面积是900平方米如果重采样后分辨率不再是30米要用transform里的a和e参数重新计算像素宽高。4.3 洼地填充与河网提取水文分析是DEM最核心的应用场景。湖北的江汉平原地势极为平坦DEM里充满无数伪洼地——它们不是真实地形而是数据噪声。直接用原始DEM提取河网会得到一组断断续续的乱线必须先填洼。这里提供一个用Python pysheds库的完整流程比在QGIS里点选菜单更适合批量处理和参数复现from pysheds.grid import Grid grid Grid.from_raster(hubei_dem_utm.tif) dem grid.read_raster(hubei_dem_utm.tif) # 第一步填充所有洼地 pit_filled grid.fill_depressions(dem) # 第二步解决平坦区域的流向不定问题 flats grid.resolve_flats(pit_filled) # 第三步D8算法计算流向 dirmap (64, 128, 1, 2, 4, 8, 16, 32) grid.flowdir(flats, out_namedir, dirmapdirmap) # 第四步汇流累积量 grid.accumulation(flats, dirmapdirmap, out_nameacc) acc grid.view(acc) # 第五步阈值提取河网输出为栅格 river acc 300 grid.clip_to(river) grid.to_raster(river, hubei_river.tif)fill_depressions把高程低于周围最低出口的洼地抬高到溢出点高度resolve_flats处理被填平后的平坦区域给微地形加一个单调递增的梯度否则流向计算在平地上会随机乱指。阈值300的物理含义是“上游累积超过300个像素的地区算河道”对应约27万平方米的汇水面积。湖北丘陵地区这个阈值能抓住常年流水河平原区要降到100左右才不会漏掉小支流。流域提取则在river栅格基础上选择出口点利用pysheds的catchment函数输入出口坐标和流向栅格即可得到单一流域范围。5. ASTER GDEM V3的DSM陷阱与空洞修复技巧5.1 先搞清楚它更像DSM还是DEMASTER GDEM V3虽然名字叫DEM但从反演原理上它更接近DSM数字表面模型。光学立体像对匹配的是地物的顶面森林冠层和建筑物屋顶都会拉高高程值。鄂西山区森林覆盖密集同一位置用V3和用LiDAR测得的裸地高程普遍有几米到十几米的偏差坡向朝阴面的森林尤其明显。如果需要的是工程填挖方级别的裸地高程V3的绝对值不能直接用于设计断面如果做的是区域尺度的汇流分析和坡度分级这个精度完全够用。从DSM生成真正意义上的DEM常规做法是用形态学开运算去除凸起地物或配合植被高度模型CHM差分。但ASTER V3的分辨率只有30米单棵树的冠层被混入单个像元形态学滤波的效果有限我一般不建议在V3上做这类处理最多报告里注明“存在植被冠层影响”。5.2 空洞探测与gdal_fillnodata修复V3虽然比V2少了很多空洞但长江沿岸和恩施山区的云层覆盖区仍偶发无值像素。用GDAL自带模块快速定位空洞gdalinfo -stats hubei_dem_clip.tif | grep -E Minimum|Maximum|NoData如果Minimum远小于周边正常高程比如出现-32768说明存在未统一编码的空洞。修复用gdal_fillnodatagdal_fillnodata.py -md 10 -si 0 \ hubei_dem_clip.tif \ hubei_dem_filled.tif-md 10限制搜索半径10个像素距离空洞边缘超过300米的位置不做插值-si 0不启用平滑迭代避免把正常地形抹平。修复完成后重新跑gdalinfo确认Minimum回到合理范围。这一步必须在裁剪和重投影之后做否则插值会跨越边界向外扩。5.3 一个趁手习惯三步gdalinfo帮你排错处理完每个阶段我都会跑一条组合命令检查输出栅格的完整性和投影正确性比打开GIS界面快得多gdalinfo hubei_dem_utm.tif | grep -E Coordinate System is|Size is|Pixel Size|NoData|Minimum|MaximumCoordinate System确认投影参数没丢Size和Pixel Size检查重采样后的行列数是否合理NoData和极值确认没有黑边污染。如果Maximum超过区域最高峰3105米说明有异常噪声未被滤除回到上一步检查原始数据的问题区域。本文还有配套的精品资源点击获取
GDAL处理30米DEM：裁剪、重投影与河网提取实战
简介一份适用于地理信息、测绘、环境与城市规划等领域的湖北省30米分辨率DEM数据包源自ASTER GDEM V3影像以GeoTIFF格式存储并附带WGS84坐标系统可直接在ArcGIS、QGIS等软件中加载。栅格影像提供地形高程、坡度、坡向等关键信息可用于地表稳定性评估、洪涝灾害分析、基础设施选址、生态保护及地质灾害风险预测等应用配合附带的湖北省行政边界Shapefile可快速完成省域/市县级裁剪、叠加与分区统计便于区域专题制图和建模。压缩包内共11个文件包括.tif栅格DEM、.shp省界矢量、.prj投影参数、.tfw坐标参考及.dbf属性表等整体大小约233.94MB。已有760人学习/下载适合GIS初学者及需要开展湖北地形分析、科研或项目应用的工程技术人员使用。1. 一份30米DEM文件为什么值得存五年拿到湖北的ASTER GDEM V3数据时很多人第一反应是拖进ArcGIS看一眼山形然后关掉等真正做水文分析或坡度分级时才发现像素尺寸、投影和空值都没确认结果算出来的汇流累积量完全不可用。这份30米分辨率的湖北省DEM原始文件命名里已经透露了它的底细ASTGTMV003是2019年8月发布的全球第三版ASTER高程产品比SRTM 1弧秒更稳定比ALOS AW3D30更省事。对做地质灾害评估、流域划分、基建选址的人来说它够用但要盯住几个坑才能用。2. GeoTIFF与Shapefile湖北省压缩包里的文件分工2.1 .tif与.tfw像素怎么对到地球坐标压缩包里的主文件是HuBei_DEM_30m_ASTGTMV003.tif。GeoTIFF在普通TIFF的基础上内嵌了地理空间元数据GDAL读取后就能识别坐标系和像素尺寸。真正决定“影像的每个像素落在哪里”的是紧随其后的.tfw世界文件。tfw是一个六参数仿射变换文本内容大致是0.0002777778 0 0 -0.0002777778 112.0001388889 33.0001388889六个数依次表示像素宽度、行旋转、列旋转、像素高度负值表示影像从左上角向下排列、左上角X坐标、左上角Y坐标。这里有个容易误读的地方像素宽度和高度是0.0002777778度约等于30米说明该文件采用WGS84经纬度坐标组织而不是UTM投影。很多人在ArcGIS里直接量距离发现“单位是度”问题就出在这里。GeoTIFF内嵌元数据足够时tfw可以省略但ASTER分发的文件习惯保留它因为部分老GIS软件优先读tfw而非内嵌信息。.vat.dbf是把栅格值映射为属性说明的数据库表对连续高程栅格来说通常没有实际内容可以忽略。2.2 那一堆.dbf/.shx/.sbx是什么压缩包里的湖北省.shp是行政边界的矢量数据但Shapefile不是单文件格式而是至少由.shp、.shx、.dbf三件套组成。.shp记录几何坐标.shx是几何索引.dbf存属性字段.sbn和.sbx是ArcGIS生成的空间索引.prj记录坐标系.shp.xml是元数据。这套边界文件的价值在于DEM是按规则网格组织的栅格边界是矢量多边形两者在坐标系一致的前提下才能正确裁剪。处理前先用一个文件管理器确认六个基础文件齐全缺了.sbn不影响读取但缺了.dbf会导致属性字段全部丢失。我一般会在拿到数据后先跑一次ogrinfo确认边界坐标系与DEM一致再决定是否需要重投影。2.3 ASTER GDEM V3为什么它是DEM数据下载的首选V3版本号称对V2做了一次系统性修复核心变化是补掉了大量水体区域高程异常和云层导致的空洞并在低纬度区域替换了新的立体像对数据。对于中国中部省份来说V3的可用性明显优于V2尤其是长江干流两侧的平缓地带V2有时候会出现异常突起的“气泡”V3要干净得多。缺点也有ASTER GDEM的本质是光学立体像对反演的地表高度它描述的是“地表以上首次反射面”不是裸地高程。森林覆盖区和城区会偏高这个特性在第5章详细展开。这也是为什么严谨的处理流程里拿到30米DEM后必须先做质量检查再做裁剪和分析而不是直接出图。3. 裁剪DEM数据的GDAL实操从整片到湖北省界3.1 先看清元数据gdalinfo读坐标系与nodata不管是用QGIS还是纯命令行第一步永远是检查源数据的描述信息。在终端里进入解压目录执行gdalinfo HuBei_DEM_30m_ASTGTMV003.tif输出里重点看四类内容。Coordinate System字段标明坐标系可能显示GEOGCS[WGS 84]也可能显示PROJCS的UTM投影Size是像素行列数结合Pixel Size能估算覆盖范围NoData Value如果没写说明该tif直接用0度或负值表示空值统计时会污染结果最后一行的Corner Coordinates给出四个角点经纬度。ASTER GDEM V3官方分发的tif不统一设置NoData空值区域可能以-9999或0填充这在后续坡度计算中会造成明显的“沟壑”伪影。因此我的处理习惯是第一步就用gdal_translate统一赋空值gdal_translate -a_nodata -9999 -ot Int16 \ HuBei_DEM_30m_ASTGTMV003.tif \ hubei_dem_nodata.tif-a_nodata把-9999写为NoData-ot Int16将高程存储从Float32压到16位有符号整型因为中国东部省份的高程一般不超过9000米Int16足够文件体积直接减半。注意这一步只改元数据不动像素值如果原空值是0而真实地形里也有海拔0米的点请谨慎使用对湖北而言问题不大。3.2 gdalwarp裁剪cutline与crop_to_cutline拿到湖北省.shp后用gdalwarp裁剪这是“裁剪dem数据”的通用做法。不要用ArcGIS的“提取掩膜”因为命令行可复现且不会生成庞大的临时文件gdalwarp -cutline 湖北省.shp \ -crop_to_cutline \ -dstnodata -9999 \ -of GTiff \ hubei_dem_nodata.tif \ hubei_dem_clip.tif参数含义-cutline指定矢量裁剪边界-crop_to_cutline让输出栅格的范围严格贴合多边形外接矩形而不是保留原始数据的完整矩形范围。“湖北省.shp”这里的路径如果有中文在Windows上建议先复制到纯英文目录再执行否则部分GDAL版本会抛出dataset access failed。裁剪完成后必须检查边界处的黑边。如果-dstnodata与源数据nodata不一致边界外的填充黑边会在后续坡度计算中被当作真实高程产生一圈异常陡坡。我习惯在裁剪前先跑每条边的直方图gdalinfo -hist的NoData统计可以快速判断空值占比是否异常。3.3 重投影计算坡度前必须先做投影变换这是很多教程跳过的一步湖北省DEM原始数据如果是WGS84经纬度坐标像素宽度是0.0002778度而纬度方向1度约111公里经度方向1度在湖北纬度约96公里x、y两个方向的单位不等长。直接在这份栅格上算坡度得到的角度是错的因为三角函数的前提是x和y量纲一致。常见做法是转成投影坐标系。湖北位于东经109度至116度跨UTM 49N和50N两个分带直接用单带投影会把省份切成两段。我通常采用自定义中央经线的横轴墨卡托投影或用全国统一的Albers等积圆锥投影。GDAL里临时定义投影gdalwarp -t_srs projtmerc lat_00 lon_0111 k1 x_0500000 y_00 ellpsWGS84 unitsm \ -r cubic \ -dstnodata -9999 \ hubei_dem_clip.tif \ hubei_dem_utm.tiflon_0111把中央经线设在湖北中心附近unitsm保证输出单位为米-r cubic用三次卷积重采样。ASTER的原始30米在重采样后空间分辨率略受影响但对坡度这类派生参数三次卷积比双线性更平滑。如果后续要与其它数据叠加记得同时把矢量边界用ogr2ogr转成同一投影。4. 湖北DEM的地形参数计算与河网提取实战4.1 用gdal.DEMProcessing算坡度与坡向投影确认无误后坡度计算交给GDAL自带的DEMProcessing比ArcGIS的Slope工具更快且结果完全一致。Python调用方式如下from osgeo import gdal gdal.DEMProcessing( hubei_slope.tif, hubei_dem_utm.tif, slope, formatGTiff, slopeFormatdegree, computeEdgesTrue )slopeFormatdegree输出0到90度的坡角如果地面分析需要百分比坡度改成percent。computeEdgesTrue用来补算影像边缘的坡道否则最外圈一圈像素的坡度是NoData镶嵌到更大的图幅里会留下细缝。坡向计算同样一条命令把slope换成aspect即可。坡向输出0到360度其中-1或0表示平地。GDAL的aspect以正北为0度顺时针增加和气象上习惯的方位角一致ArcGIS默认输出也是这个约定。4.2 高程分级与面积统计分析湖北地形时高程分级比裸DEM更直观。我习惯用rasterio做分段统计避免在GIS里反复目估import numpy as np import rasterio with rasterio.open(hubei_dem_utm.tif) as src: dem src.read(1) meta src.profile # 按湖北实际地形定义分级区间 bins [0, 100, 300, 800, 1500, 3000] labels [平原, 丘陵, 低山, 中山, 高山] reclass np.digitize(dem, bins) - 1 reclass np.where(dem 0, 0, reclass) with rasterio.open( hubei_elev_class.tif, w, driverGTiff, heightdem.shape[0], widthdem.shape[1], count1, dtypeuint8, crssrc.crs, transformsrc.transform ) as dst: dst.write(reclass.astype(uint8), 1)高程区间划分没有唯一标准。湖北的平原高程大多在100米以下鄂西山地最高峰神农顶约3105米中间跨越2000米以上的梯度。bins数组的边界设计直接影响分类面积统计给规划部门出图时我一般用[0, 50, 100, 300, 800, 1500]六档能更清楚反映江汉平原的平坦特征。分段结果用uint8存储栅格文件体积极小。分类完成后统计各级面积可以用numpy的bincount再乘以单像素面积。30米分辨率下每像素面积是900平方米如果重采样后分辨率不再是30米要用transform里的a和e参数重新计算像素宽高。4.3 洼地填充与河网提取水文分析是DEM最核心的应用场景。湖北的江汉平原地势极为平坦DEM里充满无数伪洼地——它们不是真实地形而是数据噪声。直接用原始DEM提取河网会得到一组断断续续的乱线必须先填洼。这里提供一个用Python pysheds库的完整流程比在QGIS里点选菜单更适合批量处理和参数复现from pysheds.grid import Grid grid Grid.from_raster(hubei_dem_utm.tif) dem grid.read_raster(hubei_dem_utm.tif) # 第一步填充所有洼地 pit_filled grid.fill_depressions(dem) # 第二步解决平坦区域的流向不定问题 flats grid.resolve_flats(pit_filled) # 第三步D8算法计算流向 dirmap (64, 128, 1, 2, 4, 8, 16, 32) grid.flowdir(flats, out_namedir, dirmapdirmap) # 第四步汇流累积量 grid.accumulation(flats, dirmapdirmap, out_nameacc) acc grid.view(acc) # 第五步阈值提取河网输出为栅格 river acc 300 grid.clip_to(river) grid.to_raster(river, hubei_river.tif)fill_depressions把高程低于周围最低出口的洼地抬高到溢出点高度resolve_flats处理被填平后的平坦区域给微地形加一个单调递增的梯度否则流向计算在平地上会随机乱指。阈值300的物理含义是“上游累积超过300个像素的地区算河道”对应约27万平方米的汇水面积。湖北丘陵地区这个阈值能抓住常年流水河平原区要降到100左右才不会漏掉小支流。流域提取则在river栅格基础上选择出口点利用pysheds的catchment函数输入出口坐标和流向栅格即可得到单一流域范围。5. ASTER GDEM V3的DSM陷阱与空洞修复技巧5.1 先搞清楚它更像DSM还是DEMASTER GDEM V3虽然名字叫DEM但从反演原理上它更接近DSM数字表面模型。光学立体像对匹配的是地物的顶面森林冠层和建筑物屋顶都会拉高高程值。鄂西山区森林覆盖密集同一位置用V3和用LiDAR测得的裸地高程普遍有几米到十几米的偏差坡向朝阴面的森林尤其明显。如果需要的是工程填挖方级别的裸地高程V3的绝对值不能直接用于设计断面如果做的是区域尺度的汇流分析和坡度分级这个精度完全够用。从DSM生成真正意义上的DEM常规做法是用形态学开运算去除凸起地物或配合植被高度模型CHM差分。但ASTER V3的分辨率只有30米单棵树的冠层被混入单个像元形态学滤波的效果有限我一般不建议在V3上做这类处理最多报告里注明“存在植被冠层影响”。5.2 空洞探测与gdal_fillnodata修复V3虽然比V2少了很多空洞但长江沿岸和恩施山区的云层覆盖区仍偶发无值像素。用GDAL自带模块快速定位空洞gdalinfo -stats hubei_dem_clip.tif | grep -E Minimum|Maximum|NoData如果Minimum远小于周边正常高程比如出现-32768说明存在未统一编码的空洞。修复用gdal_fillnodatagdal_fillnodata.py -md 10 -si 0 \ hubei_dem_clip.tif \ hubei_dem_filled.tif-md 10限制搜索半径10个像素距离空洞边缘超过300米的位置不做插值-si 0不启用平滑迭代避免把正常地形抹平。修复完成后重新跑gdalinfo确认Minimum回到合理范围。这一步必须在裁剪和重投影之后做否则插值会跨越边界向外扩。5.3 一个趁手习惯三步gdalinfo帮你排错处理完每个阶段我都会跑一条组合命令检查输出栅格的完整性和投影正确性比打开GIS界面快得多gdalinfo hubei_dem_utm.tif | grep -E Coordinate System is|Size is|Pixel Size|NoData|Minimum|MaximumCoordinate System确认投影参数没丢Size和Pixel Size检查重采样后的行列数是否合理NoData和极值确认没有黑边污染。如果Maximum超过区域最高峰3105米说明有异常噪声未被滤除回到上一步检查原始数据的问题区域。本文还有配套的精品资源点击获取
