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.cpg

ogrinfo -so输出的Extent那两行要重点看：数值落在 -180~180 和 -90~90 之间，说明是地理坐标（单位是度）；数值是六位、七位数，说明已经是投影坐标（单位是米）。这一步判断错了，后面所有面积统计都会差出好几个数量级。gdalsrsinfo会打印 EPSG 编号或 WKT，国内 2022 年前后的区划数据常见做法是 CGCS2000 地理坐标，EPSG 为 4490，但不要默认，一定要读出来确认。

如果.cpg缺失或者内容写的是 UTF-8 而实际数据是 GBK，属性表里的中文会变成问号或方块。这种情况不要在原文件上反复折腾，直接转成 GeoPackage 一次性解决编码和字段长度的双重限制：

# -lco ENCODING=UTF-8 强制以 UTF-8 落库，-nln 指定新图层名 ogr2ogr -f GPKG ./js_2022.gpkg ./js_2022/county.shp \ -nln county -lco ENCODING=UTF-8 -a_srs EPSG:4490

-a_srs是给那些.prj丢失或写错的文件强行指定坐标系，用之前必须确认数据的真实基准，否则等于给错误答案盖了个正确的章。GeoPackage 支持任意长度字段名和 UTF-8，DBF 那两个坑（字段名最多 10 字节、字符字段最多 254 字节）在 GPKG 里都不存在，长文本备注类字段一定要走这一步。

2.2 属性字段怎么读：区划代码是唯一可靠的主键

三级区划图层最容易踩的坑，是把“名称”当关联键。名称在不同版本里可能有简称、全称、带“市”不带“市”的差异，甚至同一个名字在不同层级重复出现。真正稳的是行政区划代码：省级 2 位补 6 位、地市级 4 位补 6 位、县级 6 位完整码。代码本身有层级含义，前 2 位是省，前 4 位是市，做上下级聚合时直接截字符串就行，不用做任何名称匹配。

典型字段结构大致如下，字段名以实际.dbf为准，不同来源命名习惯差别很大：

表格里最后两行是最容易被误用的。数据源方给的面积字段，你不知道它是在什么投影下算的、有没有扣掉水域、是不是投影变形区的结果。自己用等积投影重算一遍，几十行代码的事，比事后解释数据对不上要省心得多。

读入后第一件事是规范化主键：

import geopandas as gpd county = gpd.read_file("./js_2022.gpkg", layer="county", encoding="utf-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 里那一行决定后面所有数字对不对

国内区划和路网数据常见的坐标系有三类，用错一个，结果就废：

第三行是全省统计的正解。江苏东西跨度大，用单一 3 度带（比如中央经线 117°E 或 120°E）覆盖全省，边缘地带的长度和面积都会有可观偏差；用 Albers 等积投影，面积是严格保持的，长度虽然有变形但全省范围内可以接受，做密度类指标完全够用。

# 中国标准 Albers 等积投影，GRS80 椭球对应 CGCS2000 基准 ALBERS = ("+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 " "+x_0=0 +y_0=0 +ellps=GRS80 +units=m +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_0=105是中央经线，ellps=GRS80保证和 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, layer=layer) 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 = ("+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 " "+x_0=0 +y_0=0 +ellps=GRS80 +units=m +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) # how="intersection" 只保留有交集的部分，线被面边界切断 clip = gpd.overlay( road_m[["geometry"]], county_m[["adcode", "name", "geometry"]], how="intersection", keep_geom_type=True, ) # 切割后重新算长度，这一步是关键，不能用切割前的长度 clip["seg_km"] = clip.length / 1000.0 stat = clip.groupby(["adcode", "name"], as_index=False).agg(road_km=("seg_km", "sum"))

make_valid是必须的，overlay内部走的是 GEOS 的拓扑运算，遇到自相交多边形会直接抛异常，报错信息里带TopologyException和一对坐标，那对坐标就是问题位置。keep_geom_type=True是为了过滤掉面线相交产生的点状碎片，不加这个参数结果里会混进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"]], on="adcode", how="left", validate="one_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())

validate="one_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(dropna=False))

