SHP区划与路网GIS分析：坐标、叠加与密度避坑指南

简介面向GIS学习与空间分析场景这份江苏省地理数据包整合了省级、地市级、区县级三级行政区划边界以及道路网和铁路网矢量图层适合从事国土空间规划、交通网络研究、地图制图与课程实验的师生及技术人员使用可解决区域边界底图缺失、路网数据分散等问题。压缩包共33个文件约41.72MB以shp、dbf、shx、prj等Shapefile核心文件为主并配有xml元数据与sbn、sbx空间索引文件便于在ArcGIS、QGIS等平台直接加载、查询与叠加分析。已有1973人学习下载具备一定参考热度。数据覆盖省、市、县多级尺度道路与铁路要素可与行政边界组合用于制图表达、可达性分析、路网密度统计及空间分布展示规范的文件结构也方便按图层检索、裁剪和二次处理适合作为区域研究的基础数据底稿。1. 一套 SHP 顶三个数据源省域三级区划加双路网先解决什么做江苏范围内的 GIS 分析最耗时间的从来不是算法而是底图。行政区划要省、市、县三级道路网要分等级公路网经常还要单独成层三个来源、三套坐标系、三份更新时间光对齐口径就能耗掉半天。这份 2022.07 的 SHP 包把省、市、县三级行政区划和道路网、公路网放在同一批文件里最直接的价值是省掉了“对不齐”的活区划面自带行政代码可以当主键路网线可以直接和面对叠。它适合做区域统计、路网密度核算、专题制图底图、空间连接类练习也适合当教学或内部演示的固定样本。不适合的场合也说清楚需要月度级路网更新、需要精细到车道或收费站节点的这批数据不够用得另找源。下面按拆包、加载、叠加、排错、固化五步走把每一处能翻车的地方都标出来。2. 拆包先看结构SHP 的五个伴生文件与三级区划属性表拿到压缩包先别急着往软件里拖。Shapefile 是上世纪留下来的格式规范它对“一个图层”的定义和大多数人的直觉不一样对字段名、字符长度、单文件体积都有一堆硬约束。先把结构摸清楚后面九成的玄学问题都能提前避开。2.1 一个图层不是一条 .shp伴生文件缺一个就出事解压后你会看到同一前缀下跟着一串后缀。真正承载体的是.shp但只靠它打不开完整的图层。.shx是空间索引缺了它绝大多数库会直接拒绝读取或者全表扫描.dbf存属性表缺了它图形还在但一个字段都没有.prj记录坐标参考缺了它软件默认按未知坐标系处理你后面算出来的面积长度全是错的.cpg声明字符编码这个文件是中文乱码的头号元凶很多资源包干脆不带。先跑一次元数据体检把要素数、几何类型、字段清单、坐标范围一次性打出来# -so 只输出摘要-al 覆盖所有图层避免逐层敲命令 ogrinfo -so -al ./js_2022/county.shp # 单独确认 .prj 声明的是什么坐标系 gdalsrsinfo ./js_2022/county.shp # 看字符编码声明有 .cpg 才会输出 cat ./js_2022/county.cpgogrinfo -so输出的Extent那两行要重点看数值落在 -180~180 和 -90~90 之间说明是地理坐标单位是度数值是六位、七位数说明已经是投影坐标单位是米。这一步判断错了后面所有面积统计都会差出好几个数量级。gdalsrsinfo会打印 EPSG 编号或 WKT国内 2022 年前后的区划数据常见做法是 CGCS2000 地理坐标EPSG 为 4490但不要默认一定要读出来确认。如果.cpg缺失或者内容写的是 UTF-8 而实际数据是 GBK属性表里的中文会变成问号或方块。这种情况不要在原文件上反复折腾直接转成 GeoPackage 一次性解决编码和字段长度的双重限制# -lco ENCODINGUTF-8 强制以 UTF-8 落库-nln 指定新图层名 ogr2ogr -f GPKG ./js_2022.gpkg ./js_2022/county.shp \ -nln county -lco ENCODINGUTF-8 -a_srs EPSG:4490-a_srs是给那些.prj丢失或写错的文件强行指定坐标系用之前必须确认数据的真实基准否则等于给错误答案盖了个正确的章。GeoPackage 支持任意长度字段名和 UTF-8DBF 那两个坑字段名最多 10 字节、字符字段最多 254 字节在 GPKG 里都不存在长文本备注类字段一定要走这一步。2.2 属性字段怎么读区划代码是唯一可靠的主键三级区划图层最容易踩的坑是把“名称”当关联键。名称在不同版本里可能有简称、全称、带“市”不带“市”的差异甚至同一个名字在不同层级重复出现。真正稳的是行政区划代码省级 2 位补 6 位、地市级 4 位补 6 位、县级 6 位完整码。代码本身有层级含义前 2 位是省前 4 位是市做上下级聚合时直接截字符串就行不用做任何名称匹配。典型字段结构大致如下字段名以实际.dbf为准不同来源命名习惯差别很大常见字段名类型含义使用建议ADCODE / CODE / XZQDM字符行政区划代码首选主键读入后统一补零到 6 位NAME / 名称 / XZQMC字符行政区名称仅用于展示不用于关联LEVEL / 层级字符或数值省/市/县标识三级混在一个文件时用它拆层面积类字段数值数据源方计算的面积口径不明建议自己重算SHAPE_Leng / SHAPE_Area数值制图软件生成的几何量单位跟随坐标系不可直接使用表格里最后两行是最容易被误用的。数据源方给的面积字段你不知道它是在什么投影下算的、有没有扣掉水域、是不是投影变形区的结果。自己用等积投影重算一遍几十行代码的事比事后解释数据对不上要省心得多。读入后第一件事是规范化主键import geopandas as gpd county gpd.read_file(./js_2022.gpkg, layercounty, encodingutf-8) # 候选字段名按优先级找找不到直接报错避免静默用错列 def pick(gdf, candidates): for c in candidates: if c in gdf.columns: return c raise KeyError(f字段缺失候选: {candidates}实际: {list(gdf.columns)}) code_col pick(county, [ADCODE, CODE, XZQDM]) name_col pick(county, [NAME, 名称, XZQMC]) # 补零到 6 位ArcGIS 导出时前导零常被吃掉变成整数 county[adcode] county[code_col].astype(str).str.strip().str.zfill(6) county[name] county[name_col].astype(str).str.strip() # 检查一遍长度分布非 6 位的行要单独看 print(county[adcode].str.len().value_counts())astype(str)必须在zfill之前做因为很多资源包在导出时把代码列存成了数值型110000会变成110000.0直接补零会得到一串带小数点的废码。str.len().value_counts()打印出来的分布如果全是 6说明主键干净混进 5 位、7 位或者nan的多半是空几何行或者汇总行要单独剔掉再往下走。2.3 坐标系判断.prj 里那一行决定后面所有数字对不对国内区划和路网数据常见的坐标系有三类用错一个结果就废坐标系典型 EPSG单位适用场景陷阱CGCS2000 地理坐标4490度数据存储、跨带拼接直接算面积长度完全无意义CGCS2000 3 度带高斯克吕格4527 / 4528 等米单带范围内的精细量算跨带后东侧变形急剧放大中国 Albers 等积投影自定义 proj 串米全省级面积统计、密度核算需要手动构造无标准 EPSG第三行是全省统计的正解。江苏东西跨度大用单一 3 度带比如中央经线 117°E 或 120°E覆盖全省边缘地带的长度和面积都会有可观偏差用 Albers 等积投影面积是严格保持的长度虽然有变形但全省范围内可以接受做密度类指标完全够用。# 中国标准 Albers 等积投影GRS80 椭球对应 CGCS2000 基准 ALBERS (projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 ellpsGRS80 unitsm no_defs) county_m county.to_crs(ALBERS) county_m[area_km2] county_m.area / 1e6 # 交叉验证和粗略已知量级对比差一个数量级就是选错了投影 print(county_m[area_km2].describe())lat_1和lat_2是双标准纬线取 25°N 和 47°N 是中国全境通用的分带参数lon_0105是中央经线ellpsGRS80保证和 CGCS2000 椭球一致。这三个值不要随手改改了面积就不等了。/ 1e6是把平方米换成平方公里describe()出来的最大值应该在几千的量级如果出现几十万说明to_crs没生效或者原数据本来就是投影坐标又被二次投影了。3. 三条加载路径QGIS 体检、GeoPandas 批处理、PostGIS 建库同一批数据不同角色用不同工具。做一次性出图的用 QGIS 最快做批量统计的走 GeoPandas要多人共享或者数据量大到内存放不下的落 PostGIS。三条路径的坑各不相同分开说。3.1 QGIS 里先做一次完整性体检QGIS 的价值不在于画图而在于它能把几何问题可视化地暴露出来。把区划和路网一起拖进去之后先做三件事。第一件看几何有效性。菜单里跑一次有效性检查或者直接用字段计算器调is_valid()。多边形自相交在区划数据里非常常见尤其是沿河道、海岸线的边界源数据做过简化或者人工编辑后容易产生自相交和零面积碎片。这类几何在 QGIS 里显示完全正常但一旦做叠加分析就会抛拓扑异常。第二件看几何类型是否混杂。区划图层理论上应该是MultiPolygon但实际数据里经常多部件和单部件混存甚至混进GeometryCollection。路网图层要确认是LineString还是MultiLineString混合类型在部分空间库的索引构建阶段会直接报错。第三件看属性表的空值和重复主键。按主键字段排序肉眼扫一遍有没有空行用“按字段分组统计”看有没有一个代码出现两次。重复主键在后续 join 时会静默产生笛卡尔积式的行数膨胀几百行的表变几万行不查日志根本发现不了。QGIS 里还有个好用的动作把图层导出为 GeoPackage 时勾选“强制多部件”能一次性把几何类型统一掉省去后面在代码里explode的麻烦。3.2 GeoPandas 批量读入与字段规范化多个图层一起处理时别写一堆重复的read_file。用一个字典把路径和图层名管起来循环读入、循环体检代码短且不容易漏层。from pathlib import Path import geopandas as gpd src Path(./js_2022.gpkg) targets { province: province, city: city, county: county, road: road, highway: highway, } gdf {} for key, layer in targets.items(): g gpd.read_file(src, layerlayer) gdf[key] g # 一次性打印要素数、坐标系、几何类型种类异常层立刻能看出来 print(f{key:9s} 要素{len(g):6d} CRS{g.crs} 类型{sorted(g.geom_type.dropna().unique())})g.crs如果打印出None说明.prj在转换过程中丢了必须回到上一步补-a_srs。geom_type.dropna()是为了跳过空几何sorted()让输出可比一眼能看出某个层的类型集合是不是单元素。字段规范化单独抽一个函数把所有层的名称列和代码列统一成name、adcode后面做关联时就不用每处都判断字段名def normalize_admin(g, code_cands, name_cands, level): g g.copy() g[adcode] g[pick(g, code_cands)].astype(str).str.strip().str.zfill(6) g[name] g[pick(g, name_cands)].astype(str).str.strip() g[level] level # 只保留必要列避免不同层字段名冲突 return g[[adcode, name, level, geometry]] city normalize_admin(gdf[city], [ADCODE,CODE], [NAME,名称], city) county normalize_admin(gdf[county], [ADCODE,CODE], [NAME,名称], county)只保留必要列这一步看着浪费实际很关键。三级区划图层往往带了几十个统计字段字段名在 DBF 里已经被截断成NAME_1、NAME_12这类跨层合并时极易撞名或者被 pandas 自动加后缀排查起来非常费劲。统一裁剪到adcode、name、level、geometry四列后面所有操作都在这个干净结构上做。3.3 PostGIS 落库shp2pgsql 参数与空间索引数据量上到几十万条线、或者需要多人共享查询时PostGIS 是更合适的选择。导入用shp2pgsql参数选择上有几个关键点。# -s 指定源 SRID-W 指定源编码-I 建空间索引-D 用 dump 模式比默认的 copy 更快 shp2pgsql -s 4490 -W GBK -I -D ./js_2022/road.shp gis.road \ | psql -h localhost -p 5432 -U gis -d gisdb shp2pgsql -s 4490 -W GBK -I -D ./js_2022/county.shp gis.county \ | psql -h localhost -p 5432 -U gis -d gisdb-s 4490必须和源数据的真实基准一致写错了 PostGIS 会认为数据就在那个坐标系下之后做ST_Transform会算出完全错误的结果而且不会报错。-W GBK对应源文件的编码如果源文件带.cpg且声明为 UTF-8这里就改成 UTF-8不确定的话先用ogrinfo看一眼或者干脆先转成 GeoPackage 再导入把编码问题前置解决。-I建 GiST 空间索引没有索引的ST_Intersects在大表上会退化成全表扫描百万级要素能跑上几十分钟。-D用 PostgreSQL 的 dump 格式比默认的 insert 模式快一个数量级代价是不生成可读的 SQL 文本。落库之后立刻验证 SRID 和要素数-- 确认 SRID 与几何类型以及要素数是否和源文件对得上 SELECT ST_SRID(geom) AS srid, GeometryType(geom) AS gtype, COUNT(*) AS n FROM gis.county GROUP BY 1, 2; -- 补一个 GIST 索引如果 -I 没生效 CREATE INDEX IF NOT EXISTS idx_county_geom ON gis.county USING GIST (geom); CREATE INDEX IF NOT EXISTS idx_road_geom ON gis.road USING GIST (geom); -- 主键索引join 时用得上 CREATE INDEX IF NOT EXISTS idx_county_adcode ON gis.county (adcode);Group by 1, 2那个查询能同时暴露两个问题SRID 不是 4490 说明导入参数写错了gtype出现GEOMETRYCOLLECTION或者两种类型并存说明源数据几何混杂后面做叠加前必须先修。这两个问题在导入阶段解决成本远低于分析跑到一半再回头排查。4. 区划叠路网归属判定、密度计算与分级统计的三个口径数据加载好之后真正的工作量才刚开始。区划和路网叠加这块同一个问题至少有三种写法结果能差出百分之几十。差别不在代码水平在口径。4.1 线落哪个县空间连接的两种写法与差异最直觉的做法是空间连接把每条路和它相交的县连起来。GeoPandas 里一行sjoin搞定。但这里有个隐蔽的问题——一条横穿三个县的国道在sjoin之后会变成三行每行的几何都是那条完整国道长度字段也是全长。按县分组求和这条路的长度被完整计入了三个县总数虚高。正确的做法是先做几何切割把线按面边界打断再统计每段的长度import geopandas as gpd from shapely import make_valid ALBERS (projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 ellpsGRS80 unitsm no_defs) road_m gdf[road].to_crs(ALBERS) county_m gdf[county].to_crs(ALBERS) # 叠加前必修几何否则 overlay 大概率抛 TopologyException road_m[geometry] road_m.geometry.apply(make_valid) county_m[geometry] county_m.geometry.apply(make_valid) # howintersection 只保留有交集的部分线被面边界切断 clip gpd.overlay( road_m[[geometry]], county_m[[adcode, name, geometry]], howintersection, keep_geom_typeTrue, ) # 切割后重新算长度这一步是关键不能用切割前的长度 clip[seg_km] clip.length / 1000.0 stat clip.groupby([adcode, name], as_indexFalse).agg(road_km(seg_km, sum))make_valid是必须的overlay内部走的是 GEOS 的拓扑运算遇到自相交多边形会直接抛异常报错信息里带TopologyException和一对坐标那对坐标就是问题位置。keep_geom_typeTrue是为了过滤掉面线相交产生的点状碎片不加这个参数结果里会混进Point后面算长度全部是 0。clip.length必须在overlay之后再算因为overlay返回的是切割后的新几何长度字段不会自动更新沿用切割前的值等于白切。如果只是想快速看个大概sjoin也不是不能用但要清楚它统计的是“途经长度”而不是“境内长度”做汇报时必须说明口径否则和省里公布的里程数对不上解释起来很尴尬。4.2 路网密度面积口径、长度口径与单位换算路网密度最常见的定义是每平方公里国土面积上的公路里程单位 km/km²。这个指标看起来简单实际有三个口径要统一。面积口径。分母用的是行政区总面积还是扣除水域、山体后的可建设面积结果能差一倍以上。区划图层给的是行政辖区面积那就老老实实用这个在报告里注明“行政面积口径”别自己脑补扣除。长度口径。分子是全部路网还是只算等级公路差别巨大。乡道、村道在城市边缘能占一半以上的里程。分级统计时必须明确参与求和的道路等级。投影口径。这一条最容易被忽略。在 EPSG:4326 下算.length得到的是度数数值在 0.01 量级除以面积也是平方度之后看似无量纲实际上毫无意义。必须投影到米制之后再做除法。# 面积在 Albers 投影下计算单位换算为平方公里 county_m[area_km2] county_m.area / 1e6 stat stat.merge( county_m[[adcode, area_km2]], onadcode, howleft, validateone_to_one, # 加校验主键重复会立刻报错而不是静默膨胀 ) # km/km²若要换算成常见的 km/100km²乘 100 stat[density_km_per_km2] stat[road_km] / stat[area_km2] stat[density_km_per_100km2] stat[density_km_per_km2] * 100 print(stat[density_km_per_km2].describe())validateone_to_one是 merge 时最值得养成的习惯。缺了它主键有重复就会产生多对多笛卡尔积行数悄悄翻倍后面的均值、最大值全部失真而且不会报任何错。加上之后重复主键立刻抛异常问题在源头就被拦住。密度值describe()出来中位数落在 1 到 2 的区间比较合理如果出现 20 以上多半是面积单位没换算或者分母被覆盖成了小数值。4.3 按道路等级分层等级字段清洗与透视路网图层一般带一个等级字段取值可能是“高速、国道、省道、县道、乡道”这类中文也可能是G/S/X/Y这类单字母编码还可能是数字码。清洗的第一步不是映射是先看清实际有哪些值import pandas as pd road gdf[road].copy() grade_col pick(road, [GRADE, 等级, DLDJ, TYPE]) # 先看真实取值分布不要凭经验假设 print(road[grade_col].astype(str).str.strip().value_counts(dropnaFalse))打印出来经常能看到这种情况同一个等级有“高速公路”“高速”“G”三种写法还有带前后空格的、带全角字符的、以及None和字符串NULL混在一起。这时候做映射要先把它们收敛# 先做字符串规范化去空格、去全角、统一大小写 road[grade_raw] ( road[grade_col] .astype(str) .str.replace(\u3000, , regexFalse) # 全角空格 .str.strip() .str.upper() ) # 映射表按实际取值补全这里只列出常见情形 grade_map { 高速公路: expressway, 高速: expressway, G: expressway, 国道: national, S: provincial, 省道: provincial, 县道: county, 乡道: township, } road[grade] road[grade_raw].map(grade_map).fillna(other).map()之后fillna(other)把没匹配上的全部归到“其他”而不是丢掉。丢掉会让总里程对不上归到“其他”至少总量守恒后续还能回过头看有多少没识别出来。如果“其他”占比超过一成说明映射表不全回去把value_counts的输出逐条过一遍。分层统计用透视表最省事行是县、列是等级一次出全clip[grade] clip.index.map(road[grade]) # 按索引对齐需先 reset_index 保证唯一 pivot ( clip.pivot_table(indexname, columnsgrade, valuesseg_km, aggfuncsum) .fillna(0) .round(2) ) # 每行补一个合计方便和总量核对 pivot[total] pivot.sum(axis1) print(pivot.sort_values(total, ascendingFalse).head(10))pivot_table的aggfuncsum必须显式写默认是求均值不写会得到一堆莫名其妙的平均值。fillna(0)是因为某些县可能完全没有高速透视后是NaN不填零的话后续做加减会污染整列。补total列之后和 4.1 里算出的road_km对一下两个数应该完全一致不一致说明中间的索引对齐或过滤步骤出了问题。4.4 结果出图与导出统计结果最终要落到图上。出图这块不用写得复杂关键是分级方式要说清楚。等间距分级在密度这种右偏分布上会让绝大多数区域挤在最低档看起来一片同色分位数分级视觉效果好但每档的边界值不规整汇报时不好解释。我一般用自然断点同时把断点值导出来存档。# 自然断点需要 mapclassify先装pip install mapclassify stat.plot( columndensity_km_per_km2, schemenatural_breaks, k5, legendTrue, edgecolorwhite, linewidth0.3, figsize(9, 7), ) # 导表时把分级边界一并写出方便复现 stat.drop(columnsgeometry).to_csv(./density_by_county.csv, indexFalse, encodingutf-8-sig)k5是经验值超过 7 档图例就开始挤。encodingutf-8-sig是为了 Excel 直接打开不乱码这个细节在做交付时很省事。导出 CSV 前drop(columnsgeometry)几何列没法序列化成文本留着会报错或者变成一坨 WKT 字符串。5. 避坑排查五个把好数据用成废数据的现场记录这一节全是实际踩出来的。每一条都按现象、原因、解决三段写遇到同类问题时可以直接对号入座。5.1 五条踩坑记录第一条中文属性全变问号。现象图层能打开图形正常但属性表里所有中文变成???或者方块用ogrinfo看倒是有内容。原因.cpg文件缺失或者内容声明为UTF-8而实际字节是 GBK。Shapefile 本身不带编码信息全靠.cpg或者调用方指定两者对不上就乱码。解决先用ogrinfo -so看要素能不能读出能读说明只是显示编码问题。用ogr2ogr转一次 GeoPackage 并显式指定-lco ENCODINGUTF-8如果转出来还是乱码说明源数据确实是 GBK那就把读入时的encoding参数显式写成gb18030比gbk覆盖字符更全不要依赖自动探测。第二条面积算出来差了几万倍。现象gdf.area求出来的最大值是零点几或者几百万。原因在 EPSG:4326 地理坐标下直接算面积单位是“平方度”一个度在赤道约 111 公里平方之后量级完全对不上。反过来如果数据本来就是投影坐标又被to_crs(4326)转了一次再算面积同样离谱。解决算之前先print(gdf.crs)确认基准然后无条件to_crs到米制投影。全省级统计用 Albers单市精细量算用本地 3 度带。算完之后用describe()交叉验证量级县级面积正常应该在几百到几千平方公里。第三条全省用单带投影东西两头数据变形。现象用中央经线 120°E 的 3 度带投影全省做出来的密度图西侧几个县明显偏低和常识不符。原因3 度带的带宽只有 3 个经度覆盖范围约 1.5 个经度江苏东西跨度远超这个范围离中央经线越远长度变形越大。长度变形直接进入密度计算的分子。解决全省级统计一律换 Albers 等积投影面积严格保持长度变形在可接受范围。只有在做单个地级市内部、跨度小于一个带的范围时才用 3 度带高斯克吕格。这条不是理论洁癖是实际做出来的图能不能看的问题。第四条overlay 抛 TopologyException。现象gpd.overlay跑到一半报错信息里带TopologyException和一对坐标。原因参与运算的多边形存在自相交、悬挂边、零面积碎片或者几何类型混杂了GeometryCollection。区划数据经过多次编辑和格式转换后这类问题非常普遍QGIS 里看着完全正常。解决叠加前对两侧几何都跑make_valid并且用explode(index_partsFalse)把多部件拆成单部件。如果还有个别要素过不去用buffer(0)兜底但要说清楚buffer(0)会轻微改变边界形状对精度要求高的场景应该定位到具体要素手工修而不是全表一刀切。第五条字段名被截断、长文本丢失。现象导出的.dbf里字段名变成NAME_1、ADCODE_2超过 254 字节的备注列后半截没了。原因Shapefile 的 DBF 规范硬限制字段名最多 10 个字节UTF-8 下中文只够 3 个字字符字段最长 254 字节。这是格式本身的约束不是软件 bug。解决分析阶段一律用 GeoPackage 或者 PostGIS不用.dbf。确实需要交付 SHP 时把长字段拆成多个短字段或者另附一份 CSV 用主键关联。别试图在 SHP 里塞长文本那是跟格式规范较劲。5.2 一个批量自检脚本上面五条里前四条都可以用脚本自动查出来。每次拿到新数据先跑一遍比事后一个个排查省时间import geopandas as gpd from shapely.validation import explain_validity def audit(path, layer, name_col_candidates(NAME, 名称)): g gpd.read_file(path, layerlayer) print(f {layer} ) print(f要素数: {len(g)}) print(f坐标系: {g.crs}) print(f几何类型: {sorted(g.geom_type.dropna().unique())}) print(f空几何: {g.geometry.isna().sum()}) # 几何有效性返回第一个无效要素的原因 invalid g[~g.geometry.is_valid] print(f无效几何: {len(invalid)}) if len(invalid): print(首个无效原因:, explain_validity(invalid.geometry.iloc[0])) # 主键检查空值 重复 for c in name_col_candidates: if c in g.columns: s g[c].astype(str).str.strip() print(f字段 {c}: 空值{s.isin([, nan, None]).sum()} 重复{s.duplicated().sum()}) break audit(./js_2022.gpkg, county) audit(./js_2022.gpkg, road)explain_validity会返回一句人话比如“Self-intersection at or near point (…)”或者“Ring Self-intersection”比异常的堆栈信息好读得多。主键检查那段只查第一个命中的候选字段因为三级区划的名称字段名可能不一样但每个图层内部认一个就够。isin([, nan, None])是为了兜住astype(str)把缺失值转成字符串nan的情况这个坑在 pandas 里非常常见。6. 进阶用法把一次性 SHP 固化成可复用的分析底图同一个项目里底图往往要被反复用几十次。每次都从 SHP 读、每次转投影、每次修几何累积起来的时间很可观而且不同人不同时间跑的版本还可能不一致。更稳的做法是第一次拿到数据就固化成一份“分析底图包”统一坐标系、预计算面积长度、统一主键、附一份元数据。import json import geopandas as gpd from shapely import make_valid ALBERS (projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 ellpsGRS80 unitsm no_defs) def bake(src_path, layer, out_path, code_cands, name_cands, geom_fixTrue): g gpd.read_file(src_path, layerlayer) g[adcode] g[pick(g, code_cands)].astype(str).str.strip().str.zfill(6) g[name] g[pick(g, name_cands)].astype(str).str.strip() g g[[adcode, name, geometry]] # 统一到等积投影预计算几何量后续不再重复转换 g g.to_crs(ALBERS) if geom_fix: g[geometry] g.geometry.apply(make_valid) g[area_km2] g.geometry.apply(lambda x: x.area / 1e6) g[len_km] g.geometry.apply(lambda x: x.length / 1000.0) # GPKG 支持同库多图层一个文件装下整套底图 g.to_file(out_path, layerlayer, driverGPKG) return {layer: layer, count: len(g), crs: ALBERS} meta {} for layer, codes, names in [ (county, [ADCODE, CODE], [NAME, 名称]), (road, [ROADID, ID], [NAME, 名称]), ]: meta[layer] bake(./js_2022.gpkg, layer, ./base.gpkg, codes, names) meta[source] 2022.07 江苏省行政区划与路网 SHP meta[note] 面积与长度在 Albers 等积投影下预计算单位分别为 km2 和 km with open(./base_meta.json, w, encodingutf-8) as f: json.dump(meta, f, ensure_asciiFalse, indent2)g.to_file追加同一路径时驱动必须显式写driverGPKG不写的话fiona可能按扩展名推断失败或者覆盖已有图层。meta字典里那份 JSON 是给三个月后的自己看的特别是note字段把投影参数和单位写明白比在聊天记录里翻找当初怎么定的口径要靠谱得多。还有一个细节值得单独提路网图层在bake阶段就把len_km算好后面做密度时直接groupby求和就行不用每次重新投影求长度。代价是这份预计算的长度是“线本身长度”跨县的那条线在按县统计时仍然需要切割重算这一点不能省。预计算的价值在于快速出总里程、做数据校验而不是替代切片。从那以后我每拿到一批 SHP第一件事都是跑一遍audit加ogrinfo -so把编码、坐标系、几何有效性、主键重复这四项过一遍确认干净了再动手。这个习惯至少省掉了我三次“图做完了发现数据是错”的返工。希望帮到你。本文还有配套的精品资源点击获取