简介这套中国河流及流域矢量数据包面向GIS初学者、水利或地理相关专业学生以及需要基础水系底图进行空间分析的技术人员。资源聚焦河流分级与流域划分直接对应日常GIS操作中常见的线状河流和面状流域数据需求可服务于河流网络制图、流域范围界定、洪水淹没风险分析等任务。压缩包共46个文件核心为SHP矢量文件同时包含shx几何索引、dbf属性表、prj坐标投影等配套文件以及sbn/sbx空间索引和jpg预览图整体大小约9.7MB解压后可按图层直接导入ArcGIS、QGIS等主流平台使用。数据覆盖一级至五级河流以及三级流域等多个分层既能支撑水资源分布研究也便于课堂讲解矢量数据组织方式与河流分级体系。目前已有138人学习下载相比从零自行绘制这套数据能明显减少数据预处理时间适合作为水文GIS项目的起步资料。1. “河流及流域.zip”在实战中到底是哪路数据收到一个叫“河流及流域.zip”的压缩包我一般不会急着解压双击先想清楚包里的东西将要拿来干什么。这类命名在项目里太常见了多是外业调查团队、水利设计院或协作方发来的地理数据集合它不是软件而是一套待加工的底料可能是河流线、湖泊面、流域边界和断面位置也可能连带覆盖研究区的数字高程模型。用途决定流程有人要拿它做淹没分析有人要划分子流域做水文计算还有人只是想出一张带河流和流域边界的规范图。这三种诉求对数据质量的要求完全不同但共用同一条预处理路线。这篇文章就是按我自己的工程习惯把“收到包之后怎么做”完整捋一遍先看家底再统一坐标系接着做河网提取和子流域划分最后把最容易翻车的地方按踩坑记录写清楚。2. 拆包先拆“身份”格式、坐标系与文件命名的门道几乎每个刚拿到“河流及流域.zip”的人都会犯同一个错解压之后只认名字。“河流.shp”“流域边界.shp”看起来很清楚但真正决定数据能不能用的是文件后缀、坐标系统和几何类型。名字起得再漂亮坐标系不对、几何类型不匹配每一步都会卡住。2.1 后缀比文件名诚实shapefile 与栅格的识别压缩包解压时我建议先建一个纯英文、不带空格的目录比如D:/hydro_data/river_basin/。中文目录在部分 GIS 组件和 Python 路径处理中会触发编码问题尤其老旧的 GDAL 版本中文路径能直接把读取变成乱码。用 bash 的解压命令加文件类型识别能快速看清家底。mkdir -p /d/hydro_data/river_basin unzip 河流及流域.zip -d /d/hydro_data/river_basin cd /d/hydro_data/river_basin ls -la file *.shp *.tif *.prj 2/dev/null解压后用file命令看二进制头能区分矢量文件shp与栅格文件tif。ls -la的作用是看文件大小和数量很多包不会只有一个 shp而是带着同名但后缀不同的四个兄弟。要注意shapefile 的完整形态是至少三个文件存放几何的.shp、存放属性表的.dbf、存放索引的.shx如果缺少.dbf属性字段会全部丢失只剩一个没有名字的几何壳。.prj文件也不可小看它记录坐标系统但不少第三方提供的数据里根本没有这个文件。2.2 坐标系“黑匣子”没有 .prj 时默认按 WGS84 猜.prj缺失是最常见也最坑的情况。绝大多数桌面软件在缺.prj时会默认按 EPSG:4326 也就是 WGS84 经纬度来渲染但如果这份数据实际是用高斯-克吕格投影做的坐标值八位数你看到的就是一堆“漂到海里去”的线。判断方法很直接看坐标数值的量级。import geopandas as gpd gdf gpd.read_file(/d/hydro_data/river_basin/rivers.shp) print(gdf.geometry.x.min(), gdf.geometry.x.max()) print(gdf.geometry.y.min(), gdf.geometry.y.max()) print(gdf.crs)这段代码读入河流线并打印横纵坐标范围。如果 x 在百万量级、y 在几十万量级几乎可以断定是投影坐标系如果 x 在 70 到 140 之间、y 在 0 到 60 之间则是经纬度坐标。gdf.crs在缺.prj文件时会返回 None此时不要盲目赋值建议通过坐标量级反推投影带号再找甲方确认。投影带号猜错的后果是整条流域边界在空间上错位几百米后续所有分析全部白做。2.3 一个小时内跑完的数据体检清单拿到数据后的第一个小时我一般按下面这张表做体检不做完不进入下一步。表格里的检查项不是走形式每个都对应一个实际会发生的故障。数据项检查内容通过标准河流线几何类型是否 LineString是否存在空几何全部有效线要素连续流域边界几何类型是否面要素是否闭合Polygon 且无边界开口坐标系统有无.prj坐标量级是否与研究区匹配与 DEM 一致或可转换属性表是否有河流名称、等级字段关键字段无全空高程数据栅格有无有效值有无明显负值空洞数值范围真实可信配合表内检查我还会跑下面这段代码来抽查看几何问题。这一步花不了多少时间却能把“到分析阶段才发现的坑”提前到起点。import geopandas as gpd rivers gpd.read_file(/d/hydro_data/river_basin/rivers.shp) basin gpd.read_file(/d/hydro_data/river_basin/basin.shp) print(空几何数量, rivers.geometry.is_empty.sum()) print(几何类型, rivers.geometry.geom_type.unique()) invalid rivers.geometry.is_valid.sum() print(有效几何, invalid, /, len(rivers))逻辑说明is_empty是判断几何对象是否存在is_valid是判断拓扑是否合法。常见的不合法包括自相交线、未闭合的环。对河流线来说自相交不一定致命但对流域边界来说未闭合面会造成后续裁剪面积偏差。看到不合法结果后先不要动手改往下看坐标系和裁剪很多问题会在投影转换后自动暴露或自动修复。3. 把“河流及流域”变成一套能用的底图坐标统一与数据清洗体检完成后下一步不是立即做水文分析而是把数据整理成一个干净、稳定、可复用的底图。这一步有三分统一坐标系、按照流域边界裁剪、清洗河网的几何毛病。很多团队跳过这段直接跑模型结果模型报错后回头排查才发现是输入数据本身的问题。3.1 统一坐标系先定义再投影当包里有 DEM 栅格又有河流矢量时两者坐标系大概率不一致这是项目的常态。栅格常用 WGS84 或国家 2000 地理坐标系矢量常被转成高斯投影或 UTM 投影。统一坐标系的最低原则是让矢量和栅格最终落在同一个投影坐标系上单位必须是米。面积、坡度、累积流量统统依赖米的单位用经纬度算面积是错的误差随纬度变大。import geopandas as gpd import rasterio rivers gpd.read_file(/d/hydro_data/river_basin/rivers.shp) dem rasterio.open(/d/hydro_data/river_basin/dem.tif) print(矢量原始坐标系, rivers.crs) print(栅格原始坐标系, dem.crs) # 如果矢量没有坐标系先按坐标量级手动指定再转换 if rivers.crs is None: rivers rivers.set_crs(EPSG:4326) target EPSG:32650 # 以研究区所在的 UTM 投影带为示例 rivers_proj rivers.to_crs(target) dem_proj_path /d/hydro_data/river_basin/dem_utm.tif # 用 rasterio 把 DEM 重投影到同一条坐标系 with rasterio.open(dem) as src: profile src.profile.copy() profile.update(epsg32650, driverGTiff) transformed src.read( out_shape(src.count, src.height // 2, src.width // 2), resamplingrasterio.enums.Resampling.bilinear, ) with rasterio.open(dem_proj_path, w, **profile) as dst: dst.write(transformed)参数说明EPSG:32650是 UTM 50 带的代码覆盖东经 114 度到 120 度区域你要根据自己的研究区经度换带号换算规则是经度除以 6 取整加 31。set_crs只在缺少坐标系时用它不改变坐标数值只给数据打上标签to_crs才真正做坐标转换。栅格重投影时我把out_shape缩小一半这是为了测试提速正式计算时不要缩直接用原始分辨率。定义和投影的区别容易混淆一句话讲清set_crs是告诉软件“我这个数是什么坐标系”to_crs是让软件“帮我把这个数换算到另一个坐标系”。数据源没有.prj文件时先定义再投影顺序反了会产生把经纬度数值当米去做投影的灾难。3.2 用流域边界把范围裁干净裁剪不只是为了减小数据量更是为了让后续的河网提取和子流域划分不被流域外的地形干扰。比如河流数据里常有干流上游伸到流域边界之外或者流域边界本身是设计院手工勾绘的和 DEM 的边缘不对齐。正确的裁剪姿势是先看边界与 DEM 是否在空间上咬合再执行裁剪。import geopandas as gpd from shapely.geometry import mapping basin gpd.read_file(/d/hydro_data/river_basin/basin.shp) rivers gpd.read_file(/d/hydro_data/river_basin/rivers.shp) # 先修复几何buffer 0 能处理大部分拓扑错误 basin_clean basin.geometry.buffer(0).unary_union # 按流域裁剪河流 rivers_clipped gpd.clip(rivers, basin_clean) rivers_clipped.to_file(/d/hydro_data/river_basin/rivers_clipped.shp) print(裁剪前河流数量, len(rivers)) print(裁剪后河流数量, len(rivers_clipped))这段代码里buffer(0)很关键。它不改变面形状但会重建几何对象把面内自相交的环和微小的拓扑裂隙修掉。unary_union的作用是把多部件面合并成一个整体保证后续clip不会因为多个面互相重叠而造成重复裁剪。裁剪后输出文件可仔细检查河流数量如果断崖式减少很可能是流域边界和河流本身在空间上就没有真正贴合需要回到坐标系检查那一环。3.3 河网数据清洗悬空边和“回头河”第三方给的河网数据里最常见的问题是“挂断头线”和“方向混乱”。挂断头线指河网在分叉处没有完全连通看起来断了一小截方向混乱则是指线要素的起点终点没有统一规则有的从上游到下游有的从下游到上游。水文分析时这两类问题都会污染后续的处理尤其当你依赖矢量河网去确定出口点位置时。import geopandas as gpd rivers gpd.read_file(/d/hydro_data/river_basin/rivers_clipped.shp) # 检查悬挂点线端点不在任何其他线上 from shapely.geometry import Point, LineString coords [] for line in rivers.geometry: if isinstance(line, LineString): coords.append(Point(line.coords[0])) coords.append(Point(line.coords[-1])) # 统计端点重复次数端点在整条河网中只出现一次的即为悬挂点 from collections import Counter count Counter([(p.x, p.y) for p in coords]) dangling [p for p in count if count[p] 1] print(疑似悬挂端点数量, len(dangling))逻辑说明Counter统计所有线要素的端点坐标河流交汇处会被两条以上的线共享因此正常交汇点出现次数大于 1只有出现次数为 1 的点才是断头。这个脚本只是初筛不能直接删除这些线残缺河网可能是数据采集缺漏也可能是一侧是支流源头。处理方式通常是保留主河道对支流断裂处做人工复核而不是批量删数据。“回头河”则是线方向不统一的问题判断方法和端点统计类似但要看线首尾坐标是否落在流域上游。我自己的习惯是矢量河网只用来做位置对照不作为水流方向依据真正的水流方向全部交给 DEM 计算这也是后面第 4 章能跑通的前提。4. 从 DEM 到河网与子流域一套最小可复现的流域分析流程有了干净的底图接下来是整篇最核心的部分用 DEM 提取河网、划分子流域。这套流程是流域分析的“标准动作”不管包里的矢量数据多全最终能够算流路、算汇水面积的必须是一套连续、无裂缝的栅格计算链填洼、算流向、算累积流量、提河网、定出口、生成子流域。4.1 填洼为什么是第一步阈值不是越大越好DEM 数据里存在大量“假洼地”成因包括原始数据噪声、插值算法残留和地形过度平滑。如果不填洼水流会在洼地里打转累积流量结果里会出现一道道不合理的黑色斑块河流在平地上断成好几截。填洼的常用工具是 pysheds它读入 DEM 后直接对洼地像素进行抬升处理让地形表面处处向流域边界倾斜。import numpy as np import rasterio from pysheds.grid import Grid grid Grid.from_raster(/d/hydro_data/river_basin/dem_utm.tif) dem grid.read_raster(/d/hydro_data/river_basin/dem_utm.tif) # 填洼pysheds 默认只处理真实洼地不强行抹平平地 filled_dem grid.fill_depressions(dem, apply_small_depression_filterFalse) grid.write_raster( /d/hydro_data/river_basin/dem_filled.tif, filled_dem, dtypefloat32, overwriteTrue, )逻辑说明fill_depressions的作用是找到所有无法排水的低洼像元把它们提升到周边能够排水的最低高度。apply_small_depression_filter控制是否先对极小的洼地做滤波我一般设为 False原因在踩坑章里细说。禁用滤波器时填洼只动真正不连通的洼地而不是把整片缓坡推平。填洼的参数不是越大越好。很多教程建议把apply_small_depression_filter打开并调大滤波窗口这会把天然洼地、喀斯特漏斗、人工坑塘全部抹掉得到的“平坦面”看似平滑实际上让后续坡度失真。正确思路是先用原始分辨率算一遍看填洼后生成的流向在关键河道处是否与真实河流吻合不吻合再考虑微调。4.2 流向算法选 D8 还是 D-inf两个标准和一个折中填洼完成后计算流向。流向算法有很多工程里用得最多的两个是 D8 和 D-inf。D8 把水流分配给 8 邻域中坡度最陡的一个方向逻辑简单、结果稳定但在地形平缓区域会形成大量平行直线河道。D-inf 将水流按坡度比例分配给相邻两个方向更适合模拟真实地表漫流但子流域边界会更破碎。from pysheds.grid import Grid grid Grid.from_raster(/d/hydro_data/river_basin/dem_filled.tif) filled_dem grid.read_raster(/d/hydro_data/river_basin/dem_filled.tif) # 计算流向method 可切换 d8 或 dinf flow_dir grid.flowdir(filled_dem, methodd8) # 根据流向计算累积流量d8 用 bool 型dinf 用对应权重 accum grid.accumulation(flow_dir, methodd8) grid.write_raster( /d/hydro_data/river_basin/accumulation.tif, accum, dtypefloat32, overwriteTrue, )参数说明method控制算法类型只有两种选择d8 与 dinf通俗写法是d8和dinf。选择标准我总结成两条第一研究区是山地丘陵、河道单一稳定选 D8后续划出的子流域边界干净第二研究区是平原、水网密集、人工沟渠纵横选 D-inf 更接近真实汇流。折中方案是先用 D8 跑通整条流程再用 D-inf 做一次敏感性对比如果两种算法下河网和子流域面积差异小于 10%说明地形控制强结果可信差异过大则需要检查 DEM 分辨率是否够细。4.3 河网提取累积流量阈值的三个参考系累积流量图算完后需要设置一个阈值来确定哪些栅格算河道、哪些不算。阈值是汇流网格数高于阈值的网格才能成为河网的一部分。这个参数没有标准答案属于典型的“玄学”参数但我总结出三个可落地的参考系按集水面积反推、按 Strahler 河级约束、按河网密度试算。import numpy as np from pysheds.grid import Grid grid Grid.from_raster(/d/hydro_data/river_basin/dem_filled.tif) accum grid.read_raster(/d/hydro_data/river_basin/accumulation.tif) # 第一种参考按最小汇水面积反推阈值 resolution grid.cell_size # 栅格分辨率单位米 min_area_km2 5 # 最小河流对应汇水面积单位平方公里 threshold (min_area_km2 * 1e6) / (resolution ** 2) print(阈值网格数, threshold) # 第二种参考按干流分叉位置人工选定截取河道 mask accum threshold grid.write_raster( /d/hydro_data/river_basin/river_mask.tif, mask.astype(int8), dtypeint8, overwriteTrue, )参数说明min_area_km2是你要保留的最小河流对应的上游汇水面积山区取 1 到 5 平方公里合适平原区要取 10 平方公里以上否则会把田间排水沟全部识别成河道。第二个参考系是用 Strahler 河级卡河道即提取河网后只保留大于等于二级的河段此法能有效去掉单支短河道。第三种做法是“试算对比”分别用阈值的零点五倍、一倍、二倍提取河网叠到卫星影像上看密度找到与影像水系最接近的那个值。阈值偏低河网密集成刺猬阈值偏高只剩干流宁可先偏低再人工删也不要一开始就高到丢主河。4.4 划分子流域从出口点到集水区河网确定后子流域划分需要先给定出口点。出口点是水流离开研究区的位置也就是流域边界上的最低点通常选在主河道与流域边界交点处。如果没有实测出口坐标一种办法是计算累积流量最大值的位置作为自然出口点。from pysheds.grid import Grid import geopandas as gpd from shapely.geometry import Point grid Grid.from_raster(/d/hydro_data/river_basin/dem_filled.tif) accum grid.read_raster(/d/hydro_data/river_basin/accumulation.tif) # 找累积流量最大值点作为出口需要转回地理坐标 max_accum_idx np.unravel_index(np.argmax(accum, axisNone), accum.shape) x, y grid.affine * (max_accum_idx[1], max_accum_idx[0]) outlet gpd.GeoDataFrame( {id: [0]}, geometry[Point(x, y)], crsEPSG:32650, ) # 生成集水区即流域出口对应的上游汇水范围 catch grid.catchment(flow_dir, xoutlet.geometry.x[0], youtlet.geometry.y[0]) grid.write_raster( /d/hydro_data/river_basin/catchment.tif, catch, dtypeint8, overwriteTrue, ) # 将集水区栅格直接转为矢量面 catch_poly grid.polygonize(catch) gpd.GeoDataFrame.from_features(catch_poly).to_file( /d/hydro_data/river_basin/catchment.shp )逻辑说明np.argmax是快速定位最大累积流量格点位置配合grid.affine将数组行列号转成投影坐标这一步常用但最容易忘记很多人拿着行列号就去打点了。grid.catchment的实质是从出口点开始沿流向反查上游所有能流入该点的像元它生成的栅格是一个大的 0-1 掩膜。polygonize将掩膜转成矢量面但输出的多边形容易出现锯齿状边界建议在 GIS 里做一次平滑或通用化后再出图。当研究区内有多条支流时更精细的做法是给每个河网栅格交点都打上出口点逐个生成子流域。实际操作中我通常用“先划分、后合并”的策略自动划分结果按面积排序面积小于设定值的小斑块并给相邻大斑块避免子流域数量过多参考文献和汇报材料很难对几十个碎片化的子流域一一描述。5. 五个高频翻车点河流与流域处理的避坑清单流程跑通只是起点工程里真正耗时间的是排查日志里一长串莫名其妙的错误。以下五条是我在多次项目中攒下的血泪经验每条都按照“现象、原因、解决”来展开对照你自己的情况排查能省下半天到一天时间。5.1 面积“差一位”算法没变单位变了现象算出的子流域面积数值诡异要么全部只有零点几要么全部上千跟流域实际范围对不上。原因栅格数据本身是经纬度坐标系矢量已经转成投影坐标系但栅格被直接读入分析流程分辨率此时是度而不是米面积公式中混入了度与米的换算错误。解决进入填洼前先检查栅格的 CRS 是否与矢量一致直接打印src.crs与grid.cell_size。如果cell_size数值在 0.0001 这个量级说明分辨率单位是度需要先栅格重投影再分析不要图省事在分析结果上乘一个“大概的系数”。5.2 填洼阈值开太大河道被抹成平地现象河网提取结果只剩干流一条粗线所有支流全部消失累积流量图上一片均匀色。原因在fill_depressions时打开了apply_small_depression_filter且滤波窗口设置过大把真实的小型河谷都当成洼地填平了。解决回退到不滤波只填 True Depression如果研究中确实需要滤掉微地形滤波窗口从 3 开始逐渐增加每加一次都要对比提取的河网与影像上的真实支流。记住填洼是让“积水能流出去”不是让“地形变平坦”。5.3 河流在流域边界处凭空断开现象提取出的河网在流域边界附近戛然而止下游几百米又能看到河道中间被挖掉一段。原因裁剪时用了完整的流域边界直接裁 DEM边界处的栅格因为部分落入边界外而被设为 NoData流向计算时这些区域不参与汇流河网自然断开。解决裁剪 DEM 前先对流域边界做缓冲区将边界向外扩出一到两个栅格分辨率分析完成后在出图时再叠回精确边界。这个缓冲区也顺带解决了边界处像素因四舍五入产生的锯齿问题。5.4 “回头河”叠加对比时让人怀疑人生现象将提取的河网与原始矢量河流叠在一张图上发现部分河段走向相反支流像是从下游流向上游。原因原始矢量数据中的线方向在实际生产过程中没有按统一规范采集有的从河口往源头顶点绘制有的完全相反。解决方向上以 DEM 计算为准矢量河网只用于抽稀和名称标注。如果你必须要修正矢量方向可按线首尾点的高程判断高程高的一端应作为上游起点但这一方法在平原区会失效此时不要强改直接隐藏方向符号出图时用河流名称字段区分即可。5.5 大面积流域跑不动内存被一张全景图压垮现象处理整条大江全流域时DEM 文件不足 1GB但代码运行到accumulation那步突然卡死内存占用飙升到 90% 以上。原因基于网格的算法在内部会维护多份中间数组dem、filled_dem、flow_dir、accum同时驻留内存数据量乘以倍率后远超实际文件大小。解决把大流域按支流流域范围切成几块分别计算导出结果后再拼接。切块时要保证相邻块之间有至少一个流向长度的重叠否则拼接处会出现一条人为河界。如果不想做切块可考虑将栅格降采样一半分辨率跑一次概算验证流程正确后再用全分辨率跑最终成果。6. 用“后悔药”心态做验证不只出图还出可对比的产物流程全部跑完、子流域边界也画出来了很多人的习惯是导出一张图就收工。我的习惯不一样任何一次批量处理前我会先把中间产物全部留档并给每个产物打上参数标签。理由是分析参数反复调了很久最后确认用的参数组合可能与最初版本相差很大没有留档就没有后悔药可吃。我会为每个处理阶段保留三样东西参数记录表、结果为 tif 栅格或 shp 矢量、以及一张包含结果和影像底图的对比图。参数记录表只需要六行写清楚阈值、算法、分辨率、填洼设置、坐标系、版本日期。对比图用于人工判断这一步无法自动化需要盯着看十分钟确认河道没有断头、支流没有异常密集或稀疏、边界没有跑出流域界。一个简单但有效的验证是用提取出的河网与原包里的矢量河流做叠置率计算。将矢量河流栅格化与提取河网做相交分析统计提取河网中落在原始河网缓冲区内的长度比例。比例高于百分之八十说明提取结果可靠低于百分之六十则说明阈值或填洼设置需要调整。此验证办法在人工沟渠密集的区域会失效因为提取结果更加精确而原始矢量漏绘较多此时应人工抽查几个子流域与影像的匹配度。最后补一条经验不要把数据直接覆盖在原文件名上。我吃过一次亏反复调参后顺手用原文件名保存了结果第二天要回退发现原始解压数据已经被覆盖只能重新解压重跑一遍。现在我的输出目录按20250125_basin_area5km这种命名日期加关键参数任何中间版本都可追溯。尽管理论上所有问题都可以重新解决但一次重跑的成本往往远超你预测的时间。希望帮到你。本文还有配套的精品资源点击获取