简介面向GIS学习、流域水文分析与生态环境保护研究黄河流域空间数据集将流域面要素、一级区边界线、DEM栅格与矢量河网整合在同一压缩包内适合在ArcGIS中开展叠加分析、网络分析、流域划分与地形可视化也可用于土地利用变化、洪水监测和水资源管理的课程练习与专题制图。包内共40个文件约2.15MB以shp、shx、dbf、prj、sbn、sbx等Shapefile矢量及索引投影文件为主另有adf、dat、nit等DEM栅格文件、xml元数据与log日志坐标系统和空间索引信息较完整便于直接加载与检索。目前已有267人学习。使用者可据此搭建黄河流域基础地理数据库结合河网与高程数据模拟水文循环、评估土地利用变化对河流流量的影响并利用边界与面要素辅助流域规划、洪水预警和生态变化分析适合作为中高级GIS练习与科研前期底图。1. 面要素、边界线、DEM、矢量河网这套 shp 资源到底能省掉哪几步一套按 shp 文件组织、准备在 ArcGIS 里使用的大型流域面要素、边界线、DEM、矢量河网资源最怕的不是不会打开而是打开后直接开始裁切。面要素决定统计单元边界线决定裁切范围DEM 决定高程与坡度语境矢量河网决定水系连通关系四者坐标系、范围、字段任一处对不上十分钟的操作能换来两小时返工。这份资源适合做流域分区统计、河网密度核算、DEM 水文派生、专题制图的人也适合刚接触 ArcGIS 空间分析、需要一套完整分层数据练手的人。它不是单个图层而是一组需要先体检、再使用的空间底座。2. 先看坐标系与数据分层ArcGIS 打开前要核对的四件事2.1 面要素与边界线的几何角色谁负责统计谁负责裁切面要素通常是流域分区、子流域或行政统计单元字段里往往带分区编码、名称、面积这类属性适合做汇总和制图。边界线则是外边界或分段界线几何上是线不是面直接拿它去裁面容易遇到“裁剪要素不是面”的报错。我一般先把边界线转成面再检查闭合性和悬挂点最后才用面去裁面要素。这样做的原因是 ArcGIS 的裁剪工具对裁剪几何有要求线要素即使看起来围成一圈也不等于闭合多边形。边界线转面后还要看属性是否保留如果边界线带分段编码转面时容易把多条线合成一个面后续统计会丢分段信息。常见做法是保留一份原始边界线做核对另一份转面只用于裁切不要在同一条线上反复改。2.2 DEM 的像元、Z 单位与 NoData别把高程当装饰DEM 是栅格不是普通图片。打开后先看像元大小、Z 单位、NoData 值和统计范围比看颜色更重要。像元大小决定河网提取的粒度Z 单位决定高程是米还是其他单位NoData 决定边缘会不会被当成 0 参与填洼。很多“填洼后出现大片平地”的问题根子不是填洼算法而是 NoData 被错误处理。ArcGIS 里可以在栅格属性里看像元大小和 NoData也可以用脚本批量读。DEM 与矢量河网联动时还要看两者范围是否一致DEM 比边界线小一圈河网在边缘就会断头DEM 比边界线大很多填洼和流量累积会多算外围汇水。我的习惯是先按边界线外接矩形裁一版工作 DEM保留原始 DEM 不动所有派生都放在工作目录里。2.3 矢量河网的拓扑与属性流向、等级、名称字段怎么读矢量河网是线要素重点看三件事线方向、等级字段、名称字段。线方向影响后续网络分析和符号化箭头等级字段影响分级设色名称字段影响标注和查错。很多 shp 河网的线方向是乱的直接做流向箭头会一半朝上一半朝下这时不要手工翻线先用要素折点转点或空间查询判断上下游关系再决定是否统一方向。河网属性里如果有“长度”字段要确认它是地理坐标系下的度还是投影坐标系下的米度值不能直接拿来做河网密度。拓扑方面重点看悬挂点支流没有接到干流、河道在边界处断头、两条线重复叠加都会让后面的密度统计偏高或偏低。先把这些几何问题定位清楚再谈符号化。2.4 用 ArcGIS 做一次只读体检目录、投影、范围、字段拿到资源后先跑只读体检不添加字段、不裁剪、不投影只读元数据。下面这段脚本会递归列出目录下的要素类和栅格打印几何类型、坐标系、范围和字段。路径换成你自己的资源目录即可。import arcpy import os # 资源根目录换成你自己的路径脚本只读不修改任何数据 root rF:\模拟项目X\流域数据 # 递归遍历要素类shp 也会被 Walk 识别 for dirpath, dirnames, filenames in arcpy.da.Walk(root, datatypeFeatureClass): for fc in filenames: path os.path.join(dirpath, fc) desc arcpy.Describe(path) sr desc.spatialReference print(要素类:, path) print( 几何类型:, desc.shapeType) print( 坐标系:, sr.name if sr and sr.name else 未定义) print( 范围:, desc.extent.XMin, desc.extent.YMin, desc.extent.XMax, desc.extent.YMax) # 排除几何和 OID 字段重点看属性字段名与类型 fields [(f.name, f.type, f.length) for f in arcpy.ListFields(path) if f.type not in (Geometry, OID)] print( 字段:, fields) # 单独列栅格重点看像元大小、坐标系和范围 for dirpath, dirnames, filenames in arcpy.da.Walk(root, datatypeRasterDataset): for ras in filenames: path os.path.join(dirpath, ras) desc arcpy.Describe(path) sr desc.spatialReference print(栅格:, path) print( 像元大小:, desc.meanCellWidth, desc.meanCellHeight) print( 坐标系:, sr.name if sr and sr.name else 未定义) print( 范围:, desc.extent.XMin, desc.extent.YMin, desc.extent.XMax, desc.extent.YMax)这段脚本的逻辑是先走目录再读 Describe 元数据不打开编辑会话。arcpy.da.Walk的datatype参数控制只列要素类还是栅格避免把无关文件混进来。desc.spatialReference.name为空时说明数据没有.prj或坐标系未定义不要急着用“定义投影”硬套先确认原始来源。字段列表里重点看中文名、长度和类型后面导出 shp 时字段名超过 10 个字符会被截断这一步提前暴露问题。注意只读体检阶段不要运行“修复几何”或“投影”覆盖原文件先复制一份工作副本所有修改落在副本上。3. 面要素与边界线落地裁剪、拓扑检查、面积统计的完整链路3.1 用边界线裁面要素Clip 与 Pairwise Clip 的差别边界线不能直接当裁剪面用常见做法是先把闭合边界线转成面再用面裁面。ArcGIS 里有 Clip 和 Pairwise Clip 两类工具前者老牌后者对大数据和并行处理更友好参数语义也略有差别。我的做法是边界线先FeatureToPolygon转面检查转出来的面数量是否与预期一致再用PairwiseClip裁面要素。下面这段脚本先转面再裁切输出到工作目录不覆盖原始面要素。import arcpy arcpy.env.overwriteOutput True # 原始边界线与待裁面要素 boundary_line rF:\模拟项目X\流域数据\边界线.shp in_fc rF:\模拟项目X\流域数据\面要素.shp # 工作目录与输出 work_gdb rF:\模拟项目X\工作.gdb boundary_poly rF:\模拟项目X\工作.gdb\边界线_转面 out_fc rF:\模拟项目X\工作.gdb\面要素_裁剪 # 线转面容差按数据坐标单位给投影坐标下常用 0.001 米级 arcpy.management.FeatureToPolygon( boundary_line, boundary_poly, 0.001, ATTRIBUTES, NO_LABEL ) # 用转出来的面裁面要素 arcpy.analysis.PairwiseClip(in_fc, boundary_poly, out_fc) print(裁剪完成:, out_fc)FeatureToPolygon的第三个参数是聚类容差单位跟输入坐标系一致投影坐标下给 0.001 表示 1 毫米级容差用来合并几乎重合的折点给太大可能把该分开的线并到一起。ATTRIBUTES表示把边界线属性带到面上方便后面按分段统计。PairwiseClip的裁剪要素必须是面输出会继承输入面要素的字段结构。裁剪后先看要素数量有没有异常减少再看边缘有没有被裁出细碎多边形细碎面往往是边界线不闭合或容差过小导致的。3.2 拓扑检查悬挂点、自相交、缝隙怎么定位面要素裁剪完不要直接算面积先建拓扑。面拓扑常用规则是“不能有重叠”“不能有缝隙”“不能自相交”边界线拓扑常用“不能有悬挂点”“不能有伪节点”。拓扑不是玄学它就是把几何关系显式化让错误要素落在一个可导出的错误表里。下面脚本在文件地理数据库里建一个拓扑把面要素加进去并加三条规则最后验证并导出错误。import arcpy arcpy.env.overwriteOutput True gdb rF:\模拟项目X\检查.gdb fc rF:\模拟项目X\工作.gdb\面要素_裁剪 topo rF:\模拟项目X\检查.gdb\面拓扑 err_fc rF:\模拟项目X\检查.gdb\面拓扑错误 # 建库、建拓扑、加要素类 arcpy.management.CreateFileGDB(rF:\模拟项目X, 检查.gdb) arcpy.management.CreateTopology(gdb, 面拓扑) arcpy.management.AddFeatureClassToTopology( topo, fc, 1, 1, 面要素_裁剪 ) # 三条面拓扑规则不重叠、无缝隙、不自相交 arcpy.management.AddRuleToTopology(topo, Must Not Overlap (Area), fc) arcpy.management.AddRuleToTopology(topo, Must Not Have Gaps (Area), fc) arcpy.management.AddRuleToTopology(topo, Must Not Self-Intersect (Area), fc) # 验证并导出错误要素 arcpy.management.ValidateTopology(topo) arcpy.management.ExportTopologyErrors(topo, gdb, 面拓扑错误) print(拓扑错误已导出:, err_fc)AddFeatureClassToTopology的 XY 容差和 Z 容差按数据单位给脚本里用 1 和 1 是占位写法实际要按你的坐标单位调整投影坐标下 1 米容差可能过大地理坐标下 1 度更离谱。Must Not Have Gaps会把相邻面之间的缝隙标出来但边界处的开口不一定是错误要结合边界线判断。ExportTopologyErrors导出的错误要素保留拓扑规则和错误类型后面可以按RuleType筛出重叠或缝隙逐条处理。拓扑错误不要一键全删先看错误位置是否落在真实分区边界上。3.3 面积与周长字段投影坐标系下重算不在地理坐标系里算地理坐标系下算出来的面积单位是“平方度”没有业务意义。面积统计前必须确认数据在投影坐标系下再用CalculateGeometryAttributes重算面积和周长。下面脚本先加字段再按平方米重算字段名用短名避免 shp 截断。import arcpy arcpy.env.overwriteOutput True fc rF:\模拟项目X\工作.gdb\面要素_裁剪 # 加短字段名方便后续导出 shp for name, ftype in [(AREA_M2, DOUBLE), (PERIM_M, DOUBLE)]: if name not in [f.name for f in arcpy.ListFields(fc)]: arcpy.management.AddField(fc, name, ftype) # 面积按平方米周长按米要求当前要素类已处于投影坐标系 arcpy.management.CalculateGeometryAttributes( fc, [[AREA_M2, AREA], [PERIM_M, PERIMETER]], area_unitSQUARE_METERS, length_unitMETERS ) print(面积与周长已重算)CalculateGeometryAttributes的area_unit和length_unit只在当前坐标系支持对应单位时生效如果要素类还是地理坐标系脚本会给出警告或结果不可信。字段名用AREA_M2、PERIM_M是为了避开 shp 的 10 字符限制同时保留单位语义。重算后抽查最大面和最小面跟边界线范围做一次量级核对如果面积总和远大于边界线外接矩形面积说明面有重叠或重复记录。3.4 导出成果与字段别名Shapefile 字段名限制的绕法Shapefile 的字段名最多 10 个字符字段别名不保存DBF 编码还可能让中文属性乱码。我的习惯是中间过程全放文件地理数据库交付时再按需导出 shp如果对方必须用 shp字段名提前改成短名中文属性用 UTF-8 或按对方环境约定编码。下面脚本检查导出后的字段长度和中文可读性避免交付后才发现字段被截断。import arcpy fc rF:\模拟项目X\交付\面要素_裁剪.shp for f in arcpy.ListFields(fc): if f.type in (Geometry, OID): continue print(字段:, f.name, | 长度:, len(f.name), | 类型:, f.type) if len(f.name) 10: print( 警告字段名超过 10 字符shp 可能截断)这段脚本只读字段表不修改数据。len(f.name)判断的是实际字段名长度不是别名长度shp 不保存别名所以不要指望在 shp 里保留长中文别名。字段名超过 10 字符时导出前先在 GDB 里改成短名再复制到 shp。中文属性乱码则要看 DBF 代码页必要时用“表转表”指定代码页或者直接交付文件地理数据库减少编码环节。4. DEM 与矢量河网联动填洼、流向、流量累积、河网栅格转矢量4.1 DEM 预处理填洼、投影、裁剪的顺序DEM 预处理的顺序会影响结果。我的习惯是先把 DEM 投影到与面要素一致的投影坐标系再按边界线外接矩形裁一版工作 DEM最后填洼。先投影再填洼是为了让流向和流量累积的像元面积有明确单位先裁剪再填洼是为了避免外围汇水把边界内的洼地填平。填洼不是越狠越好填洼量过大可能抹掉真实洼地填洼量过小又会让流向在洼地内部打转。下面脚本按“投影—裁剪—填洼”的顺序走输出全部放工作目录。import arcpy from arcpy.sa import Fill arcpy.env.overwriteOutput True dem rF:\模拟项目X\流域数据\dem.tif boundary_poly rF:\模拟项目X\工作.gdb\边界线_转面 dem_proj rF:\模拟项目X\工作.gdb\dem_proj dem_clip rF:\模拟项目X\工作.gdb\dem_clip dem_fill rF:\模拟项目X\工作.gdb\dem_fill # 以边界线转面的坐标系为目标投影 target_sr arcpy.Describe(boundary_poly).spatialReference # 投影栅格像元大小沿用输入避免二次重采样 arcpy.management.ProjectRaster(dem, dem_proj, target_sr) # 按边界线转面裁剪保持栅格范围与矢量范围一致 arcpy.management.Clip(dem_proj, , dem_clip, boundary_poly, NONE, , NO_MAINTAIN_EXTENT) # 填洼输出浮点栅格 Fill(dem_clip).save(dem_fill) print(DEM 预处理完成:, dem_fill)ProjectRaster的目标坐标系直接取边界线转面的坐标系避免手工选错。Clip的裁剪几何用面NO_MAINTAIN_EXTENT表示输出范围贴合裁剪面不保留输入栅格的完整外接矩形。Fill只处理洼地不改变像元大小如果填洼后出现大片同值平地回头检查 NoData 是否被当成 0或者裁剪边缘是否把外围汇水切断了。4.2 流向与流量累积阈值不是拍脑袋流向用 D8 还是多流向取决于后续河网提取的精度要求。D8 简单、稳定适合矢量河网复核多流向更平滑但参数更多。流量累积的阈值不是固定值常见做法是用“最小河网汇水面积 ÷ 像元面积”换算成像元数再拿几组阈值做对比。下面脚本先算流向和流量累积再用阈值提取栅格河网。import arcpy from arcpy.sa import FlowDirection, FlowAccumulation, Con arcpy.env.overwriteOutput True arcpy.CheckOutExtension(Spatial) dem_fill rF:\模拟项目X\工作.gdb\dem_fill flow_dir rF:\模拟项目X\工作.gdb\flow_dir flow_acc rF:\模拟项目X\工作.gdb\flow_acc stream_ras rF:\模拟项目X\工作.gdb\stream_ras # D8 流向NORMAL 表示按常规地表流向编码 FlowDirection(dem_fill, NORMAL).save(flow_dir) # 流量累积默认按像元数累积 FlowAccumulation(flow_dir).save(flow_acc) # 阈值示例800 个像元实际按最小汇水面积除以像元面积换算 threshold 800 Con(arcpy.Raster(flow_acc) threshold, 1).save(stream_ras) print(栅格河网完成:, stream_ras)FlowDirection的NORMAL和FORCE影响边缘流向处理常规流域用NORMAL即可。FlowAccumulation默认输出像元数如果要换算成汇水面积再乘像元面积。threshold是像元数不是面积假设像元 30 米800 个像元约 0.72 平方公里这个量级要跟你的制图比例尺匹配。阈值调低河网变密伪河道增多阈值调高支流丢失。我一般至少跑三组阈值叠在 DEM 阴影上目视对比再定最终值。4.3 栅格河网转矢量栅格转线、平滑、去伪河道栅格河网转矢量后会出现锯齿和短碎线直接交付很难看。先用StreamToFeature把栅格河网和流向一起转成线保留流向信息再做平滑和长度筛选。下面脚本完成转线、加长度字段、按长度删短碎线。import arcpy arcpy.env.overwriteOutput True arcpy.CheckOutExtension(Spatial) stream_ras rF:\模拟项目X\工作.gdb\stream_ras flow_dir rF:\模拟项目X\工作.gdb\flow_dir stream_line rF:\模拟项目X\工作.gdb\stream_line # 栅格河网转线SIMPLIFY 会压缩折点 arcpy.sa.StreamToFeature(stream_ras, flow_dir, stream_line, SIMPLIFY) # 加长度字段并按米重算 arcpy.management.AddField(stream_line, LEN_M, DOUBLE) arcpy.management.CalculateGeometryAttributes( stream_line, [[LEN_M, LENGTH]], length_unitMETERS ) # 删除小于 500 米的短碎线阈值按出图比例尺调整 arcpy.management.SelectLayerByAttribute( stream_line, NEW_SELECTION, LEN_M 500 ) arcpy.management.DeleteRows(stream_line) print(矢量河网完成:, stream_line)StreamToFeature的第四个参数SIMPLIFY会简化线形减少折点如果后续要做精确长度统计可以改用NO_SIMPLIFY再单独平滑。LEN_M按米重算前提是当前坐标系为投影坐标系。删除短碎线前先备份一份未筛选的线避免把真实短支流删掉500 米只是示例实际按最小制图单元和河网密度调整。4.4 矢量河网与 DEM 一致性检查高程回读与坡度抽查矢量河网转出来后要确认它跟 DEM 对得上。常见做法是把 DEM 高程回读到河网折点上看同一河段是否沿流向下降如果出现大量逆坡说明填洼、流向或阈值有问题。下面脚本把河网转点并提取 DEM 值输出检查点。import arcpy arcpy.env.overwriteOutput True arcpy.CheckOutExtension(Spatial) stream_line rF:\模拟项目X\工作.gdb\stream_line dem_fill rF:\模拟项目X\工作.gdb\dem_fill check_pts rF:\模拟项目X\工作.gdb\stream_check # 线转点保留线属性便于按河段核对 arcpy.management.FeatureVerticesToPoints( stream_line, check_pts, ALL ) # 从填洼 DEM 回读高程到点 arcpy.sa.ExtractValuesToPoints( check_pts, dem_fill, check_pts, INTERPOLATE, VALUE_ONLY ) print(一致性检查点完成:, check_pts)FeatureVerticesToPoints的ALL会把所有折点转成点点太多时可以先抽稀。ExtractValuesToPoints的INTERPOLATE表示按像元中心插值VALUE_ONLY只把栅格值写入RASTERVALU字段。检查时按河段分组看高程变化出现连续逆坡就回到流向和填洼参数不要直接在矢量河网上改高程。5. 避坑与排查shp 中文乱码、投影错位、栅格 NoData、拓扑报错的现场处理5.1 五类高频翻车现场从现象倒推到参数第一类shp 属性表中文变问号或乱码。现象是字段名还能看字段值变成乱码原因是 DBF 代码页与当前环境不一致或者导出时没有指定编码。解决方式不是反复改字段名而是把中间数据放进文件地理数据库交付前再用“表转表”或“复制要素”指定代码页如果对方环境只认 shp先给一份 UTF-8 编码的 CSV 做对照再给 shp。第二类面要素与边界线错位几百米。现象是图层能叠在一起但边界整体偏移原因是坐标系未定义、误定义或者投影时少了地理变换参数。解决方式是先看.prj是否存在再用Project做坐标转换不要用DefineProjection硬改。已经错位的数据先找回原始坐标系再决定是否重投影地理变换参数要按数据所在区域选择不能随手默认。第三类DEM 裁剪后边缘出现黑边或异常低值。现象是填洼后边缘一大片平地流量累积在边缘堆高原因是 NoData 被当成 0 参与计算或者裁剪范围把外围汇水切断。解决方式是在栅格计算时用SetNull保留 NoData不要用 0 填充裁剪时用边界线外接矩形留一点缓冲填洼后再按边界裁最终成果。第四类拓扑报“must not overlap”但目视看不出重叠。现象是错误要素集中在分区边界原因是折点容差和坐标精度导致微重叠或者两条记录重复。解决方式是先跑“修复几何”再按坐标单位设置拓扑容差最后用错误表逐条核对微重叠可以用“整合”工具在容差内合并但整合前必须备份。第五类栅格转线后河道破碎、断头。现象是河网一段一段支流接不上干流原因是流量累积阈值过高、填洼不足、DEM 像元太大。解决方式是降低阈值做对比检查填洼是否覆盖整个汇水区必要时换更高分辨率 DEM转线后不要急着平滑先看断头位置是否落在真实分水岭上。5.2 排查顺序从范围到字段的十分钟清单数据出问题时我一般按“范围—坐标系—几何—字段—拓扑—栅格参数”的顺序排查不先从符号化找原因。范围对不上后面全是白搭坐标系没定义裁剪和投影都会偏几何类型不对工具直接报错字段名和编码不对属性统计会丢拓扑不干净面积和密度都不可信栅格像元、NoData、Z 单位不对水文派生全盘偏。下面这张表是我常用的十分钟清单按顺序跑一遍大多数问题能定位到具体环节。检查项工具或脚本通过标准失败处理范围Describe 读四至面、线、DEM 范围互相覆盖重投影或按外接矩形重裁坐标系看.prj与 Describe三者坐标系一致或可转换用 Project 转换不硬定义几何几何类型与修复几何面是面、线是线、无空几何修复几何后重建拓扑字段ListFields 看名称长度中文可读、字段名不超过 10 字符转 GDB 或改短名拓扑建拓扑并验证无重叠、无缝隙、无悬挂按错误表逐条修栅格像元、NoData、Z 单位像元一致、NoData 明确重采样或 SetNull这张表的价值在于顺序不在工具本身。先跑范围与坐标系能排除一半“图层对不上”的问题再跑几何与字段能排除裁剪报错和属性乱码最后跑拓扑与栅格参数才进入水文派生和面积统计。每次修改前留一份工作副本修改后用同一张表复验不要凭目视感觉“应该没问题”。6. 交付前验证与进阶用法用一套检查脚本把数据锁死6.1 交付前跑三项验证范围一致、字段可读、河网回落交付前我不看符号化效果只看三项硬指标范围是否落在边界线外接矩形内字段是否可读且长度合规矢量河网回读 DEM 后是否沿流向下降。三项都过再谈制图。下面这段短脚本只做范围与字段核对河网回落沿用第 4 章的检查点思路。import arcpy fc rF:\模拟项目X\交付\面要素_裁剪.shp boundary rF:\模拟项目X\工作.gdb\边界线_转面 # 范围核对面要素四至不应超出边界线外接矩形 f_ext arcpy.Describe(fc).extent b_ext arcpy.Describe(boundary).extent print(面范围:, f_ext.XMin, f_ext.YMin, f_ext.XMax, f_ext.YMax) print(边界范围:, b_ext.XMin, b_ext.YMin, b_ext.XMax, b_ext.YMax) # 字段核对中文可读、字段名不超过 10 字符 for f in arcpy.ListFields(fc): if f.type not in (Geometry, OID): print(f.name, len(f.name), f.type)Describe().extent读的是数据实际范围不是图层框字段循环排除几何和 OID只看业务字段。范围超出边界时先查裁剪要素是否用了转面后的边界而不是原始线。字段名超过 10 字符时回到 GDB 改短名再导出。6.2 一个习惯先只读体检再动裁剪我吃过一次亏拿到一套流域 shp 后直接裁面结果边界线没闭合裁出来的面缺了一块面积统计整整少了一个分区。从那以后我每次拿到这类流域 shp 数据都先跑只读体检和交付验证再动裁剪拓扑错误、字段长度、河网回落三项不过绝不进入制图。希望帮到你。本文还有配套的精品资源点击获取
