简介面向地理信息从业者、规划研究人员与高校师生这份江苏省GIS矢量数据包以Shapefile标准格式整理覆盖河流、湖泊、县域、地级市、高速公路、国道、铁路、湿地、丘陵及经纬网等要素可用于地图制图、空间叠加分析与区域研究帮助读者快速搭建江苏地理底图、省去逐项搜集与格式转换的麻烦。压缩包共83个文件约974KB以shp几何数据、dbf属性表、shx索引及prj投影文件为主另含sbx、sbn空间索引与xml元数据各要素配套齐全导入ArcGIS、QGIS即可直接读取。已有2386人学习下载说明其在教学与项目实践中具备一定参考价值。数据涵盖线状河流与路网、面状湖泊与行政边界便于开展洪水风险叠加、地区发展差异对比、交通网络优化等分析也可作为空间分析课程与规划项目的底图素材。1. 拿到一份省级水系与行政边界 shp 数据先想清楚它能干什么上周有个做水利信息化的朋友丢给我一个压缩包解压后是一堆.shp、.dbf、.shx、.prj文件目录名写着河流、湖泊、县域、高速几大类。他问的第一句话不是「怎么打开」而是「这玩意儿能直接进我的分析流程吗」。这个问题其实问到了点子上——省级尺度的河流、湖泊、县域、高速 shp 数据本质上是一套空间底图能不能用、好不好用取决于你对坐标系、属性字段、几何质量这三件事的判断而不是取决于文件能不能在 GIS 软件里显示出来。这类数据最常见的落地场景有四类一是做流域或行政区的水系连通性统计比如某个县境内有多少条河流、湖泊水面面积占比多少二是做高速路网与河湖水系的交叉分析比如桥梁选址、洪涝风险路段识别三是作为其他业务数据的空间参照底图比如把监测站点、排污口、工程点位叠上去做可视化四是做空间数据处理的练手素材练投影转换、裁剪、拓扑检查这些基本功。适合的人群也很明确做水利、交通、国土空间规划、应急管理相关 GIS 分析的人以及需要一份真实省级尺度数据来跑通处理流程的开发者。不适合指望它做高精度工程测量的人——省级公开尺度的数据几何精度和现势性都有天花板这一点后面会反复提到。2. 省级河流湖泊县域高速 shp 的坐标系与字段先读懂再动手2.1 为什么坐标系判断是第一步而不是最后一步很多人拿到 shp 的第一反应是双击打开看图形图形出来了就以为没问题。但省级尺度的数据几乎不可能是单一坐标系打包的——河流、湖泊这类自然要素常用地理坐标系比如 CGCS2000 经纬度县域边界可能也是经纬度而高速路网有时候会被处理成投影坐标系以便算长度。如果你不先确认每个图层的坐标系后面做面积统计、长度计算、缓冲区分析结果会错得离谱而且错得很隐蔽。判断方法很直接看.prj文件内容或者在 GIS 软件里看图层属性里的坐标参考信息。.prj是纯文本用记事本就能打开里面会写明是GCS_China_Geodetic_Coordinate_System_2000这类地理坐标系还是CGCS2000_3_Degree_GK_Zone_XX这类投影坐标系。我一般会写个小脚本批量读一遍比一个个点开快得多。import os from pathlib import Path # 批量读取 shp 对应的 prj 文件快速判断坐标系类型 def scan_prj(root_dir): results [] for prj_path in Path(root_dir).rglob(*.prj): shp_name prj_path.stem with open(prj_path, r, encodingutf-8, errorsignore) as f: content f.read() # 简单关键词判断GEOGCS 为地理坐标系PROJCS 为投影坐标系 if PROJCS in content: crs_type 投影坐标系 elif GEOGCS in content: crs_type 地理坐标系 else: crs_type 未知 results.append((shp_name, crs_type, content[:120])) return results for name, crs_type, snippet in scan_prj(./data): print(f{name}: {crs_type}) print(f 片段: {snippet}\n)这段代码的逻辑很朴素遍历目录下所有.prj读文本用PROJCS和GEOGCS两个关键字区分投影坐标系和地理坐标系。参数上唯一需要注意的是编码——有些老数据的 prj 文件不是 UTF-8所以加了errorsignore兜底。跑完之后你会得到一张清单哪些图层是经纬度、哪些是投影一目了然。这一步花五分钟能省掉后面几个小时的返工。2.2 属性字段怎么读dbf 里藏着数据的可用性边界shp 的几何只是骨架真正决定数据能不能用的是.dbf里的属性字段。河流图层通常会有名称、等级、长度这类字段湖泊图层会有名称、面积县域图层会有行政区代码和名称高速图层会有路线名称或编号。但省级公开数据的字段命名往往不统一可能是拼音缩写可能是中文也可能字段名被截断成 10 个字符以内dbf 格式的字段名长度限制。读字段最稳的方式是用geopandas或pyshp不要靠 GIS 软件里肉眼扫。下面这段代码把每个图层的字段名、类型、非空数量打出来顺便看看有没有关键字段缺失。import geopandas as gpd from pathlib import Path # 逐个读取 shp输出字段结构和记录数 def inspect_fields(shp_path): gdf gpd.read_file(shp_path) print(f图层: {Path(shp_path).name}) print(f记录数: {len(gdf)}, 几何类型: {gdf.geom_type.unique()}) print(f坐标系: {gdf.crs}) for col in gdf.columns: if col geometry: continue non_null gdf[col].notna().sum() print(f 字段 {col}: 类型{gdf[col].dtype}, 非空{non_null}) print(- * 50) for shp in Path(./data).rglob(*.shp): inspect_fields(shp)逻辑说明gpd.read_file直接吃 shp返回 GeoDataFramegeom_type.unique()能告诉你这个图层是线、面还是点混合几何类型在省级数据里并不罕见遇到要单独处理notna().sum()统计非空值如果某个关键字段非空数量远小于记录数说明这个字段大面积缺失不能作为分析主键。参数上没有太多可调的重点是养成「先看字段再动手」的习惯。我见过有人直接拿一个空字段去做分组统计跑完发现全是空值分组白忙一场。2.3 几何质量初检自相交、空几何、重复要素省级数据在合并、裁剪过程中很容易产生几何问题。最常见的是面要素自相交湖泊边界打结、空几何属性有记录但图形是空的、重复要素同一条河被录了两遍。这些问题在可视化时看不出来但一做空间叠加或面积计算就会报错或者出脏结果。用geopandas配合shapely做一轮初检成本很低import geopandas as gpd def check_geometry(shp_path): gdf gpd.read_file(shp_path) empty_count gdf.geometry.is_empty.sum() null_count gdf.geometry.isna().sum() invalid_count (~gdf.geometry.is_valid).sum() dup_count gdf.geometry.duplicated().sum() print(f{shp_path}) print(f 空几何: {empty_count}, 空值: {null_count}) print(f 无效几何: {invalid_count}, 重复几何: {dup_count}) for shp in [./data/rivers.shp, ./data/lakes.shp, ./data/county.shp, ./data/highway.shp]: check_geometry(shp)is_valid判断的是 OGC 几何有效性自相交的面会被标为无效duplicated()对几何做重复判断能抓出完全重合的要素。发现无效几何后常见做法是用buffer(0)或者make_valid修复但修复前一定要先搞清楚是数据本身的问题还是投影转换引入的问题否则修完还是错的。这一步的输出建议留档作为后续处理的基线。3. 用 geopandas 把河流湖泊县域高速跑通一次空间叠加3.1 统一坐标系投影转换的选与坑做任何跨图层的空间分析之前必须把所有图层统一到同一个坐标系。省级尺度的面积和长度统计我一般统一转到以米为单位的投影坐标系比如 CGCS2000 的高斯克吕格投影或者 Albers 等积投影。选哪个取决于你要算什么算长度用等距或高斯投影都行算面积优先用等积投影Albers因为等积投影能保证面积比例不失真。import geopandas as gpd # 读取各图层 rivers gpd.read_file(./data/rivers.shp) lakes gpd.read_file(./data/lakes.shp) county gpd.read_file(./data/county.shp) highway gpd.read_file(./data/highway.shp) # 统一转到 CGCS2000 3度带高斯投影以中央经线120度为例EPSG:4549 为常见带号 target_crs EPSG:4549 rivers_p rivers.to_crs(target_crs) lakes_p lakes.to_crs(target_crs) county_p county.to_crs(target_crs) highway_p highway.to_crs(target_crs) print(转换后坐标系:, rivers_p.crs) print(河流总长度(米):, rivers_p.length.sum()) print(湖泊总面积(平方米):, lakes_p.area.sum())逻辑说明to_crs是原地转换返回新的 GeoDataFramelength和area只有在投影坐标系下才有物理意义经纬度下算出来的是度数没有参考价值。参数上EPSG:4549只是举例实际用哪个带号要看数据覆盖的经度范围跨带的话要么分带处理要么换 Albers 这类大范围投影。这里最容易翻车的地方是有人把经纬度数据直接.length求和得到一个几十万的数字还以为是米实际是度差了十万八千里。3.2 县域裁剪河流湖泊空间叠加的标准流程统一坐标系之后最典型的操作是用县域边界去裁剪河流和湖泊统计每个县境内的水系情况。标准流程是「空间连接 分组统计」或者「裁剪 按县聚合」前者保留要素完整性后者会切断跨县要素。做统计我一般用空间连接因为一条河跨两个县按县统计长度时应该各算各的。# 河流与县域做空间连接统计每县河流条数与总长度 rivers_join gpd.sjoin(rivers_p, county_p[[县名, geometry]], howinner, predicateintersects) # 按县聚合 county_river_stats rivers_join.groupby(县名).agg( 河流条数(geometry, count), ).reset_index() # 长度需要先算到要素上再聚合 rivers_join[seg_length] rivers_join.length length_stats rivers_join.groupby(县名)[seg_length].sum().reset_index() length_stats.columns [县名, 河流总长度_米] result county_river_stats.merge(length_stats, on县名) print(result.head())逻辑说明sjoin的predicateintersects表示只要河流与县域有交集就关联上howinner只保留能匹配上的记录。注意这里统计的长度是整条河流的长度不是县内那一段的长度——如果要算县内实际长度得先做overlay裁剪再算。这是一个常见的口径差异做报告前一定要跟需求方确认清楚。参数上predicate还可以选within、contains含义不同选错会导致漏统计或重复统计。3.3 高速与河湖交叉分析缓冲区与相交判断高速路网和水系的交叉分析典型需求是找出高速跨越河流湖泊的位置用于桥梁核查或洪涝风险排查。做法是给河流湖泊做缓冲区再和高速做相交判断。缓冲区半径根据业务定一般河流取 50 到 200 米湖泊可以大一些。# 河流做100米缓冲区与高速相交找出交叉路段 rivers_buffer rivers_p.copy() rivers_buffer[geometry] rivers_p.buffer(100) cross gpd.sjoin(highway_p, rivers_buffer[[geometry]], howinner, predicateintersects) cross_unique cross.drop_duplicates(subset[geometry]) print(f与河流100米缓冲区相交的高速要素数: {len(cross_unique)})逻辑说明buffer(100)生成 100 米缓冲面sjoin判断高速与缓冲面是否相交。drop_duplicates是因为一条高速可能同时与多条河流的缓冲区相交会产生重复记录。参数上缓冲区半径是核心太小会漏掉实际交叉太大会引入大量误判建议先用几个半径值跑一遍看结果数量变化再定最终值。这里还有个坑如果高速和河流坐标系不一致就做 buffer半径单位会错所以前面统一坐标系这一步不能省。4. 省级 shp 数据处理的避坑清单五个真实翻车现场4.1 现象面积统计结果大得离谱原因经纬度下直接算面积第一次拿到湖泊数据直接gdf.area.sum()得到一个几万的数字以为是平方公里实际是平方度。经纬度坐标系下area返回的是度平方没有任何物理意义。解决办法是先to_crs转到等积投影再算算完再按需换算单位。这个坑几乎每个新手都会踩一次血泪经验就是看到面积或长度结果先问一句「当前坐标系是什么」。4.2 现象sjoin 之后记录数暴涨原因一对多关联没去重用县域去关联河流时如果一条河流与多个县相交sjoin会为每个县生成一条记录导致河流被重复计数。做条数统计时这没问题本来就该各算各的但做「河流总数」这种全局统计时就会翻倍。解决办法是统计全局时用原始图层统计分县时用关联结果两个口径分开算别混用。4.3 现象裁剪后要素消失原因几何无效导致 overlay 失败用overlay做裁剪时如果参与运算的面要素有自相交结果可能为空或者报拓扑异常。排查方法是先跑一遍is_valid检查对无效几何做make_valid修复再裁剪。修复后建议对比修复前后的面积差异差异过大说明原数据问题严重要谨慎使用。4.4 现象中文属性乱码原因dbf 编码不是 UTF-8打开 dbf 发现中文字段全是乱码这是编码问题。dbf 常见编码有 GBK、GB2312、UTF-8 几种geopandas读取时可以指定encoding参数。如果读进来是乱码试试encodinggbk或encodingutf-8还不行就用pyshp读原始字节判断。这个坑不影响几何但会让属性关联全部失效属于必须提前确认的事项。4.5 现象不同来源图层叠加错位原因坐标系标注错误而非转换错误两个图层都显示是 CGCS2000但叠在一起偏了几百米。这种情况往往是某个图层的.prj标注和实际坐标不符——数据被处理过但 prj 没更新。排查方法是找一个已知地物比如某个明显的湖泊轮廓在两个图层里对比位置偏差规律如果是固定平移基本可以确认是标注问题。解决办法是手动指定正确的坐标系再转换而不是盲目信任 prj 文件。5. 把省级 shp 做成可复用底图字段标准化与批量导出5.1 字段标准化给四个图层定一套统一命名原始数据的字段命名往往五花八门直接拿去做分析每换一个图层就要改一次代码。我的习惯是先做一轮字段标准化把四个图层的核心字段统一成英文短名比如name、code、length_m、area_m2这样后续脚本可以复用同一套逻辑。标准化之前先做字段映射表明确哪个原字段对应哪个新字段映射关系写进配置不要硬编码在代码里。import geopandas as gpd # 字段映射配置图层名 - {原字段: 新字段} field_map { rivers: {河流名称: name, 河流等级: grade}, lakes: {湖泊名称: name, 水面面积: area_raw}, county: {县名: name, 行政区代码: code}, highway: {路线名称: name, 路线编号: code}, } def standardize_fields(shp_path, mapping): gdf gpd.read_file(shp_path) # 只保留映射表中存在的字段避免 KeyError rename_dict {k: v for k, v in mapping.items() if k in gdf.columns} gdf gdf.rename(columnsrename_dict) return gdf rivers_std standardize_fields(./data/rivers.shp, field_map[rivers]) print(rivers_std.columns.tolist())逻辑说明映射表用字典配置rename前先过滤不存在的字段避免因为某个图层缺字段导致整个流程中断。参数上没什么可调的关键是映射表要跟业务方对齐别自己拍脑袋改名。标准化之后四个图层的name字段含义一致后续做联合查询或者统一出图就顺了。5.2 批量导出 GeoPackage比一堆 shp 更好管shp 格式有一堆伴生文件四个图层就是十几个文件传输和版本管理都麻烦。我一般会把标准化后的图层批量写进一个 GeoPackage.gpkg单文件、支持多图层、字段名长度不受限比 shp 省心。导出时注意图层名不要用中文避免某些工具兼容问题。import geopandas as gpd from pathlib import Path layers { rivers_std: rivers_std, lakes_std: standardize_fields(./data/lakes.shp, field_map[lakes]), county_std: standardize_fields(./data/county.shp, field_map[county]), highway_std: standardize_fields(./data/highway.shp, field_map[highway]), } out_gpkg ./output/base_map.gpkg if Path(out_gpkg).exists(): Path(out_gpkg).unlink() # 覆盖前先删除避免追加导致图层重复 for layer_name, gdf in layers.items(): gdf.to_file(out_gpkg, layerlayer_name, driverGPKG) print(f已写入图层: {layer_name}, 要素数: {len(gdf)})逻辑说明to_file指定driverGPKG写入 GeoPackagelayer参数决定图层名。覆盖写入前先删文件是因为 GPKG 的追加模式容易产生同名图层冲突直接重建最干净。参数上driver不要写错写成ESRI Shapefile就退回多文件模式了。导出后建议用fiona.listlayers或者重新读一遍确认图层完整别导完就不管了。5.3 一个可复用的处理入口把流程串成函数把前面的步骤串成一个入口函数输入原始目录输出标准化后的 GeoPackage 和一份质检报告。这样下次拿到同类型数据改一下路径就能跑不用从头写。质检报告至少包含每个图层的记录数、坐标系、无效几何数、字段缺失情况出问题时有据可查。def build_base_map(raw_dir, out_gpkg, target_crsEPSG:4549): report [] for shp in Path(raw_dir).rglob(*.shp): gdf gpd.read_file(shp) invalid (~gdf.geometry.is_valid).sum() gdf_proj gdf.to_crs(target_crs) report.append({ layer: shp.stem, records: len(gdf), invalid_geom: int(invalid), crs: str(gdf.crs), }) gdf_proj.to_file(out_gpkg, layershp.stem, driverGPKG) return report report build_base_map(./data, ./output/base_map_v2.gpkg) for r in report: print(r)逻辑说明这个函数做了三件事——读、检、转、写质检信息收集在report里返回。参数target_crs留了默认值但可覆盖方便不同项目切换投影。注意这里没有做字段标准化因为标准化依赖映射表属于业务层配置不适合硬塞进通用函数。实际使用时我会把标准化和这个入口分开通用流程管几何和坐标系业务配置管字段职责清晰改起来不互相牵连。5.4 验证导出结果三个必查项导出之后别急着交付至少查三件事。第一图层数量和要素数是否和原始一致防止写入过程中丢图层第二坐标系是否统一为目标投影防止某个图层漏转第三随机抽几个要素对比原始和导出后的几何范围total_bounds偏差应该在容差范围内。这三项查完基本可以确认导出没有引入新问题。我一般会把这套检查写成脚本每次导出自动跑比人工点开看靠谱。省级数据体量不小人工抽查容易漏脚本化是唯一可持续的办法。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/11/4 16:41:01 · 阅读 1.2k 省级河流湖泊县域高速shp数据处理：坐标系、字段与空间分析实战
简介面向地理信息从业者、规划研究人员与高校师生这份江苏省GIS矢量数据包以Shapefile标准格式整理覆盖河流、湖泊、县域、地级市、高速公路、国道、铁路、湿地、丘陵及经纬网等要素可用于地图制图、空间叠加分析与区域研究帮助读者快速搭建江苏地理底图、省去逐项搜集与格式转换的麻烦。压缩包共83个文件约974KB以shp几何数据、dbf属性表、shx索引及prj投影文件为主另含sbx、sbn空间索引与xml元数据各要素配套齐全导入ArcGIS、QGIS即可直接读取。已有2386人学习下载说明其在教学与项目实践中具备一定参考价值。数据涵盖线状河流与路网、面状湖泊与行政边界便于开展洪水风险叠加、地区发展差异对比、交通网络优化等分析也可作为空间分析课程与规划项目的底图素材。1. 拿到一份省级水系与行政边界 shp 数据先想清楚它能干什么上周有个做水利信息化的朋友丢给我一个压缩包解压后是一堆.shp、.dbf、.shx、.prj文件目录名写着河流、湖泊、县域、高速几大类。他问的第一句话不是「怎么打开」而是「这玩意儿能直接进我的分析流程吗」。这个问题其实问到了点子上——省级尺度的河流、湖泊、县域、高速 shp 数据本质上是一套空间底图能不能用、好不好用取决于你对坐标系、属性字段、几何质量这三件事的判断而不是取决于文件能不能在 GIS 软件里显示出来。这类数据最常见的落地场景有四类一是做流域或行政区的水系连通性统计比如某个县境内有多少条河流、湖泊水面面积占比多少二是做高速路网与河湖水系的交叉分析比如桥梁选址、洪涝风险路段识别三是作为其他业务数据的空间参照底图比如把监测站点、排污口、工程点位叠上去做可视化四是做空间数据处理的练手素材练投影转换、裁剪、拓扑检查这些基本功。适合的人群也很明确做水利、交通、国土空间规划、应急管理相关 GIS 分析的人以及需要一份真实省级尺度数据来跑通处理流程的开发者。不适合指望它做高精度工程测量的人——省级公开尺度的数据几何精度和现势性都有天花板这一点后面会反复提到。2. 省级河流湖泊县域高速 shp 的坐标系与字段先读懂再动手2.1 为什么坐标系判断是第一步而不是最后一步很多人拿到 shp 的第一反应是双击打开看图形图形出来了就以为没问题。但省级尺度的数据几乎不可能是单一坐标系打包的——河流、湖泊这类自然要素常用地理坐标系比如 CGCS2000 经纬度县域边界可能也是经纬度而高速路网有时候会被处理成投影坐标系以便算长度。如果你不先确认每个图层的坐标系后面做面积统计、长度计算、缓冲区分析结果会错得离谱而且错得很隐蔽。判断方法很直接看.prj文件内容或者在 GIS 软件里看图层属性里的坐标参考信息。.prj是纯文本用记事本就能打开里面会写明是GCS_China_Geodetic_Coordinate_System_2000这类地理坐标系还是CGCS2000_3_Degree_GK_Zone_XX这类投影坐标系。我一般会写个小脚本批量读一遍比一个个点开快得多。import os from pathlib import Path # 批量读取 shp 对应的 prj 文件快速判断坐标系类型 def scan_prj(root_dir): results [] for prj_path in Path(root_dir).rglob(*.prj): shp_name prj_path.stem with open(prj_path, r, encodingutf-8, errorsignore) as f: content f.read() # 简单关键词判断GEOGCS 为地理坐标系PROJCS 为投影坐标系 if PROJCS in content: crs_type 投影坐标系 elif GEOGCS in content: crs_type 地理坐标系 else: crs_type 未知 results.append((shp_name, crs_type, content[:120])) return results for name, crs_type, snippet in scan_prj(./data): print(f{name}: {crs_type}) print(f 片段: {snippet}\n)这段代码的逻辑很朴素遍历目录下所有.prj读文本用PROJCS和GEOGCS两个关键字区分投影坐标系和地理坐标系。参数上唯一需要注意的是编码——有些老数据的 prj 文件不是 UTF-8所以加了errorsignore兜底。跑完之后你会得到一张清单哪些图层是经纬度、哪些是投影一目了然。这一步花五分钟能省掉后面几个小时的返工。2.2 属性字段怎么读dbf 里藏着数据的可用性边界shp 的几何只是骨架真正决定数据能不能用的是.dbf里的属性字段。河流图层通常会有名称、等级、长度这类字段湖泊图层会有名称、面积县域图层会有行政区代码和名称高速图层会有路线名称或编号。但省级公开数据的字段命名往往不统一可能是拼音缩写可能是中文也可能字段名被截断成 10 个字符以内dbf 格式的字段名长度限制。读字段最稳的方式是用geopandas或pyshp不要靠 GIS 软件里肉眼扫。下面这段代码把每个图层的字段名、类型、非空数量打出来顺便看看有没有关键字段缺失。import geopandas as gpd from pathlib import Path # 逐个读取 shp输出字段结构和记录数 def inspect_fields(shp_path): gdf gpd.read_file(shp_path) print(f图层: {Path(shp_path).name}) print(f记录数: {len(gdf)}, 几何类型: {gdf.geom_type.unique()}) print(f坐标系: {gdf.crs}) for col in gdf.columns: if col geometry: continue non_null gdf[col].notna().sum() print(f 字段 {col}: 类型{gdf[col].dtype}, 非空{non_null}) print(- * 50) for shp in Path(./data).rglob(*.shp): inspect_fields(shp)逻辑说明gpd.read_file直接吃 shp返回 GeoDataFramegeom_type.unique()能告诉你这个图层是线、面还是点混合几何类型在省级数据里并不罕见遇到要单独处理notna().sum()统计非空值如果某个关键字段非空数量远小于记录数说明这个字段大面积缺失不能作为分析主键。参数上没有太多可调的重点是养成「先看字段再动手」的习惯。我见过有人直接拿一个空字段去做分组统计跑完发现全是空值分组白忙一场。2.3 几何质量初检自相交、空几何、重复要素省级数据在合并、裁剪过程中很容易产生几何问题。最常见的是面要素自相交湖泊边界打结、空几何属性有记录但图形是空的、重复要素同一条河被录了两遍。这些问题在可视化时看不出来但一做空间叠加或面积计算就会报错或者出脏结果。用geopandas配合shapely做一轮初检成本很低import geopandas as gpd def check_geometry(shp_path): gdf gpd.read_file(shp_path) empty_count gdf.geometry.is_empty.sum() null_count gdf.geometry.isna().sum() invalid_count (~gdf.geometry.is_valid).sum() dup_count gdf.geometry.duplicated().sum() print(f{shp_path}) print(f 空几何: {empty_count}, 空值: {null_count}) print(f 无效几何: {invalid_count}, 重复几何: {dup_count}) for shp in [./data/rivers.shp, ./data/lakes.shp, ./data/county.shp, ./data/highway.shp]: check_geometry(shp)is_valid判断的是 OGC 几何有效性自相交的面会被标为无效duplicated()对几何做重复判断能抓出完全重合的要素。发现无效几何后常见做法是用buffer(0)或者make_valid修复但修复前一定要先搞清楚是数据本身的问题还是投影转换引入的问题否则修完还是错的。这一步的输出建议留档作为后续处理的基线。3. 用 geopandas 把河流湖泊县域高速跑通一次空间叠加3.1 统一坐标系投影转换的选与坑做任何跨图层的空间分析之前必须把所有图层统一到同一个坐标系。省级尺度的面积和长度统计我一般统一转到以米为单位的投影坐标系比如 CGCS2000 的高斯克吕格投影或者 Albers 等积投影。选哪个取决于你要算什么算长度用等距或高斯投影都行算面积优先用等积投影Albers因为等积投影能保证面积比例不失真。import geopandas as gpd # 读取各图层 rivers gpd.read_file(./data/rivers.shp) lakes gpd.read_file(./data/lakes.shp) county gpd.read_file(./data/county.shp) highway gpd.read_file(./data/highway.shp) # 统一转到 CGCS2000 3度带高斯投影以中央经线120度为例EPSG:4549 为常见带号 target_crs EPSG:4549 rivers_p rivers.to_crs(target_crs) lakes_p lakes.to_crs(target_crs) county_p county.to_crs(target_crs) highway_p highway.to_crs(target_crs) print(转换后坐标系:, rivers_p.crs) print(河流总长度(米):, rivers_p.length.sum()) print(湖泊总面积(平方米):, lakes_p.area.sum())逻辑说明to_crs是原地转换返回新的 GeoDataFramelength和area只有在投影坐标系下才有物理意义经纬度下算出来的是度数没有参考价值。参数上EPSG:4549只是举例实际用哪个带号要看数据覆盖的经度范围跨带的话要么分带处理要么换 Albers 这类大范围投影。这里最容易翻车的地方是有人把经纬度数据直接.length求和得到一个几十万的数字还以为是米实际是度差了十万八千里。3.2 县域裁剪河流湖泊空间叠加的标准流程统一坐标系之后最典型的操作是用县域边界去裁剪河流和湖泊统计每个县境内的水系情况。标准流程是「空间连接 分组统计」或者「裁剪 按县聚合」前者保留要素完整性后者会切断跨县要素。做统计我一般用空间连接因为一条河跨两个县按县统计长度时应该各算各的。# 河流与县域做空间连接统计每县河流条数与总长度 rivers_join gpd.sjoin(rivers_p, county_p[[县名, geometry]], howinner, predicateintersects) # 按县聚合 county_river_stats rivers_join.groupby(县名).agg( 河流条数(geometry, count), ).reset_index() # 长度需要先算到要素上再聚合 rivers_join[seg_length] rivers_join.length length_stats rivers_join.groupby(县名)[seg_length].sum().reset_index() length_stats.columns [县名, 河流总长度_米] result county_river_stats.merge(length_stats, on县名) print(result.head())逻辑说明sjoin的predicateintersects表示只要河流与县域有交集就关联上howinner只保留能匹配上的记录。注意这里统计的长度是整条河流的长度不是县内那一段的长度——如果要算县内实际长度得先做overlay裁剪再算。这是一个常见的口径差异做报告前一定要跟需求方确认清楚。参数上predicate还可以选within、contains含义不同选错会导致漏统计或重复统计。3.3 高速与河湖交叉分析缓冲区与相交判断高速路网和水系的交叉分析典型需求是找出高速跨越河流湖泊的位置用于桥梁核查或洪涝风险排查。做法是给河流湖泊做缓冲区再和高速做相交判断。缓冲区半径根据业务定一般河流取 50 到 200 米湖泊可以大一些。# 河流做100米缓冲区与高速相交找出交叉路段 rivers_buffer rivers_p.copy() rivers_buffer[geometry] rivers_p.buffer(100) cross gpd.sjoin(highway_p, rivers_buffer[[geometry]], howinner, predicateintersects) cross_unique cross.drop_duplicates(subset[geometry]) print(f与河流100米缓冲区相交的高速要素数: {len(cross_unique)})逻辑说明buffer(100)生成 100 米缓冲面sjoin判断高速与缓冲面是否相交。drop_duplicates是因为一条高速可能同时与多条河流的缓冲区相交会产生重复记录。参数上缓冲区半径是核心太小会漏掉实际交叉太大会引入大量误判建议先用几个半径值跑一遍看结果数量变化再定最终值。这里还有个坑如果高速和河流坐标系不一致就做 buffer半径单位会错所以前面统一坐标系这一步不能省。4. 省级 shp 数据处理的避坑清单五个真实翻车现场4.1 现象面积统计结果大得离谱原因经纬度下直接算面积第一次拿到湖泊数据直接gdf.area.sum()得到一个几万的数字以为是平方公里实际是平方度。经纬度坐标系下area返回的是度平方没有任何物理意义。解决办法是先to_crs转到等积投影再算算完再按需换算单位。这个坑几乎每个新手都会踩一次血泪经验就是看到面积或长度结果先问一句「当前坐标系是什么」。4.2 现象sjoin 之后记录数暴涨原因一对多关联没去重用县域去关联河流时如果一条河流与多个县相交sjoin会为每个县生成一条记录导致河流被重复计数。做条数统计时这没问题本来就该各算各的但做「河流总数」这种全局统计时就会翻倍。解决办法是统计全局时用原始图层统计分县时用关联结果两个口径分开算别混用。4.3 现象裁剪后要素消失原因几何无效导致 overlay 失败用overlay做裁剪时如果参与运算的面要素有自相交结果可能为空或者报拓扑异常。排查方法是先跑一遍is_valid检查对无效几何做make_valid修复再裁剪。修复后建议对比修复前后的面积差异差异过大说明原数据问题严重要谨慎使用。4.4 现象中文属性乱码原因dbf 编码不是 UTF-8打开 dbf 发现中文字段全是乱码这是编码问题。dbf 常见编码有 GBK、GB2312、UTF-8 几种geopandas读取时可以指定encoding参数。如果读进来是乱码试试encodinggbk或encodingutf-8还不行就用pyshp读原始字节判断。这个坑不影响几何但会让属性关联全部失效属于必须提前确认的事项。4.5 现象不同来源图层叠加错位原因坐标系标注错误而非转换错误两个图层都显示是 CGCS2000但叠在一起偏了几百米。这种情况往往是某个图层的.prj标注和实际坐标不符——数据被处理过但 prj 没更新。排查方法是找一个已知地物比如某个明显的湖泊轮廓在两个图层里对比位置偏差规律如果是固定平移基本可以确认是标注问题。解决办法是手动指定正确的坐标系再转换而不是盲目信任 prj 文件。5. 把省级 shp 做成可复用底图字段标准化与批量导出5.1 字段标准化给四个图层定一套统一命名原始数据的字段命名往往五花八门直接拿去做分析每换一个图层就要改一次代码。我的习惯是先做一轮字段标准化把四个图层的核心字段统一成英文短名比如name、code、length_m、area_m2这样后续脚本可以复用同一套逻辑。标准化之前先做字段映射表明确哪个原字段对应哪个新字段映射关系写进配置不要硬编码在代码里。import geopandas as gpd # 字段映射配置图层名 - {原字段: 新字段} field_map { rivers: {河流名称: name, 河流等级: grade}, lakes: {湖泊名称: name, 水面面积: area_raw}, county: {县名: name, 行政区代码: code}, highway: {路线名称: name, 路线编号: code}, } def standardize_fields(shp_path, mapping): gdf gpd.read_file(shp_path) # 只保留映射表中存在的字段避免 KeyError rename_dict {k: v for k, v in mapping.items() if k in gdf.columns} gdf gdf.rename(columnsrename_dict) return gdf rivers_std standardize_fields(./data/rivers.shp, field_map[rivers]) print(rivers_std.columns.tolist())逻辑说明映射表用字典配置rename前先过滤不存在的字段避免因为某个图层缺字段导致整个流程中断。参数上没什么可调的关键是映射表要跟业务方对齐别自己拍脑袋改名。标准化之后四个图层的name字段含义一致后续做联合查询或者统一出图就顺了。5.2 批量导出 GeoPackage比一堆 shp 更好管shp 格式有一堆伴生文件四个图层就是十几个文件传输和版本管理都麻烦。我一般会把标准化后的图层批量写进一个 GeoPackage.gpkg单文件、支持多图层、字段名长度不受限比 shp 省心。导出时注意图层名不要用中文避免某些工具兼容问题。import geopandas as gpd from pathlib import Path layers { rivers_std: rivers_std, lakes_std: standardize_fields(./data/lakes.shp, field_map[lakes]), county_std: standardize_fields(./data/county.shp, field_map[county]), highway_std: standardize_fields(./data/highway.shp, field_map[highway]), } out_gpkg ./output/base_map.gpkg if Path(out_gpkg).exists(): Path(out_gpkg).unlink() # 覆盖前先删除避免追加导致图层重复 for layer_name, gdf in layers.items(): gdf.to_file(out_gpkg, layerlayer_name, driverGPKG) print(f已写入图层: {layer_name}, 要素数: {len(gdf)})逻辑说明to_file指定driverGPKG写入 GeoPackagelayer参数决定图层名。覆盖写入前先删文件是因为 GPKG 的追加模式容易产生同名图层冲突直接重建最干净。参数上driver不要写错写成ESRI Shapefile就退回多文件模式了。导出后建议用fiona.listlayers或者重新读一遍确认图层完整别导完就不管了。5.3 一个可复用的处理入口把流程串成函数把前面的步骤串成一个入口函数输入原始目录输出标准化后的 GeoPackage 和一份质检报告。这样下次拿到同类型数据改一下路径就能跑不用从头写。质检报告至少包含每个图层的记录数、坐标系、无效几何数、字段缺失情况出问题时有据可查。def build_base_map(raw_dir, out_gpkg, target_crsEPSG:4549): report [] for shp in Path(raw_dir).rglob(*.shp): gdf gpd.read_file(shp) invalid (~gdf.geometry.is_valid).sum() gdf_proj gdf.to_crs(target_crs) report.append({ layer: shp.stem, records: len(gdf), invalid_geom: int(invalid), crs: str(gdf.crs), }) gdf_proj.to_file(out_gpkg, layershp.stem, driverGPKG) return report report build_base_map(./data, ./output/base_map_v2.gpkg) for r in report: print(r)逻辑说明这个函数做了三件事——读、检、转、写质检信息收集在report里返回。参数target_crs留了默认值但可覆盖方便不同项目切换投影。注意这里没有做字段标准化因为标准化依赖映射表属于业务层配置不适合硬塞进通用函数。实际使用时我会把标准化和这个入口分开通用流程管几何和坐标系业务配置管字段职责清晰改起来不互相牵连。5.4 验证导出结果三个必查项导出之后别急着交付至少查三件事。第一图层数量和要素数是否和原始一致防止写入过程中丢图层第二坐标系是否统一为目标投影防止某个图层漏转第三随机抽几个要素对比原始和导出后的几何范围total_bounds偏差应该在容差范围内。这三项查完基本可以确认导出没有引入新问题。我一般会把这套检查写成脚本每次导出自动跑比人工点开看靠谱。省级数据体量不小人工抽查容易漏脚本化是唯一可持续的办法。希望帮到你。本文还有配套的精品资源点击获取