打印出来经常能看到这种情况：同一个等级有“高速公路”“高速”“G”三种写法，还有带前后空格的、带全角字符的、以及None和字符串"NULL"混在一起。这时候做映射要先把它们收敛：

# 先做字符串规范化：去空格、去全角、统一大小写 road["grade_raw"] = ( road[grade_col] .astype(str) .str.replace("\u3000", "", regex=False) # 全角空格 .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(index="name", columns="grade", values="seg_km", aggfunc="sum") .fillna(0) .round(2) ) # 每行补一个合计，方便和总量核对 pivot["total"] = pivot.sum(axis=1) print(pivot.sort_values("total", ascending=False).head(10))

pivot_table的aggfunc="sum"必须显式写，默认是求均值，不写会得到一堆莫名其妙的平均值。fillna(0)是因为某些县可能完全没有高速，透视后是NaN，不填零的话后续做加减会污染整列。补total列之后，和 4.1 里算出的road_km对一下，两个数应该完全一致，不一致说明中间的索引对齐或过滤步骤出了问题。

4.4 结果出图与导出

统计结果最终要落到图上。出图这块不用写得复杂，关键是分级方式要说清楚。等间距分级在密度这种右偏分布上会让绝大多数区域挤在最低档，看起来一片同色；分位数分级视觉效果好但每档的边界值不规整，汇报时不好解释。我一般用自然断点，同时把断点值导出来存档。

# 自然断点需要 mapclassify，先装：pip install mapclassify stat.plot( column="density_km_per_km2", scheme="natural_breaks", k=5, legend=True, edgecolor="white", linewidth=0.3, figsize=(9, 7), ) # 导表时把分级边界一并写出，方便复现 stat.drop(columns="geometry").to_csv("./density_by_county.csv", index=False, encoding="utf-8-sig")

k=5是经验值，超过 7 档图例就开始挤。encoding="utf-8-sig"是为了 Excel 直接打开不乱码，这个细节在做交付时很省事。导出 CSV 前drop(columns="geometry")，几何列没法序列化成文本，留着会报错或者变成一坨 WKT 字符串。

5. 避坑排查：五个把好数据用成废数据的现场记录

这一节全是实际踩出来的。每一条都按现象、原因、解决三段写，遇到同类问题时可以直接对号入座。

5.1 五条踩坑记录

第一条，中文属性全变问号。

现象：图层能打开，图形正常，但属性表里所有中文变成???或者方块，用ogrinfo看倒是有内容。

原因：.cpg文件缺失，或者内容声明为UTF-8而实际字节是 GBK。Shapefile 本身不带编码信息，全靠.cpg或者调用方指定，两者对不上就乱码。

解决：先用ogrinfo -so看要素能不能读出，能读说明只是显示编码问题。用ogr2ogr转一次 GeoPackage 并显式指定-lco ENCODING=UTF-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_parts=False)把多部件拆成单部件。如果还有个别要素过不去，用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, layer=layer) 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 = ("+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 " "+x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs") def bake(src_path, layer, out_path, code_cands, name_cands, geom_fix=True): g = gpd.read_file(src_path, layer=layer) 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, layer=layer, driver="GPKG") 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", encoding="utf-8") as f: json.dump(meta, f, ensure_ascii=False, indent=2)

g.to_file追加同一路径时，驱动必须显式写driver="GPKG"，不写的话fiona可能按扩展名推断失败或者覆盖已有图层。meta字典里那份 JSON 是给三个月后的自己看的，特别是note字段，把投影参数和单位写明白，比在聊天记录里翻找当初怎么定的口径要靠谱得多。

还有一个细节值得单独提：路网图层在bake阶段就把len_km算好，后面做密度时直接groupby求和就行，不用每次重新投影求长度。代价是这份预计算的长度是“线本身长度”，跨县的那条线在按县统计时仍然需要切割重算，这一点不能省。预计算的价值在于快速出总里程、做数据校验，而不是替代切片。

从那以后我每拿到一批 SHP，第一件事都是跑一遍audit加ogrinfo -so，把编码、坐标系、几何有效性、主键重复这四项过一遍，确认干净了再动手。这个习惯至少省掉了我三次“图做完了发现数据是错”的返工。希望帮到你。

本文还有配套的精品资源，点击获取