简介贵州省30米分辨率数字高程模型DEM数据包采用GeoTiff格式与WGS84坐标系数据源自ASTER GDEM V3版本适合GIS、遥感、测绘、环境与城市规划等领域的研究者和工程师用于地形分析与空间建模。压缩包共11个文件约265.65MB核心为贵州全省DEM栅格.tif同时包含贵州省行政边界Shapefile.shp及配套的投影、属性、索引与坐标配准文件.prj、.dbf、.shx、.tfw等可直接在主流GIS平台中加载使用。已有416人学习下载。DEM数据经拼接处理覆盖全省并带有明确的坐标系和边界文件用户可直接开展提取坡度坡向、流域分析、洪水模拟、地质灾害评价或山地城镇建设适宜性评估等工作无需再自行拼接和配准大幅节省数据预处理时间。1. 为什么30米分辨率的贵州省DEM仍然值得认真处理30米分辨率的DEM确实不是新技术但放在贵州省这个尺度上它依然是最难被替代的公开数据。省内高精度DEM要么是12.5米的分块存档需要逐片下载要么是不便共享的机载LAS点云。ASTER GDEM V3是经过多次修正的全球数据本次自己拼接后的覆盖文件把整个贵州省落在了一张GeoTIFF里面还附带了省界shp。对于研究喀斯特地貌、水文分析前期准备、流域划分这类工作这份数据足够把大趋势算准。需要提醒的是30米数据在坡度和强度分析时会低估陡坡但用于宏观选址和风险分级仍然适用。下面从数据包的实际文件结构开始逐步拆解预处理、分析、与排错流程。2. 解包与结构GeoTIFF、TFW、aux.xml与贵州省.shp的分工2.1 rar包里的这些文件分别干什么解压“贵州DEM.rar”后看到的不只是两个文件。把每个文件的功能搞明白能避免后续复制数据时漏带同名辅助文件而打不开。文件说明GuiZhou_DEM_30m_ASTGTMV003.tif主数据GeoTIFF格式内含高程值栅格贵州省.shp贵州省省界多边形矢量边界贵州省.dbf / shx / sbn / sbxshapefile的属性表、索引与几何索引贵州省.prj省界shp的坐标系统描述GuiZhou_DEM_30m_ASTGTMV003.tfwTIFF的世界文件记录仿射变换参数GuiZhou_DEM_30m_ASTGTMV003.tif.aux.xmlGDAL生成的辅助信息含NODATA等GuiZhou_DEM_30m_ASTGTMV003.tif.vat.dbf栅格属性表通常由第一次统计生成很多人会把.tfw忽略掉因为GeoTIFF内部已经写了地理参考。实际上.tfw是给不支持内部几何的软件用的例如老版本图片查看器或某些CAD工具。如果这两份文件存在且内容一致说明当初拼接时已经同步写出。aux.xml里最值得关注的是NODATA值后续做坡度或填洼前必须知道这个值否则边界部分会算出一圈极端结果。2.2 用gdalinfo确认数据真实范围、分辨率与坐标系拿到数据第一件事不是直接开QGIS而是用gdalinfo把元数据读出来。无论你用的是Windows下的OSGeo4W Shell还是Linux里的GDAL环境这个命令都是标配。gdalinfo GuiZhou_DEM_30m_ASTGTMV003.tif运行后会看到类似下面的关键块Driver: GTiff/GeoTIFF Files: GuiZhou_DEM_30m_ASTGTMV003.tif GuiZhou_DEM_30m_ASTGTMV003.tfw GuiZhou_DEM_30m_ASTGTMV003.tif.aux.xml Size is 4841, 4536 Coordinate System is: GEOGCRS[WGS 84, ...] Data axis to CRS axis mapping: 2,1 Origin (103.5691667000,29.2568056000) Pixel Size (0.000277777800000,-0.000277777800000)输出里最容易被误读的是Pixel Size (0.0002777778, -0.0002777778)单位是度不是米。0.0002777778度约等于1/3600度在赤道附近对应约30.8米在贵州纬度27°左右东西方向实际距离约为31米到30米之间所以称为“30米分辨率”是行业内习惯说法。严格来说这种地理坐标系下的像元并不是正方形米制格子如果后续要做面积和坡度计算最好先投影到适合贵州的UTM 48N或Albers等积投影。检查完分辨率后还要看范围Origin和Size决定了数据的覆盖框。贵州的经纬度范围大致是东经103°36′到109°35′北纬24°37′到29°13′如果gdalinfo的Extent与此相差过大说明拼接时边缘没有裁干净。2.3 用ogrinfo检查省界属性并确认与DEM叠不齐矢量部分用ogrinfo检查重点看prj和属性字段不要只看形状。ogrinfo -so -al 贵州省.shp输出的关键部分包括Layer name: 贵州省 Geometry: Polygon Feature Count: 1 Extent: (103.567310, 24.591955) - (109.374590, 29.244900) Layer SRS: WGS_1984_UTM_Zone_48N id: 1 name: 贵州省如果prj里写的是WGS_1984_UTM_Zone_48N而DEM是WGS 84那么两者叠加后经纬度坐标下边界会出现在正确位置但当你用米制单位量算距离时DEM还没有投影叠加结果会出现比例错位。拿到这份数据后建议将省界和DEM统一转换到同一投影再做裁剪或面积统计。通常做法是先复制一份DEM到UTM 48N投影然后再用省界裁剪。3. 坡度、坡向与山体阴影自写脚本和gdaldem的对比3.1 坡度是怎么算出来的差分窗口与地形因子对DEM做地形分析基础是你得明白算法在像元邻域上怎么算。坡度通常用二阶差分法对3x3窗口计中心像元在东西和南北方向的高程变化率dZ/dX与dZ/dY坡度就是这两个方向梯度的反正切合成。坡向则是该梯度方向在平面上的方位角从正北起顺时针计算。这里的核心参数有两个一个是高程缩放因子常用1.0第二个是水平距离单位。如果DEM是地理坐标那么dX和dY的单位是度必须乘上纬度对应的弧长系数才能得到米制坡度。全程经验做法是先把数据投影到UTM让水平和垂直单位统一成米省去换算坡度计算也更稳。ASTER GDEM V3在平坦坝区仍有一些随机噪声计算坡度前可以先做一次3x3中值滤波否则算出的平均坡度在缓坡区域会比真值偏大。3.2 命令行直接生成这几种数据GDAL自带的gdaldem是工作效率最高的方案不需要自己写循环。下面这条命令生成坡度单位选“度”更适合直接看图。gdaldem slope GuiZhou_DEM_30m_ASTGTMV003.tif 贵州坡度.tif \ -p \ -s 111320.0 \ -compute_edges命令中-p表示输出单位为度若不加则输出百分数坡度-s是垂直比例当地理坐标系下水平单位是度时111320表示1度约等于111320米但这是赤道附近的值。贵州省纬度接近北纬27度严格应按111320*cos(27°)≈99140来设置不过多数情况下用111320会影响几位小数影响有限。-compute_edges让边缘像元也能参与计算而不是输出NoData。如果先把DEM投影成UTM那么-s直接用1.0即可。坡向与山体阴影分别对应gdaldem aspect和gdaldem hillshade后者还接受-az 315 -alt 45控制太阳方位和高度角。3.3 用NumPy自写坡度和坡向方便批量修参数自己写脚本的优势是能随时加滤波或修改邻域权重适合批量测试参数。下面的Python代码将GeoTIFF读成NumPy数组用np.gradient计算梯度然后合成坡度和坡向。import numpy as np from osgeo import gdal src gdal.Open(GuiZhou_DEM_30m_ASTGTMV003.tif) band src.GetRasterBand(1) dem band.ReadAsArray().astype(np.float32) nodata band.GetNoDataValue() if nodata is not None: dem[dem nodata] np.nan gt src.GetGeoTransform() # 像元宽度/高度地理坐标系下单位是度 pixel_width gt[1] pixel_height abs(gt[5]) # 纬度方向做距离校正取数据集中间纬度 center_lat np.deg2rad((gt[3] (gt[3] gt[5] * src.RasterYSize)) / 2) x_scale 111320.0 * np.cos(center_lat) y_scale 111320.0 dzdx, dzdy np.gradient(dem, pixel_width * 1.0, pixel_height * 1.0) # 将度单位梯度转为米单位 dzdx dzdx * pixel_width / x_scale dzdy dzdy * pixel_height / y_scale slope np.degrees(np.arctan(np.sqrt(dzdx**2 dzdy**2))) aspect np.degrees(np.arctan2(-dzdx, dzdy)) aspect np.where(aspect 0, aspect 360, aspect)这段代码中np.gradient(dem, hx, hy)会返回两个方向上的一阶差分。将dx以度为单位换算成米制时乘上pixel_width / x_scale这个比例恰好是每个像元对应的米距。arctan2计算坡向时用(-dzdx, dzdy)是约定坡向指向下坡方向东向为起点逆时针。如果你习惯用0度表示北还要再做一次(450 - np.degrees(...)) % 360的换算两种约定各有拥趸关键是要与后续分析工具的约定保持一致。GDAL里gdaldem aspect的输出是0度代表北采用这个约定最省事。写完坡向数组后建议和gdaldem aspect输出的结果做逐像元差值统计两者差应不超过浮点误差若差很多多半是坐标系比例参数设错了。3.4 技术选择建议工具适用场景主要参数gdaldem slope快速出图不涉及二次处理-p, -s, -compute_edgesNumPy自写需要预处理滤波或定制邻域权重pixel_width, x_scale, y_scaleQGIS地形分析在GUI里往返验证坡度单位、Z factor在贵州省这种地形起伏大的地方Z factor的影响比平缓地区更明显。若直接用地理坐标且不设-s计算出的坡度会异常偏小到几乎平的地。这是初学者最容易踩的坑。4. 水文分析填洼、流向与河网提取的完整流程4.1 为什么必须先进行填洼处理然后再分析流向未经填洼的DEM里包含许多伪洼地这些可能是真实洼地也可能是ASTER GDEM V3的噪声或拼接接缝造成的误差。如果不填洼D8单流向算法会把水流截断在洼地内部后续汇流累积量会被人为截断产生的河网会出现大量细碎断点。填洼并不是将所有洼地都填平而是把高程低于其流出路径最低点的凹陷抬升到最低流出点高程。实际操作中我一般会保留真实喀斯特漏斗区域否则会改变天然消水路径。但这里我们做全省范围的河网提取填洼是必要且合理的。常用的工具包括WhiteboxTools的FillDepressions、GRASS的r.fill.dir以及SAGA。WhiteboxTools对30米ASTER GDEM处理速度较快而且命令行与GDAL一样容易自动化。4.2 投影、填洼与流向提取命令先用GDAL把WGS84地理坐标转成UTM 48N投影这一步必须做。因为水文分析中的流向和累计量都是以米为空间尺度如果保持经纬度坐标汇流累积结果会因纬度不同而随位置变化失真。gdalwarp -t_srs EPSG:32648 -r bilinear -co COMPRESSDEFLATE \ GuiZhou_DEM_30m_ASTGTMV003.tif 贵州DEM_UTM48N.tif然后使用WhiteboxToolswhitebox_tools -rFillDepressions -v \ --input贵州DEM_UTM48N.tif \ --output贵州DEM_fill.tif \ --fix_flatsfalse填洼后可生成D8流向whitebox_tools -rD8Pointer -v \ --input贵州DEM_fill.tif \ --output贵州D8流向.tif \ --esri_pntrfalse最后计算汇流累积量whitebox_tools -rD8FlowAccumulation -v \ --input贵州DEM_fill.tif \ --pntr贵州D8流向.tif \ --output贵州累积量.tif \ --out_typeSpecific Contributing Area命令中--out_typeSpecific Contributing Area输出的单宽汇流面积单位是米便于科学意义解释如果要提取河网通常用累积量而非汇流面积。4.3 河网阈值选择提取河流时需要把汇流累积量栅格按阈值做重分类。阈值太低河网密到什么沟都算河阈值太高主干支流又连不起来。这没有绝对标准。对30米分辨率贵州地区常见阈值范围是1000到5000个像元。可以先试数据范围的0.001%到0.005%。用GDAL计算累积量的统计值再用SQL或? 直接生成矢量化河网。gdallocationinfo 贵州累积量.tif -valonly -geoloc 105 27查询贵州中部某位置的值可以快速判断量级。另一种做法是直接用WhiteboxTools的ExtractStreamswhitebox_tools -rExtractStreams -v \ --flow_accum贵州累积量.tif \ --threshold1500 \ --output贵州河网.tifthreshold1500意味着汇流累积量超过1500个像元的栅格将被划分为河网。接下来可用r.stream.to.vector或gdal_polygonize将河网栅格转为矢量再与贵州省shp叠合。4.4 河网与省界的叠加和可能遇到的问题河网提取出来后建议导入QGIS做目视检查重点看贵州西部和南部的高原边缘这些区域地形切割剧烈更容易出现平行于格网方向的假河流。ASTER GDEM V3的原始数据在陡崖区域仍可能带有条带条带位置会形成异常的直线沟槽。遇到这种情况不要急着调阈值先回到填洼前的DEM用gdaldem hillshade检查条带分布再决定是否需要对原始DEM做低通滤波。对喀斯特区域可对比已有水系里的暗河出口判断提取结果的合理性。5. 验证与排坑NODATA、条带噪声与坐标系误用的三个经验在使用这份贵州省DEM时有三处最容易让结果“看起来合理但实际出错”的环节。第一是NODATA的位置。拼接后的全省DEM边界外应该全是NODATA。运行gdalinfo并查看NoData Value需要确认这个值不是0。如果拼接时没设置NODATA文件外的高程值可能被当作0米那在坡度、填洼分析中会造成边界一圈巨大错误。可用下面命令快速将NODATA重设为原设定值gdal_translate -a_nodata -9999 \ GuiZhou_DEM_30m_ASTGTMV003.tif 贵州DEM_clean.tif执行后再用gdalinfo -stats看统计范围。贵州省海拔最低处大约105米最高海拔约2900米如果你的最大值跑到5000米或出现负值说明残留了无效像元需要先按掩膜或高程合理区间过滤。第二是条带噪声。ASTER GDEM V3虽然比V2好很多但个别区域尤其是南部潮湿区仍有横向或纵向条带。识别方法很简单用gdaldem hillshade生成山体阴影然后把图像拉伸到频率域或直接用目视“横向扫描”便能发现局部条纹。若分析区域恰好被条纹覆盖建议对该局部区域做3x3中值滤波不要全盘处理以免削弱真实地形。3x3中值滤波的副作用是损失部分真实细节30米分辨率下对一个县做宏观分析可接受。第三是坐标系统一。省界shp的prj若是UTM 48N而DEM仍为WGS84度坐标原则上QGIS动态投影能显示正确但做距离和面积计算时必须保证两端在同一投影坐标系。验证方法是对DEM执行一次gdalwarp再用gdalinfo看Corners坐标和shp范围做差值确保边界偏差在1个像元以内。按下面命令输出来检查gdalwarp -t_srs EPSG:32648 -r bilinear -cutline 贵州省.shp \ GuiZhou_DEM_30m_ASTGTMV003.tif 贵州DEM_clip_UTM48N.tif \ -crop_to_cutline -co COMPRESSDEFLATE这条命令同时完成投影、按省界裁剪和NODATA保留。若输出的面积像元数落在预期范围内再开始坡度或汇流分析才能保证后续结果是可靠的。本文还有配套的精品资源点击获取
贵州省30米DEM预处理与地形水文分析完整指南
简介贵州省30米分辨率数字高程模型DEM数据包采用GeoTiff格式与WGS84坐标系数据源自ASTER GDEM V3版本适合GIS、遥感、测绘、环境与城市规划等领域的研究者和工程师用于地形分析与空间建模。压缩包共11个文件约265.65MB核心为贵州全省DEM栅格.tif同时包含贵州省行政边界Shapefile.shp及配套的投影、属性、索引与坐标配准文件.prj、.dbf、.shx、.tfw等可直接在主流GIS平台中加载使用。已有416人学习下载。DEM数据经拼接处理覆盖全省并带有明确的坐标系和边界文件用户可直接开展提取坡度坡向、流域分析、洪水模拟、地质灾害评价或山地城镇建设适宜性评估等工作无需再自行拼接和配准大幅节省数据预处理时间。1. 为什么30米分辨率的贵州省DEM仍然值得认真处理30米分辨率的DEM确实不是新技术但放在贵州省这个尺度上它依然是最难被替代的公开数据。省内高精度DEM要么是12.5米的分块存档需要逐片下载要么是不便共享的机载LAS点云。ASTER GDEM V3是经过多次修正的全球数据本次自己拼接后的覆盖文件把整个贵州省落在了一张GeoTIFF里面还附带了省界shp。对于研究喀斯特地貌、水文分析前期准备、流域划分这类工作这份数据足够把大趋势算准。需要提醒的是30米数据在坡度和强度分析时会低估陡坡但用于宏观选址和风险分级仍然适用。下面从数据包的实际文件结构开始逐步拆解预处理、分析、与排错流程。2. 解包与结构GeoTIFF、TFW、aux.xml与贵州省.shp的分工2.1 rar包里的这些文件分别干什么解压“贵州DEM.rar”后看到的不只是两个文件。把每个文件的功能搞明白能避免后续复制数据时漏带同名辅助文件而打不开。文件说明GuiZhou_DEM_30m_ASTGTMV003.tif主数据GeoTIFF格式内含高程值栅格贵州省.shp贵州省省界多边形矢量边界贵州省.dbf / shx / sbn / sbxshapefile的属性表、索引与几何索引贵州省.prj省界shp的坐标系统描述GuiZhou_DEM_30m_ASTGTMV003.tfwTIFF的世界文件记录仿射变换参数GuiZhou_DEM_30m_ASTGTMV003.tif.aux.xmlGDAL生成的辅助信息含NODATA等GuiZhou_DEM_30m_ASTGTMV003.tif.vat.dbf栅格属性表通常由第一次统计生成很多人会把.tfw忽略掉因为GeoTIFF内部已经写了地理参考。实际上.tfw是给不支持内部几何的软件用的例如老版本图片查看器或某些CAD工具。如果这两份文件存在且内容一致说明当初拼接时已经同步写出。aux.xml里最值得关注的是NODATA值后续做坡度或填洼前必须知道这个值否则边界部分会算出一圈极端结果。2.2 用gdalinfo确认数据真实范围、分辨率与坐标系拿到数据第一件事不是直接开QGIS而是用gdalinfo把元数据读出来。无论你用的是Windows下的OSGeo4W Shell还是Linux里的GDAL环境这个命令都是标配。gdalinfo GuiZhou_DEM_30m_ASTGTMV003.tif运行后会看到类似下面的关键块Driver: GTiff/GeoTIFF Files: GuiZhou_DEM_30m_ASTGTMV003.tif GuiZhou_DEM_30m_ASTGTMV003.tfw GuiZhou_DEM_30m_ASTGTMV003.tif.aux.xml Size is 4841, 4536 Coordinate System is: GEOGCRS[WGS 84, ...] Data axis to CRS axis mapping: 2,1 Origin (103.5691667000,29.2568056000) Pixel Size (0.000277777800000,-0.000277777800000)输出里最容易被误读的是Pixel Size (0.0002777778, -0.0002777778)单位是度不是米。0.0002777778度约等于1/3600度在赤道附近对应约30.8米在贵州纬度27°左右东西方向实际距离约为31米到30米之间所以称为“30米分辨率”是行业内习惯说法。严格来说这种地理坐标系下的像元并不是正方形米制格子如果后续要做面积和坡度计算最好先投影到适合贵州的UTM 48N或Albers等积投影。检查完分辨率后还要看范围Origin和Size决定了数据的覆盖框。贵州的经纬度范围大致是东经103°36′到109°35′北纬24°37′到29°13′如果gdalinfo的Extent与此相差过大说明拼接时边缘没有裁干净。2.3 用ogrinfo检查省界属性并确认与DEM叠不齐矢量部分用ogrinfo检查重点看prj和属性字段不要只看形状。ogrinfo -so -al 贵州省.shp输出的关键部分包括Layer name: 贵州省 Geometry: Polygon Feature Count: 1 Extent: (103.567310, 24.591955) - (109.374590, 29.244900) Layer SRS: WGS_1984_UTM_Zone_48N id: 1 name: 贵州省如果prj里写的是WGS_1984_UTM_Zone_48N而DEM是WGS 84那么两者叠加后经纬度坐标下边界会出现在正确位置但当你用米制单位量算距离时DEM还没有投影叠加结果会出现比例错位。拿到这份数据后建议将省界和DEM统一转换到同一投影再做裁剪或面积统计。通常做法是先复制一份DEM到UTM 48N投影然后再用省界裁剪。3. 坡度、坡向与山体阴影自写脚本和gdaldem的对比3.1 坡度是怎么算出来的差分窗口与地形因子对DEM做地形分析基础是你得明白算法在像元邻域上怎么算。坡度通常用二阶差分法对3x3窗口计中心像元在东西和南北方向的高程变化率dZ/dX与dZ/dY坡度就是这两个方向梯度的反正切合成。坡向则是该梯度方向在平面上的方位角从正北起顺时针计算。这里的核心参数有两个一个是高程缩放因子常用1.0第二个是水平距离单位。如果DEM是地理坐标那么dX和dY的单位是度必须乘上纬度对应的弧长系数才能得到米制坡度。全程经验做法是先把数据投影到UTM让水平和垂直单位统一成米省去换算坡度计算也更稳。ASTER GDEM V3在平坦坝区仍有一些随机噪声计算坡度前可以先做一次3x3中值滤波否则算出的平均坡度在缓坡区域会比真值偏大。3.2 命令行直接生成这几种数据GDAL自带的gdaldem是工作效率最高的方案不需要自己写循环。下面这条命令生成坡度单位选“度”更适合直接看图。gdaldem slope GuiZhou_DEM_30m_ASTGTMV003.tif 贵州坡度.tif \ -p \ -s 111320.0 \ -compute_edges命令中-p表示输出单位为度若不加则输出百分数坡度-s是垂直比例当地理坐标系下水平单位是度时111320表示1度约等于111320米但这是赤道附近的值。贵州省纬度接近北纬27度严格应按111320*cos(27°)≈99140来设置不过多数情况下用111320会影响几位小数影响有限。-compute_edges让边缘像元也能参与计算而不是输出NoData。如果先把DEM投影成UTM那么-s直接用1.0即可。坡向与山体阴影分别对应gdaldem aspect和gdaldem hillshade后者还接受-az 315 -alt 45控制太阳方位和高度角。3.3 用NumPy自写坡度和坡向方便批量修参数自己写脚本的优势是能随时加滤波或修改邻域权重适合批量测试参数。下面的Python代码将GeoTIFF读成NumPy数组用np.gradient计算梯度然后合成坡度和坡向。import numpy as np from osgeo import gdal src gdal.Open(GuiZhou_DEM_30m_ASTGTMV003.tif) band src.GetRasterBand(1) dem band.ReadAsArray().astype(np.float32) nodata band.GetNoDataValue() if nodata is not None: dem[dem nodata] np.nan gt src.GetGeoTransform() # 像元宽度/高度地理坐标系下单位是度 pixel_width gt[1] pixel_height abs(gt[5]) # 纬度方向做距离校正取数据集中间纬度 center_lat np.deg2rad((gt[3] (gt[3] gt[5] * src.RasterYSize)) / 2) x_scale 111320.0 * np.cos(center_lat) y_scale 111320.0 dzdx, dzdy np.gradient(dem, pixel_width * 1.0, pixel_height * 1.0) # 将度单位梯度转为米单位 dzdx dzdx * pixel_width / x_scale dzdy dzdy * pixel_height / y_scale slope np.degrees(np.arctan(np.sqrt(dzdx**2 dzdy**2))) aspect np.degrees(np.arctan2(-dzdx, dzdy)) aspect np.where(aspect 0, aspect 360, aspect)这段代码中np.gradient(dem, hx, hy)会返回两个方向上的一阶差分。将dx以度为单位换算成米制时乘上pixel_width / x_scale这个比例恰好是每个像元对应的米距。arctan2计算坡向时用(-dzdx, dzdy)是约定坡向指向下坡方向东向为起点逆时针。如果你习惯用0度表示北还要再做一次(450 - np.degrees(...)) % 360的换算两种约定各有拥趸关键是要与后续分析工具的约定保持一致。GDAL里gdaldem aspect的输出是0度代表北采用这个约定最省事。写完坡向数组后建议和gdaldem aspect输出的结果做逐像元差值统计两者差应不超过浮点误差若差很多多半是坐标系比例参数设错了。3.4 技术选择建议工具适用场景主要参数gdaldem slope快速出图不涉及二次处理-p, -s, -compute_edgesNumPy自写需要预处理滤波或定制邻域权重pixel_width, x_scale, y_scaleQGIS地形分析在GUI里往返验证坡度单位、Z factor在贵州省这种地形起伏大的地方Z factor的影响比平缓地区更明显。若直接用地理坐标且不设-s计算出的坡度会异常偏小到几乎平的地。这是初学者最容易踩的坑。4. 水文分析填洼、流向与河网提取的完整流程4.1 为什么必须先进行填洼处理然后再分析流向未经填洼的DEM里包含许多伪洼地这些可能是真实洼地也可能是ASTER GDEM V3的噪声或拼接接缝造成的误差。如果不填洼D8单流向算法会把水流截断在洼地内部后续汇流累积量会被人为截断产生的河网会出现大量细碎断点。填洼并不是将所有洼地都填平而是把高程低于其流出路径最低点的凹陷抬升到最低流出点高程。实际操作中我一般会保留真实喀斯特漏斗区域否则会改变天然消水路径。但这里我们做全省范围的河网提取填洼是必要且合理的。常用的工具包括WhiteboxTools的FillDepressions、GRASS的r.fill.dir以及SAGA。WhiteboxTools对30米ASTER GDEM处理速度较快而且命令行与GDAL一样容易自动化。4.2 投影、填洼与流向提取命令先用GDAL把WGS84地理坐标转成UTM 48N投影这一步必须做。因为水文分析中的流向和累计量都是以米为空间尺度如果保持经纬度坐标汇流累积结果会因纬度不同而随位置变化失真。gdalwarp -t_srs EPSG:32648 -r bilinear -co COMPRESSDEFLATE \ GuiZhou_DEM_30m_ASTGTMV003.tif 贵州DEM_UTM48N.tif然后使用WhiteboxToolswhitebox_tools -rFillDepressions -v \ --input贵州DEM_UTM48N.tif \ --output贵州DEM_fill.tif \ --fix_flatsfalse填洼后可生成D8流向whitebox_tools -rD8Pointer -v \ --input贵州DEM_fill.tif \ --output贵州D8流向.tif \ --esri_pntrfalse最后计算汇流累积量whitebox_tools -rD8FlowAccumulation -v \ --input贵州DEM_fill.tif \ --pntr贵州D8流向.tif \ --output贵州累积量.tif \ --out_typeSpecific Contributing Area命令中--out_typeSpecific Contributing Area输出的单宽汇流面积单位是米便于科学意义解释如果要提取河网通常用累积量而非汇流面积。4.3 河网阈值选择提取河流时需要把汇流累积量栅格按阈值做重分类。阈值太低河网密到什么沟都算河阈值太高主干支流又连不起来。这没有绝对标准。对30米分辨率贵州地区常见阈值范围是1000到5000个像元。可以先试数据范围的0.001%到0.005%。用GDAL计算累积量的统计值再用SQL或? 直接生成矢量化河网。gdallocationinfo 贵州累积量.tif -valonly -geoloc 105 27查询贵州中部某位置的值可以快速判断量级。另一种做法是直接用WhiteboxTools的ExtractStreamswhitebox_tools -rExtractStreams -v \ --flow_accum贵州累积量.tif \ --threshold1500 \ --output贵州河网.tifthreshold1500意味着汇流累积量超过1500个像元的栅格将被划分为河网。接下来可用r.stream.to.vector或gdal_polygonize将河网栅格转为矢量再与贵州省shp叠合。4.4 河网与省界的叠加和可能遇到的问题河网提取出来后建议导入QGIS做目视检查重点看贵州西部和南部的高原边缘这些区域地形切割剧烈更容易出现平行于格网方向的假河流。ASTER GDEM V3的原始数据在陡崖区域仍可能带有条带条带位置会形成异常的直线沟槽。遇到这种情况不要急着调阈值先回到填洼前的DEM用gdaldem hillshade检查条带分布再决定是否需要对原始DEM做低通滤波。对喀斯特区域可对比已有水系里的暗河出口判断提取结果的合理性。5. 验证与排坑NODATA、条带噪声与坐标系误用的三个经验在使用这份贵州省DEM时有三处最容易让结果“看起来合理但实际出错”的环节。第一是NODATA的位置。拼接后的全省DEM边界外应该全是NODATA。运行gdalinfo并查看NoData Value需要确认这个值不是0。如果拼接时没设置NODATA文件外的高程值可能被当作0米那在坡度、填洼分析中会造成边界一圈巨大错误。可用下面命令快速将NODATA重设为原设定值gdal_translate -a_nodata -9999 \ GuiZhou_DEM_30m_ASTGTMV003.tif 贵州DEM_clean.tif执行后再用gdalinfo -stats看统计范围。贵州省海拔最低处大约105米最高海拔约2900米如果你的最大值跑到5000米或出现负值说明残留了无效像元需要先按掩膜或高程合理区间过滤。第二是条带噪声。ASTER GDEM V3虽然比V2好很多但个别区域尤其是南部潮湿区仍有横向或纵向条带。识别方法很简单用gdaldem hillshade生成山体阴影然后把图像拉伸到频率域或直接用目视“横向扫描”便能发现局部条纹。若分析区域恰好被条纹覆盖建议对该局部区域做3x3中值滤波不要全盘处理以免削弱真实地形。3x3中值滤波的副作用是损失部分真实细节30米分辨率下对一个县做宏观分析可接受。第三是坐标系统一。省界shp的prj若是UTM 48N而DEM仍为WGS84度坐标原则上QGIS动态投影能显示正确但做距离和面积计算时必须保证两端在同一投影坐标系。验证方法是对DEM执行一次gdalwarp再用gdalinfo看Corners坐标和shp范围做差值确保边界偏差在1个像元以内。按下面命令输出来检查gdalwarp -t_srs EPSG:32648 -r bilinear -cutline 贵州省.shp \ GuiZhou_DEM_30m_ASTGTMV003.tif 贵州DEM_clip_UTM48N.tif \ -crop_to_cutline -co COMPRESSDEFLATE这条命令同时完成投影、按省界裁剪和NODATA保留。若输出的面积像元数落在预期范围内再开始坡度或汇流分析才能保证后续结果是可靠的。本文还有配套的精品资源点击获取
