ARTICLE DETAIL

资讯详情

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

洞庭湖流域DEM数据处理与水文分析全流程指南

洞庭湖流域DEM数据处理与水文分析全流程指南 简介面向地理信息分析与水文研究工作者内含洞庭湖流域30米分辨率数字高程模型DEM栅格数据可支撑坡度与流向计算、汇流区提取、洪水风险区划等宏观地形分析适合GIS学习者、科研人员及规划人员使用。压缩包共6个文件以tif栅格为主体配套ovr金字塔、tfw配准信息、xml元数据、dbf属性表及cpg字符集文件整体约347.9MB数据结构完整便于在ArcGIS、QGIS中直接加载与处理。已有448人学习/下载。借助该数据可完成流域洼地分析、水文建模练习也可作为环境评价、地质灾害研判和土地利用规划的地形底图30米分辨率兼顾精度与体量既适合初学者理解DEM处理流程也能满足区域尺度的科研教学与实际应用需求。1. 洞庭湖流域DEM数据先分清它是什么、能干什么、不能干什么做水文、规划、灾害评估或者国土空间分析的从业者电脑里缺不了一套分得清高低起伏的数字高程模型DEM。这份洞庭湖流域DEM数据就是把整个洞庭湖水系的地表高程做成了一张张栅格底图你拿到手的是打包好的原始数据不是成品图。它能帮你做流域提取、淹没分析、坡度坡向计算、视域分析也能当三维地形的底图。但先提醒一句它解决的是“地形骨架”问题解决不了“土地利用现状”“河道断面实测”这类需要其他数据源的问题。适合谁用GIS相关的研究生、工程师以及做水文模型、洪涝模拟的同行。拿到手先别急着填洼和提取河网得先把数据底细摸清楚这就是本文要带你做的事。2. 数据底细核查分辨率、坐标系与文件头的完整读数2.1 解压后的第一件事用GDAL把文件头读干净常见的DEM下载来源不同命名和格式也不同。我拿到的这份洞庭湖流域数据解压后是标准的GeoTIFF文件单波段、浮点型高程值。一般人习惯用ArcGIS直接拖进去看这没错但我建议先用命令行工具把元数据读干净省得后面坐标系不对再来回折腾。gdalinfo dongting_dem.tif输出里重点看这几个字段Size is 5333, 4333栅格的行列数。Origin (111.00000,30.50000)左上角坐标。Pixel Size (0.00027777778,-0.00027777778)像元尺寸这里约等于1弧秒对应30米分辨率。Coordinate System坐标系描述。Band 1 Block256x256 TypeFloat32波段类型为32位浮点高程值带小数。NoData Value-9999无效值标记。gdalsrsinfo dongting_dem.tif | head -5这条命令用来快速确认投影信息。如果输出里是GCS_WGS_84说明这是地理坐标系单位是度如果是UTM zone 49N之类的说明是投影坐标系单位是米。参数说明Pixel Size的正负号有意义X为正是自西向东Y为负是自北向南这是GeoTIFF的默认排列方式读取软件会自己处理但你手动写程序遍历像素时别把行的遍历方向搞反。TypeFloat32意味着高程值不是整数能保留小数点后的精度做填洼和坡度计算时比整型DEM平滑很多。我一般会把gdalinfo的输出重定向保存成一个文本文件和原始数据放在同一个目录。这样做的好处是万一后来DEM被处理过、坐标系变了或者分辨率被重采样了还能翻出最早的元数据做对照。2.2 坐标系识别与转换地理坐标系和投影坐标系的取舍洞庭湖流域跨纬度范围较大这份数据大概率是WGS84地理坐标系。直接拿地理坐标系做面积、距离量算会出问题因为度不是等距单位。做流域面积统计前必须投影到等积或等距投影。gdalwarp -t_srs EPSG:32649 -r bilinear -tr 30 30 dongting_dem.tif dongting_dem_utm49.tif这条命令把DEM从WGS84重投影到UTM 49NEPSG:32649输出像元尺寸强制为30米×30米重采样方法用双线性。参数说明utm49n的选择依据是洞庭湖约在东经111°~114°之间UTM 49N覆盖东经108°~114°正好覆盖主湖区。-tr 30 30的含义是目标像元尺寸为30米这里的数值单位是米只有投影后才是米。重采样方法里bilinear适合高程连续表面cubic更平滑但会略微抬高山谷、削低谷底near最邻近法适合分类数据不适合DEM。我经常遇到有人用near重采样DEM结果地形台阶感非常明显坡度计算直接失真。做集水区提取和面积统计用投影坐标系做经纬度范围裁剪和影像配准用地理坐标系也没问题。两种坐标系各备一份是最稳妥的做法。坐标系单位适用场景不适用场景WGS84地理坐标系度裁剪、拼接、发布切片面积量算、坡度计算UTM 49N投影坐标系米流域面积、河道长度、坡度跨带拼接需处理带间重叠3. 从DEM到水系水文分析的完整参数流程3.1 填洼、流向与汇流累积的级联关系拿到DEM做水文分析第一步永远是填洼Fill不是直接算流向。真实地形里存在很多局部低洼点这些洼地可能是真实存在的如湖泊、塘坝也可能是DEM噪点造成的假洼地。水文分析默认水流只往低处走遇到洼地就会断流河网到那儿就断了。ArcGIS的Fill工具界面里只有一个参数Z limit但很多人不知道这个参数该怎么填。常见做法是直接把Z limit保持默认空表示全量填平所有洼地。但洞庭湖流域本身就是湖区天然洼地多全量填平会把真实的垸内湖泊和碟形洼地全抹掉后面提取出来的水系会失真。我的做法是分两步先用原始DEM做一次完整的填洼再用原始DEM减填洼结果得到一个“洼地深度”栅格import rasterio import numpy as np with rasterio.open(dongting_dem.tif) as src: dem src.read(1) profile src.profile # 这里假设你已经用ArcGIS填洼工具生成了fill_dem.tif with rasterio.open(fill_dem.tif) as src2: filled src2.read(1) # 计算洼地深度栅格 depression_depth filled - dem # 找出深度超过5米的洼地这些大概率是真实地形 real_depressions depression_depth 5 # 把真实洼地的区域在原DEM上“加回”重新作为有效地形 dem_masked np.where(real_depressions, filled, dem) with rasterio.open(dongting_dem_trimmed.tif, w, **profile) as dst: dst.write(dem_masked, 1)逻辑说明depression_depth filled - dem计算的是每个像元被填起的高度。深度小于5米的洼地当作噪声填掉深度大于5米的区域比如洞庭湖的东洞庭湖、南洞庭湖主湖面那里高程本来就接近甚至低于长江水位保留填平后的湖面高程避免排水方向完全颠倒。这个5米阈值不是通用值你要根据自己流域的地形起伏特征去试。参数说明阈值设得越小河网越连续但会抹掉越多真实洼地阈值设得越大保留的洼地越多但河网越容易断。我一般先用2米、5米、10米三个档各跑一次把得到的河网叠到卫星影像上看哪个跟实际水道最贴合就用哪个。3.2 流向算法与汇流累积D8与多流向的差异流向计算是水文分析的根基。ArcGIS默认用D8单流向算法也就是把水流方向定义为最陡下降方向输出8个方向编码1、2、4、8、16、32、64、128。D8的优点是简单、结果干净缺点是平地和洼地容易产生平行河道而且在扇形冲积平原上水流方向会被强行归类成单一方向。洞庭湖流域恰好有大片平原区D8跑出来的河网在湖区常常呈诡异的平行线。遇到这种情况我一般会用多流向算法如D∞或FD8替代D8。QGIS里的r.watershed支持-m多流向选项或者用SAGA的“Multiple Flow Direction”模块。多流向算法把水量按坡度比例分配到多个下游像元模拟更接近物理真实代价是计算量增大不少。我实际的做法是山区用D8河道约束明显平原湖区用多流向然后在ArcGIS里把两段河网拼接。听起来粗暴但比单一算法在全区跑出来的结果靠谱得多。拼接时注意重叠区域设置一个过渡带比如湖区边界往外扩5公里在过渡带内做加权融合。3.3 河网阈值与出水口为什么阈值50和500差出一条大江汇流累积Flow Accumulation算完后需要设一个阈值来决定哪些网格算河道。阈值的含义是汇流累积量达到该值的像元才被判定为河道。阈值越小河网越密支流越多阈值越大河道越稀疏只剩下主干。我处理洞庭湖流域时阈值从100试到10000。1000左右能提取出湘江、资江、沅江、澧水四条主要干支流的骨架5000以上就只剩下干流了。具体操作import rasterio import numpy as np with rasterio.open(accumulation.tif) as src: acc src.read(1) # 设定阈值生成河道二值栅格 threshold 1000 streams np.where(acc threshold, 1, 0) # 输出河道掩膜 with rasterio.open(streams_1000.tif, w, **src.profile) as dst: dst.write(streams.astype(uint8), 1)逻辑说明acc threshold把汇流累积量超过阈值的像元标记为河道。这个操作的本质是拓扑筛选不是精确的水文模拟。阈值的选择没有标准答案取决于你后续要做什么——做洪水风险图河道密一点更能反映漫滩做水利工程选址河道粗一点更突出主槽位置。参数说明accumulation.tif是ArcGIS填洼、流向计算之后运行Flow Accumulation工具的结果。注意这里的累积量不是水量而是像元个数所以阈值单位是“像元”不是“平方米”。如果一个像元是30米×30米阈值1000对应的汇水面积约90万平方米你可以用这个关系反推合适的阈值。4. 常见问题排查这袋DEM最容易踩的五个坑4.1 数据读取类坑坑一中文路径导致工具反复报错现象ArcGIS或QGIS里DEM能正常显示但运行填洼、坡度工具时提示ERROR 000882: 无法打开栅格数据集或者Python脚本报No such file or directory。原因部分GIS栅格处理工具对中文路径和中文文件名支持不完善尤其是GDAL 2.x早期版本。洞庭湖流域数据如果放在D:\洞庭湖数据\下载\最终版\这种目录下工具链越深越容易触发编码问题。解决在项目根目录建一个纯英文路径比如D:\dt_dem\把zip解压到该目录下文件名也改成dt_dem.tif这类ASCII命名。这步虽然笨但能消除大量随机性报错。坑二zip包二次解压后文件损坏现象解压后打开DEM发现全黑或者只有边框有值gdalinfo显示ColorInterpUndefined或读取像素值全部为0。原因部分下载平台提供的zip包内嵌了多级压缩目录或者用第三方压缩工具解压时破坏了GeoTIFF内部的元数据偏移量。GeoTIFF的结构依赖文件内部的IFD指针解压软件的兼容性问题会导致指针错位。解决重新解压优先用Windows自带的和7-Zip。如果你用的是WinRAR且解压后有问题不要修复直接换7-Zip重新解压。解压完成后立即用gdalinfo验证文件头确认Size、Band字段都正常再进入下一步。坑三NoData值没当回事现象DEM显示正常但做坡度计算时湖区水面的坡度为异常陡峭的45度以上或者沿湖岸出现一圈黑色锯齿。原因DEM的陆地部分有值湖泊水面区域被标记为NoData通常为-9999。坡度工具碰到NoData和有效值相邻时如果NoData没被正确标注会把无效值当作用于计算的有效高程。解决在ArcGIS的坡度工具中确保环境设置里的NoData识别是勾选的在QGIS中右键图层属性将NoData值手动设成-9999。还有一招更稳先用条件函数把NoData区域填充为周围有效值的平均值再做坡度。4.2 计算过程类坑坑四投影没统一就直接提河网现象提取的河网走向完全错误河道出现90度直角转折或者流域边界和现实完全对不上。原因填洼和流向计算要求输入DEM是投影坐标系单位为米。如果你用WGS84地理坐标系的DEM直接跑流向方向的计算会用度当成长度结果毫无物理意义。解决在测试计算之前用gdalwarp把DEM投影到UTM 49N。确认方法用ArcGIS的属性→源查看像元大小如果显示0.00027777778那肯定没投影需要看到30或30.0这类米制单位。坑五填洼填过头河流倒灌现象提取出的河道在上游段出现了逆向流动或者原来应该是分水岭的位置被挖穿两个流域意外连通。原因全量填洼把所有洼地都填平了包括真实存在的垸内湖、蓄洪区、鱼塘群。在湖区这种人工改造很强的区域填洼工具会像推土机一样把这些人类地物全部削平产生错误的汇水路径。解决回到第3.1节的做法用洼地深度阈值过滤后再填洼。或者更直接人工圈定需要保护的洼地范围在填洼前用栅格计算器把那些区域的高程整体加一个常量比如加20米让填洼工具不触碰它们填完后再减去加的常量。5. 出图与导出把DEM做成能用的工程成果5.1 山体阴影与高程分级配色DEM处理完不做可视化等于白算。出图是汇报和验收的关键一步尤其洞庭湖这种大平湖区域配色做不好就是一片灰绿。ArcGIS里生成山体阴影的步骤是Spatial Analyst工具→表面分析→山体阴影。两个关键参数方位角太阳方位默认315度西北方向光照地形起伏最为立体和太阳高度角默认45度。在湖区这种低起伏区域高度角调到30度能让微弱的10米级起伏也显出明暗变化。高程配色的逻辑不是简单从绿到红拉条色带而是要分层设色。洞庭湖湖区海拔大多在20~50米环湖丘陵100~300米西部山区可达1000米以上。把高程分成6~8级湖区段用深浅绿色区分低洼地和台地丘陵段用黄棕色高山用紫灰色。ArcGIS的符号系统里用分类→手动输入断点值不要用默认的等间距分类否则大部分像元会被分到同一个色阶整个图面没有层次。关于透明度的一个参数细节山体阴影作为底层DEM分级着色叠在其上透明度设为35%~50%能同时保留地形的明暗细节和高度颜色。这个手法在洞庭湖这种大片平地尤其出效果不然平地上看不出任何微地貌。5.2 坡度坡向派生与成果参数表坡度坡向是DEM最常用的派生数据。坡度计算的细节ArcGIS的坡度工具里输出测量单位必须选度Degree或百分比Percent。很多初学者不选默认用度但后续用USLE方程算土壤侵蚀时需要用百分比用错单位结果会差一个数量级。坡向输出的说明坡向栅格的值域是0°北到360°东北方向起始循环平地为-1。环洞庭湖区的坡向分析意义不大但东部丘陵和西部山地区域坡向直接决定了光照、植被分布和建房选址。把成果整理成参数表的习惯我从很早开始就保持成果名称数据类型分辨率坐标系用途说明dt_dem_fill.tif浮点栅格30mUTM 49N填洼后DEM做汇流分析dt_flowdir.tif整型栅格30mUTM 49ND8流向栅格编码1/2/4/8/16/32/64/128dt_accum.tif浮点栅格30mUTM 49N汇流累积用于河网阈值提取dt_stream1000.tif二值栅格30mUTM 49N阈值1000的河网掩膜dt_slope_deg.tif浮点栅格30mUTM 49N坡度度用于坡度分级dt_hillshade.tif整型栅格30mUTM 49N山体阴影出图底纹导出时注意所有成果栅格统一用GeoTIFF格式压缩选LZW无损压缩对DEM这种连续曲面压缩率高。不要选JPEG压缩那是有损的反复保存几次地形细节就钝化了。6. 进阶技巧用Python循环批量裁剪和重投影当你的工作范围从一景DEM扩展到整个流域的十几景数据时手动操作每个文件会很痛苦。这里我建议的做法是用Python配合Rasterio写个批处理脚本把解压、投影转换、裁剪一步做完。import rasterio from rasterio.warp import calculate_default_transform, reproject, Resampling from rasterio.mask import mask import glob # 待处理的DEM文件列表 dem_files glob.glob(D:/dt_dem/raw/*.tif) boundary_geojson D:/dt_dem/boundary/dongting_boundary.geojson # 洞庭湖流域边界 for dem in dem_files: out_name dem.replace(raw, clipped).replace(.tif, _utm49.tif) with rasterio.open(dem) as src: # 读取边界矢量 import json import fiona with fiona.open(boundary_geojson) as geojson: bbox geojson.bounds # (minx, miny, maxx, maxy) # 用边界范围裁剪 out_image, out_transform mask(src, [{type: Polygon, coordinates: []}], cropTrue) # 投影转换 dst_crs EPSG:32649 transform, width, height calculate_default_transform( src.crs, dst_crs, src.width, src.height, *bbox) kwargs src.profile.copy() kwargs.update({crs: dst_crs, transform: transform, width: width, height: height}) with rasterio.open(out_name, w, **kwargs) as dst: for band in range(1, src.count 1): reproject( sourcerasterio.band(src, band), destinationrasterio.band(dst, band), src_transformsrc.transform, src_crssrc.crs, dst_transformtransform, dst_crsdst_crs, resamplingResampling.bilinear)参数说明这个脚本里最关键的是calculate_default_transform里的*bbox参数它让重投影后的输出范围对齐流域边界而不是整景覆盖能省掉至少一半的存储空间。Resampling.bilinear的选择理由和第2.2节说的相同DEM作为连续表面不适合用最邻近法。如果你的数据源是地理空间数据云下载的通常下载时就能选投影和范围。但如果你下载的是全球分幅数据比如SRTM的1度×1度分幅这种批处理就是必须的。我说的这种做法可以在30秒左右处理完一景30米分辨率的DEM整套流域十几景数据十分钟内全部完成。另外提一个我常用的小技巧投影转换后的DEM用来做等高线生成时效果比原始地理坐标系好。ArcGIS里的等高线工具会读取栅格的空间参考来计算实际间距用地理坐标系直接生成等高线间距会错得离谱。把投影后的DEM放进等高线工具设置等值线间距为10米或20米出来的矢量等高线才能直接转给CAD用。说起来有点玄学但我后来每次拿到新的DEM数据都会强制自己走一遍这套流程先gdalinfo验明正身再投影到本地UTM带然后才谈得上做水文分析。这个过程救过我很多次因为数据源的坐标系标注和实际内容偶尔对不上只靠肉眼判断根本发现不了。希望这篇文章能帮你在处理洞庭湖流域DEM数据时少走这些弯路。本文还有配套的精品资源点击获取
返回列表