简介:资源包为洞庭湖流域30米分辨率数字高程模型(DEM)数据,面向GIS学习者、科研人员及水利环保从业者,适用于水文分析、洪水模拟、土地利用规划等宏观地形研究场景。包内共6个文件,核心为tif格式栅格高程数据,附带ovr金字塔文件用于快速显示,dbf属性表、tfw世界文件、xml元数据及cpg编码文件为配套支撑,整体压缩包约347.9MB,文件结构清晰完整。已有448人学习/下载。该数据可直接在ArcGIS、QGIS等软件中加载,利用坡度提取、流向分析、流域划分等功能,探究洞庭湖周边地形起伏特征及其对汇水过程的影响,为水资源调度、洪涝风险评估和生态保护提供基础地理数据支撑。同时也可作为区域DEM应用的教学示例,适合需要实际地形数据开展科研项目或课程设计的地学相关人士。 你拿到的这份洞庭湖流域DEM数据.zip,解压后如果只用来渲染山体阴影,那等于把一本水文地形图当成了壁纸。DEM是数字高程模型的栅格数据,每个像元存一个地面高程值,看着只是一张灰度图,但在水利、地理信息这些行当里,它是做流域河网提取、淹没范围初判、库容估算、地质灾害分析的最小底图。这篇笔记把从拿到压缩包到产出可用成果的流程拆开讲:先做数据体检,再处理坐标系和拼接裁剪,然后用pysheds提取河网,最后把常用场景里的翻车点逐个排掉。
2. 从zip到能分析的高程底图:先检查元数据,再处理坐标系与拼接裁剪
2.1 先体检:用gdalinfo看位深、NoData和坐标系统
数据解压之后的第一件事不是丢进QGIS拉伸颜色,而是用gdalinfo看一眼这个栅格到底在用什么坐标系、存的是什么类型的数值。很多初学者会忽略这一步,直接拿着经纬度的WGS84栅格去做坡度和流向计算,结果算出来的河网方向一塌糊涂,还以为填洼参数没调好。
gdalinfo 洞庭湖DEM.tif重点看输出里的三行:Driver、Size、Coordinate System和NoData Value。Driver告诉你这是不是标准GeoTIFF;Coordinate System如果显示GCS_WGS_1984,说明是经纬度地理坐标系,不能直接做水文计算;NoData Value一般是-9999、-32768或0,填洼和河网提取之前必须让工具认得这个值,否则会把空洞当成真实地形。
位深也要看一眼,Type那行会写Byte、Int16或Float32。SRTM、ALOS这类公开DEM通常用Int16存整米高程,处理起来没问题;如果遇到Byte型,说明数据可能已经做过量化分级,直接用会丢失大量地形细节。遇到这种情况我会先怀疑数据源,而不是急着算水文。还有一种麻烦是文件里同时带.tfw、.ovr和.xml,.ovr是金字塔文件,删了还会自动重建,.tfw是坐标参考信息,.xml是元数据,处理时保留它们对结果没有影响,不用管。
2.2 分幅数据怎么处理:gdalbuildvrt快速拼接再统一裁剪
如果解压出来的是几十个分幅tif,各自覆盖洞庭湖流域的一小块,第一步是把它们拼成一个完整栅格。我一般先用gdalbuildvrt生成虚拟栅格,这一步不产生真实数据拷贝,速度极快,确认范围无误后再转成正式tif。VRT是个文本文件,记录各分幅的路径和位置关系,修改起来非常方便,适合先看范围再动手。
gdalbuildvrt -o 洞庭湖.vrt 分幅目录/*.tif gdal_translate -co COMPRESS=DEFLATE -co TILED=YES 洞庭湖.vrt 洞庭湖全流域.tif用gdal_translate把VRT落成实体tif时,我习惯加COMPRESS=DEFLATE和TILED=YES。压缩可以减少磁盘占用,瓦片化让后续读取指定区域时更快,尤其当数据范围覆盖整个洞庭湖流域时,这两项能明显改善处理体验。拼完之后用gdalinfo再确认一次范围,看Origin和Pixel Size是否合理,如果像元尺寸显示成类似0.00027777778这样的度,那说明还是经纬度单位,后续投影的时候要注意重采样。
如果手头已经有流域边界shp,可以顺手做一次裁剪。用gdal.Warp比先拼全图再裁更省事,它能直接在一个命令里完成拼接、重投影和裁剪,后面的内容会细说投影参数。
2.3 为什么洞庭湖流域要投影到UTM而不是直接用经纬度
经纬度坐标系的单位是度,1度经度的实地距离在低纬度地区和高纬度地区完全不一样,坡度、坡向、流向这些依赖距离和角度的计算,必须用平面坐标。洞庭湖流域主体在东经111度到114度之间,北纬28度到30度左右,跨UTM 49N和50N两个分带,主流做法是投影到UTM 49N,也就是EPSG:32649,覆盖湖南大部分地区。用gdal.Warp可以一步完成投影和裁剪。
gdalwarp -t_srs EPSG:32649 -r cubic -tr 30 30 -cutline 流域边界.shp \ -crop_to_cutline -of GTiff \ 洞庭湖全流域.tif 洞庭湖UTM49.tif-tr 30 30把像元重采样成30米见方,和SRTM 30米原始分辨率保持一致;-r cubic用三次卷积重采样,比最邻近法平滑,比双线性法更保边缘,适合地形连续表面。正射影像或分类图用最邻近法才对,DEM必须用插值类算法,否则山峰会被抹成方块。砍到流域边界时,-crop_to_cutline直接按shp形状裁出不规则区域,而不是裁成外接矩形。裁完之后在QGIS里叠加shp检查一次,确认边界没有明显缝隙。
到这里,这份洞庭湖流域DEM数据zip才算变成可用的底图,下一步可以进入水文分析的流程。
3. 提取洞庭湖河网:填洼、流向与累积量一条龙
3.1 填洼的两个关键参数:最大填洼深度与缓冲像元
原始DEM里存在两类影响水流方向的高程异常:一类是真实地形里的封闭洼地,另一类是SRTM在湖泊、水库、陡坡峡谷里的噪点或空洞。填洼很难自动区分这两类,所以大部分工具提供最大填洼深度这个参数,限制每个洼地最多被填多深。洞庭湖湖区地形极其平缓,湖面高程一般在20到40米之间,湖盆里的微起伏很多是传感器噪声,填深限制设得太小,水流会顺着噪点画出混乱的细线;设得太大又会把真实堤岸抹平。
另一个被忽略的参数是缓冲像元。计算流向时,如果DEM边缘正好切在山脊或河岸上,边缘外没有数据,程序会把水硬生生引向边界。给流域边界向外扩一定宽度再计算,之后把边界外的结果裁掉,能避免这种伪河网。扩展宽度我习惯取50到100个像元,也就是1.5到3公里,对30米分辨率DEM来说足够消除边界效应。
from osgeo import gdal import numpy as np src_ds = gdal.Open("洞庭湖UTM49.tif") band = src_ds.GetRasterBand(1) dem_original = band.ReadAsArray() dem_original[dem_original == -9999] = np.nan读入后用NaN代替NoData,是为了让填洼工具自动忽略空洞,而不是把空洞高程当成真实地形填平。这一步在pysheds里也适用,很多奇怪的计算结果都源于NoData没有被正确识别。
3.2 用pysheds把流向和累积量算出来
pysheds是目前处理水文分析里比较顺手的Python库,底层Cython实现,跑30米分辨率的全流域数据不会等太久。它的流程是固定的:填洼、算流向、算累积量、设阈值提河网。每一步都返回一个二维数组,中间结果可以随时存成tif检查。
from pysheds.grid import Grid grid = Grid.from_raster("洞庭湖UTM49.tif") dem = grid.read_raster("洞庭湖UTM49.tif") # 填洼,限制最大填深,避免把堤岸抹平 flooded_dem = grid.fill_depressions(dem, max_depth=10) # D8流向,dirmap是八个方向的位掩码 dirmap = (1, 2, 4, 8, 16, 32, 64, 128) flow_dir = grid.flowdir(flooded_dem, dirmap=dirmap) # 累积量 accum = grid.accumulation(flow_dir, dirmap=dirmap)fill_depressions的max_depth=10意思是单个洼地的填深不超过10米,这是针对湖区平缓地形的经验值,山区流域可以放大到30米甚至更大。flowdir得到的流向数组用1到128的幂值编码D8方向,每个像元指向8邻域中坡度最陡的方向。accumulation返回的是每个像元上游汇入的像元数,不是面积,但乘以像元面积就能得到集水面积。
如果某片区域明显是湖面,但累积量出现细长的高值线,先别急着调阈值,回到填洼参数重新想。湖面在真实世界里是水平面,DEM里的湖面却有微小起伏,水流会沿这些微起伏画出一条假的“湖心河”。处理办法是手工把湖面高程置平,后面第6章讲局部修正时会提到具体做法。
3.3 河网阈值怎么定:从集水面积反推像元数
累积量数组里的值代表汇入该像元的地表水流路径数,要得到河网就得设阈值:累积量大于等于阈值的像元算河道,小于阈值的算坡面产流。阈值太小,河网密得像鱼刺;阈值太大,源头位置会向中下游收缩好多,小河全丢了。准确的做法是结合实际水文站控制断面面积来定,但快速取用可以先按这一条:阈值像元数等于目标源头集水面积除以单个像元面积。
洞庭湖流域内湘江、资水、沅江、澧水四大河流的源头集水面积差异很大,支流源头集水面积可以取3到10平方公里。30米分辨率的像元面积是900平方米,那么3平方公里对应3333个像元。实际操作时取两到三个阈值分别生成河网,叠加到遥感影像或地形阴影图上对照,看哪个跟实际河道贴合最好。
import numpy as np threshold = int(3000000 / (30 * 30)) # 3平方公里对应的像元数 river_mask = accum >= threshold output = np.where(river_mask, 1, 0) # 把河网mask落盘成tif,方便在QGIS里对照验证这里用的river_mask是布尔数组,1代表河道。落到磁盘时建议同时保存一份累积量栅格,这样不用重算就能换阈值再生一张河网。河网生成后如果想划分出各支流的子流域,常见做法是在河网节点的交汇点选择一个出口,用grid.catchment函数生成该出口的上游集水区,再在QGIS里把多个出口的集水区剪出来,效果和商业软件的流域盆地图差不多。
4. 从DEM到工程决策:坡度、库容与淹没模拟的几个实用场景
4.1 坡度和坡向:算地质灾害风险底图的注意事项
坡度是DEM最直接的应用,地表径流速度、土壤侵蚀强度、塌方隐患判别都看它。在GDAL里算坡度前有一个检查项容易被忽略:投影坐标系DEM的像元尺寸在两个方向应该一致或接近,如果x方向和y方向像元尺寸不同,坡度计算软件多数会按像元对角线边长做参考,结果系统性偏大或偏小。洞庭湖流域这类经过gdalwarp重采样后的数据,像元尺寸已经变成规则的30米乘30米,直接用即可。
from osgeo import gdal ds = gdal.Open("洞庭湖UTM49.tif") slope_ds = gdal.DEMProcessing("坡度.tif", ds, "slope", alg="Horn", slope_format="degree")DEMProcessing是GDAL内置的高程衍生算法入口,slope_format="degree"输出度数制坡度,适合工程出图;"percent"输出百分比坡度,适合水土流失方程。算法默认是Horn,它考虑邻域8个像元,比只考虑东、南两个像元的简单算法抗噪声能力强。坡向用同一入口改"aspect"就可以,输出的0到360度方向是下坡方向,和气象上习惯的上风方向容易搞混,出报告时要写清楚。
4.2 水位-库容曲线:用DEM像元面积叠加高程区间
水库或蓄洪区工程里经常需要水位和库容的关系,也就是水位每上升一米能装多少水。这个可以从DEM快速估算:把高程按照每米一个区间分段统计,每一段的面积乘以厚度1米就是该区间的体积增量,累积求和就得到水位库容曲线。相比实测水下地形,DEM库容在低水位段会偏大,因为DEM看不到水面以下的地形,湖底被当成了水面高程。
elevations = dem_original[~np.isnan(dem_original)] water_levels = np.arange(20, 40, 1) volume = [] for wl in water_levels: flooded_area = np.sum(elevations <= wl) * 30 * 30 volume.append(flooded_area)循环里每一次算的是到当前水位的累计面积,不是增量。累计面积乘以米数再累加,得到的是水位容积曲线的前半段。真正的库容曲线应该按水位分段做增量求和,也就是相邻水位的面积平均值乘以1米,这里直接用累计面积做近似,高水位段误差还在可接受范围内。想精细就用scipy.integrate.trapezoid对水位面积曲线积分,效果一样。
4.3 静态淹没范围:低于水位线的像元不是全部真相
防洪评估里最常被问的问题是这个位置淹不淹。很多第一次接触DEM的人会把“低于水位高程的像元全部标成淹没区”当成答案,这在坡地流域还行,在洞庭湖平原湖区会明显高估。原因是湖区有堤防、垸区、公路路基,这些微地形在DEM里可能只有一两米的起伏,但足以挡住水,低于水位线的像元如果和河湖水面没有连通,实际并不会进水。
正确的快速做法是先设定一个出水口或河道水位面,然后在DEM上做连通性分析,只把从河道逐像元漫出去、且高程低于水位的区域算作淹没。用QGIS的r.lake模块或者pysheds的catchment思路都能做,区别在于前者是给一个进水点,后者是给一个出口点。静态淹没图只能用来做初筛,正式方案还是要交给二维水动力模型去算,但DEM给的这张底图能帮你在项目启动阶段快速圈定重点范围。
5. 洞庭湖DEM数据处理的5个高频翻车点:现象与解法
5.1 湖区填洼后出现大片平地河网,像树枝平躺
现象:河网提取结果在湖区出现大量平行短线和环状线,河网密度明显比上游山区高,明明是一片平湖却画出无数条“河”。原因不是填洼没填平,而是湖面DEM有微幅噪声,填洼又把洼地全抹平了,水流方向几乎由浮点误差决定,累积量随机分布后超过阈值的地方就密密麻麻。解决:先用流域边界内的实测水面高程数据把湖面像元统一赋成一个常数值,再做填洼和流向计算;如果拿不到实测高程,就用分区统计把湖面区域的DEM置平。
5.2 裁剪边缘出现放射状直线,河网贴着边界拐弯
现象:河网图上的河道在流域边界处突然变成垂直边界的直线,或者所有河道在边缘汇聚。原因:直接用流域边界shp切DEM后,切掉的区域在流向计算时被当成“挡墙”,落在边缘的像元指向边界外就没有出口,程序把水流强制导向相邻像元,形成沿边界的假河道。解决:裁剪DEM时向外缓冲至少50个像元,算完河网再把边界外的栅格裁掉,或者提前把掩膜范围扩展到整个外接矩形,别让边界贴着流域线。
5.3 坡度值高得离谱,山体像被压扁
现象:同一个山区算出来的坡度大面积超过60度甚至接近90度,和实地认知明显不符。原因:DEM用了未经投影的经纬度坐标系,东西向像元宽度在纬度30度附近只有约96米,南北向是111米,坡度计算器按几何距离算自然失真,越往北偏差越大。解决:回到第2章补做gdalwarp -t_srs EPSG:32649,确认Pixel Size显示的单位是米,再重跑DEMProcessing。这类错通常不是工具坏了,是数据坐标系没到位。
5.4 NoData空洞把湖底变成0米,库容曲线突然跳变
现象:水位-库容曲线在某一区间出现折线跳变,面积突然增加几千平方公里;淹没模拟里某个湖区边缘出现大块蓝色,但地面高程远高于水位。原因:SRTM在湖面区域经常出现回波空洞,NoData值如果没被正确识别成NaN,工具会把NoData当真实的0米高程参与运算,等于凭空挖出一个大坑。解决:读数据后先强制给NoData赋值成NaN,处理完检查空洞是否被插值填补;如果空洞面积太大,用周围高程做克立金插值补洞,别依赖填洼工具顺手填平。
5.5 用默认参数重采样后,湖面面积比实测大一圈
现象:把原始30米DEM重采样成90米后,湖面边界往外扩了半公里,面积统计显著变大。原因:重采样默认用最邻近法时,湖岸像元被保留;用双线性或三次卷积时,水陆边界像元会插值出一个过渡带,低于水位阈值的部分就全被划成水面。解决:重采样加分带处理,只对陆地像元插值,水面像元保持NaN或常数;要保留湖面边界精度就干脆不重采样,直接沿用原始像元的湖陆分类结果。这个坑在面积量算和淹没模拟里影响最隐蔽,因为成果图看起来比实际更“平滑”。
6. 进阶:用手持GPS和高程基准点校验DEM,再做局部修正
处理洞庭湖流域这类地势极平、水网密布的区域时,我最看重的一步是用实测高程点校验DEM。SRTM在陆地上整体精度不错,但湖区周边人为改造强烈,坑塘、堤坝、公路路基频繁变动,而这些恰恰是水文分析里最敏感的区域。收集测区内水准点、手持GPS打点或已有工程勘察点,和DEM像元值一一对比,算平均误差和RMSE,如果偏差超过一米,就要考虑对局部区域做修正。
# 实测点文件格式:经度、纬度、高程 import pandas as pd from osgeo import gdal points = pd.read_csv("实测点.txt", sep=",") ds = gdal.Open("洞庭湖UTM49.tif") for i, row in points.iterrows(): # 把经纬度投影到UTM,再按像元坐标取高程 x, y = transform_to_utm(row["lon"], row["lat"]) result = gdal_ogr_get_elevation(ds, x, y)这只是一段示意逻辑,具体步骤是:实测点先投影到EPSG:32649,gdal读取像元值后与实测高程相减。修正时不要在全局用平移法统一加减,那样会把山体误差也带走;只针对湖面区域和河网沿线做局部置平或插值。比如湖面高程已知为一个统计值,就把湖面像元批量赋成这个值,确保流向计算和淹没模拟不在地形微幅噪声上翻车。
这套“先体检,再投影,后计算,最后验证”的流程,我现在拿到任何DEM数据包都会先走一遍,哪怕只花十分钟查坐标和NoData,也比算完一版错误河网再回头查数据靠谱得多。希望帮到你。
本文还有配套的精品资源,点击获取