行政区划shp数据这几个字看着平平无奇，但凡做过区域类项目的人都清楚，它几乎是最容易拿到、也最容易用歪的一份基础数据。下载下来解压，一堆同名不同后缀的文件摊在文件夹里，拖进桌面软件能正常显示，颜色一配、图例一挂，看起来就"能用"了。可一旦你要拿它去挂业务数据、做区域统计、做看板下钻、做面积测算，前面偷的每一步懒都会在后面变成对不上账的麻烦。我这些年在区域分析、门店选址、物流分单、看板地图这几类活儿里，反复跟这份数据打交道，踩过的坑从编码乱码到坐标系误判，从共享边不重合到区划调整导致的代码失效，几乎把能踩的都踩了一遍。这篇就把这套经验摊开讲：从拿到文件那一刻该做什么体检，到层级代码怎么设计主键，到投影怎么选、面积怎么算，再到空间连接和版本回溯的完整链路。不管你是刚接触地理数据的开发，还是被区域统计折磨过一轮的分析师，下面这些内容应该都能直接拿去用。

1. 拿到行政区划shp先别急着加载：先想清楚你要用它回答什么问题

1.1 三种用法对应三套完全不同的处理路径

很多人把行政区划shp数据当成一张"底图"，这个定位本身就限制了它的价值。它其实是一套带层级关系的空间主数据，用法的差异会直接决定后面所有技术选择。

制图展示型只关心形状好看不好看。这种场景下，坐标系随便，Web墨卡托就行，拓扑有点小缝隙也看不出来，字段乱码更是无所谓，反正只画不标。这是成本最低的用法，也是唯一的"拖进去就能用"的用法。

统计挂接型才是主流需求。比如把订单表按区域汇总、把门店地址归到某个区县、把用户数按行政区聚合。这种场景的核心矛盾是主键唯一性和名称标准化——区划代码必须干净、不能重复、不能有前后空格，名称必须能对上业务表里的写法。这时候坐标系反而次要，因为你多半是在做属性关联而不是几何运算。

分析计算型最苛刻。面积、人口密度、辐射范围、缓冲区覆盖、区域间距离，这些都要算数值。一旦涉及数值，坐标系选择就成了生死线，投影错了算出来的面积能差出一个量级；拓扑不干净，叠加分析会产生一堆碎渣多边形。

我的习惯是：在打开文件之前，先用一句话写清楚这个项目最终要输出什么。是"一张能看的地图"，还是"一张区域汇总表"，还是"一组空间测算指标"。写不出来就先别动数据，因为你根本不知道自己该清洗到什么程度。

1.2 把所有层级合并成一个图层，是最常见的第一个错误

新手最容易干的一件事，是把省级、市级、县级三个shp合并成一个大图层，然后一次性渲染，觉得这样"图上有省界也有县界，很完整"。

问题在于，行政区划是嵌套关系。一个省的面，几何上就包含了它下属所有市的面。你把省和市叠在同一个图层里，同一个空间位置会被重复绘制两到三次。视觉效果上你可能觉得"边界线变粗了"，但实际后果是：任何基于这个图层的统计都会重复计数。

举个具体场景。你想统计"全国有多少个区域单元"，合并后的图层会给你省数+市数+县数，这个数字没有任何业务含义。你想算总面积，三个层级的面叠在一起，总面积会是实际面积的两到三倍。

正确的做法是分层存放、按需调用。省级一份、市级一份、县级一份，各自独立成层，需要哪个层级就加载哪个。确实需要在一个视图里同时展示两级边界时，用样式区分（省界粗线、县界细线），而不是把数据物理合并。

1.3 先确认这份数据的时间截面

行政区划不是静态的。撤县设区、合并、新设、更名，每年都在发生。一份数据里如果混了不同年份的调整结果，会出现非常诡异的状况：某个代码在省级图层里存在，在县级图层里却找不到对应的父级。

我在一个项目里遇到过：市级图层里有一个单元，代码前缀在省级图层里查不到。排查了半小时才反应过来，这是两个不同批次的数据拼在一起的，中间正好赶上一轮调整。

所以拿到数据第一件事，翻一下有没有元数据、有没有说明文档、有没有标注数据截止时间。没有的话，就自己记一个"入库时间"，后面所有对不上的问题都从这里开始回溯。

2. Shapefile四件套与DBF字段限制，决定了你逃不掉的一轮体检

2.1 .shp、.shx、.dbf、.prj 各管什么，缺一个会怎样

Shapefile从来不是"一个文件"，它是至少四个文件组成的集合。这个设计来自很早的年代，一直沿用至今，理解它各自的职责能省掉很多莫名其妙的故障。

.prj缺失是最阴险的一种。因为软件不会报错，它会默默按某种默认坐标系渲染，图上看起来"形状是对的"，你根本发现不了。等到你把这份数据和另一份带正确坐标系的数据叠在一起，才发现整体错位了几百公里。

我见过的补救办法是：从同一批数据的其他图层里复制一份.prj过来，前提是你确认它们坐标系一致。如果连参照都没有，只能靠已知地标的实际经纬度去反推，成本很高。

还有个小细节：文件名千万别带中文和空格。某些老版本驱动在路径解析上会出问题，尤其是路径里带空格的时候，读出来是空表但不是报错，排查起来很费时间。养成习惯，解压后先重命名成英文短名。

2.2 DBF的字段名10字符上限与中文编码

DBF这个格式的字段名上限是10个字符。这意味着district_name会被截断成district_n，administrative_code会变成administ。如果你按原字段名去取值，拿到的是空。

编码问题更麻烦。DBF里的中文有两种常见存法：GBK和UTF-8。读的时候编码猜错，字段值就是一堆问号或者方块。判断方法很简单——先用UTF-8读一遍，如果中文全是乱码，换GBK再试一次。这中间还有一种情况是"部分正常部分乱码"，说明数据源头就混了编码，只能逐行修。

字段值的长度限制是254字节。UTF-8下一个汉字占3字节，理论上能放80多个汉字，一般行政区名称够用。但如果你的数据里塞了备注、别名、历史沿革说明，超长部分会被静默截断，不报错。等你发现某个名称只剩半截的时候，已经找不到原文了。

提示：DBF的文本截断是静默发生的，写入时不会给你任何警告。如果你的属性里有可能超过200字节的文本，提前拆到独立的关联表里，别硬塞进shp。

2.3 一份能直接跑的体检脚本

下面这段是我每次都用的体检流程，输出结果贴在项目文档里，后面出问题先翻它。

import geopandas as gpd import fiona import pandas as pd def inspect_shp(path, encoding='utf-8'): # 第一层：用fiona看原始元信息，不经过pandas转换 with fiona.open(path) as src: print("驱动:", src.driver) print("坐标系:", src.crs) print("几何类型:", src.schema['geometry']) print("字段清单:", list(src.schema['properties'].keys())) print("要素总数:", len(src)) # 第二层：读成GeoDataFrame做质量检查 gdf = gpd.read_file(path, encoding=encoding) print("\n各字段类型:\n", gdf.dtypes) print("空几何数量:", gdf.geometry.isna().sum()) print("无效几何数量:", (~gdf.geometry.is_valid).sum()) print("多部件要素数量:", (gdf.geometry.type == 'MultiPolygon').sum()) # 代码字段是否唯一，这是挂接成败的关键 code_col = [c for c in gdf.columns if 'code' in c.lower()] if code_col: col = code_col[0] dup = gdf[col].duplicated().sum() print(f"字段 {col} 重复值数量: {dup}") return gdf gdf = inspect_shp("admin_county.shp", encoding='utf-8')

几个判断标准，供参考：

• 无效几何数量不为0：必须修，不然后面所有几何运算都可能抛异常或者返回空结果。
• 代码字段有重复：要么是数据本身有问题，要么是你拿到了多层级合并的版本。
• 多部件要素比例很高：说明这个区域包含大量飞地或岛屿，后面算中心点要格外小心。
• 坐标系为空：.prj丢了，先解决这个再谈别的。

3. 区划代码与行政层级：让"省—市—县—乡"变成可计算的键

3.1 代码长度本身就编码了层级，别把它当普通数字

行政区划代码是一套定长数字编码，长度直接对应层级。省级是6位（后4位为0），市级6位（后2位为0），县级6位，乡镇级9位，再往下到村级是12位。

这个特性非常好用，因为它让你可以纯靠字符串截取就构建出完整的父子关系，不需要额外的层级字段。但前提是——必须按字符串处理。

我见过太多人因为把代码列读成了整数，导致前导零丢失（比如以0开头的代码变成了少一位的数），然后所有的截取逻辑全部错位。用pd.read_csv的时候记得加dtype={'code': str}，用gpd.read_file读shp时DBF里的字段一般是文本类型，但也建议显式转换一遍。

gdf['code'] = gdf['code'].astype(str).str.strip().str.zfill(6) # 派生各级父级键 gdf['prov_key'] = gdf['code'].str[:2] + '0000' gdf['city_key'] = gdf['code'].str[:4] + '00' gdf['county_key'] = gdf['code'].str[:6]

这段代码的好处是：不管你手上是县级图层还是乡镇级图层，都能一次性算出它归属的省、市、县键，后面做任何层级的汇总都不用再回头找上级图层。

3.2 名称对不上，是挂接失败的头号原因

绝大多数区域统计失败，都不是几何问题，而是名称匹配失败。常见的不一致形态：

• 空格差异：业务表里是"某某 区"，数据里是"某某区"，肉眼看不出来，字符串比较直接失败。
• 全角半角混用：括号、数字、标点，全角半角和半角全角在视觉上几乎一致。
• 后缀差异：有的业务表写"某某市"，有的写"某某"，有的写"某某地区"。
• 曾用名与新名：区划调整后，业务系统里的历史数据还是老名字。
• 行政区与功能区混淆：某些经济功能区、开发区在区划数据里根本没有独立代码，但在业务表里是一级维度。

我的处理原则是：能靠代码挂接就绝不靠名称。名称只作为展示字段和兜底匹配手段。如果业务表确实只有名称没有代码，那就必须建一张名称映射表，把业务侧的写法一对多映射到标准代码上，并且这张表要长期维护——每次遇到新的写法就补一条，而不是每次都在脚本里加if。

# 名称标准化：去空格、统一全角半角、统一括号 import re def normalize_name(s): if pd.isna(s): return s s = str(s).strip() s = s.replace('（', '(').replace('）', ')') s = re.sub(r'\s+', '', s) return s gdf['name_norm'] = gdf['name'].apply(normalize_name) biz['name_norm'] = biz['region_name'].apply(normalize_name)

标准化之后再匹配，命中率能提升一大截。剩下没命中的，逐条人工确认，补进映射表。

3.3 层级过滤与父子关系校验

构建完层级键之后，一定要做一次父子存在性校验：每个县级单元的city_key，在市级图层里是否都能找到？每个市级单元的prov_key，在省级图层里是否都存在？

prov_codes = set(prov_gdf['code']) city_codes = set(city_gdf['code']) orphan_city = city_gdf[~city_gdf['prov_key'].isin(prov_codes)] print("找不到父级的市级单元数:", len(orphan_city))

这个检查能抓出两类问题：一是数据版本不一致（前面提到过），二是代码本身有错误（比如位数不对、被截断）。孤儿单元如果不处理，后面做层级下钻的时候会出现"点进去是空的"这种诡异体验。

还有一种特殊情况要注意：直辖市和特别处理的层级。这类区域的市级和省级实际上是同一套几何，如果按常规的父子关系去做分层渲染或者面积汇总，会重复计算。处理方式是在数据里加一个level_type标记字段，把这类情况显式标注出来，汇总时排除掉。

4. 坐标系与投影：同一块面算出两个面积的根因在哪

4.1 经纬度单位是"度"，不是"米"

这是地理数据里最经典的认知陷阱。原始shp数据绝大多数用的是地理坐标系（单位是度），你在这种数据上直接调.area，得到的数字是"平方度"。

平方度是什么概念？在赤道附近，1度纬度大约111公里，1度经度也大约111公里；但到了高纬度，1度经度对应的实际距离会急剧缩小。所以同样大小的"1平方度"，在不同纬度代表的实际面积完全不同。用它做任何比较、汇总、排序，结论都是错的。

更迷惑的是，这个错误不报错、不警告，数字照常输出，只有当你拿它去和真实面积对照，才发现差了七八个数量级。

# 错误示范：在地理坐标系上直接算面积 gdf['area_wrong'] = gdf.geometry.area # 单位是平方度，无业务意义 # 正确做法：先转到等积投影 gdf_proj = gdf.to_crs("EPSG:XXXX") # 换成适合目标区域的等积投影 gdf_proj['area_m2'] = gdf_proj.geometry.area gdf_proj['area_km2'] = gdf_proj['area_m2'] / 1e6

4.2 等积、等距、等角，选错了就是白算

投影不是"随便挑一个能用的"，不同投影保留的特性完全不同。做区域统计，你需要的是等积投影；做距离测算，你需要的是等距投影；做导航和切片底图，等角投影才合适。

Web墨卡托是面积计算的重灾区。它为了保持形状和角度，把高纬度区域在纵向拉伸得越来越厉害。在它的坐标系下算面积，纬度越高误差越离谱，越靠近极地偏差越夸张。而偏偏很多在线地图服务和前端组件默认就用这个投影，好多人直接在上面算面积，结果当然是错的。

跨越大范围区域做等积投影时，参数怎么定？常见做法是：标准纬线取目标区域南北跨度的约六分之一和六分之五位置，中央经线取区域中轴经度。这样能让整个区域内的面积变形控制在很小的范围内。

4.3 面积计算的完整写法与误差校验

import geopandas as gpd from pyproj import CRS # 假设目标区域跨纬度 20~50 度，中轴经度 105 度 albers = CRS.from_proj4( "+proj=aea +lat_1=25 +lat_2=45 +lat_0=0 " "+lon_0=105 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" ) gdf_aea = gdf.to_crs(albers) gdf_aea['area_km2'] = gdf_aea.geometry.area / 1e6 print(gdf_aea['area_km2'].describe())

算完之后一定要做交叉校验，否则你不知道自己算对没算对。校验思路有几个：

一是总量对照。把所有单元面积加起来，看是否接近你已知的区域总面积。误差在百分之几以内算正常（投影变形+边界简化），差一个量级就一定是坐标系搞错了。

二是同层级横向对比。同类单元的面积的量级应该在一个合理区间内。如果某个县级单元算出来的面积和省级差不多，那要么是几何有重复部件，要么是多部件面被重复计入。

三是换一个投影再算一遍。用两个不同的等积投影各算一次，结果应该非常接近。如果两个投影算出来差很多，说明你选的投影在这个区域变形过大，或者数据有拓扑问题。

注意：面积误差在千分之几到百分之一这个范围内是正常的，属于投影变形和数据精度的合理范围。但如果两块相邻区域的面积之和明显大于它们合并后的面积，那就不是投影问题了，是几何重叠。

5. 拓扑与边界：共享边、飞地与多部件面引发的连锁反应

5.1 相邻面共享边不重合，会产生什么后果

理想状态下，两个相邻行政区的公共边应该是完全重合的一条线。但实际数据里，因为采集精度、简化算法、不同批次拼接等原因，这条边经常是"两条几乎重合但不完全重合"的线，中间留下极细的缝隙。

这在视觉上完全看不出来，但在算法层面是灾难：

• 空间叠加分析会出现大量细长的碎片多边形，属性为空，统计时要么被丢弃要么被错误计入。
• 相交运算返回的面积之和小于各自面积的总和，对不上账。
• 缓冲区分析会在缝隙处产生奇怪的断裂。

排查方法不难：把所有面做一次unary_union（合并），再和原始面的总面积对比，差值就是重叠或缝隙的规模。也可以用union_all()之后看结果里的多边形数量，如果碎片数量远大于原始要素数量，说明拓扑问题严重。

修复思路有两种。一种是用拓扑规则重建，让相邻边吸附对齐；另一种更实用——如果只是做统计，不做精确的边界运算，那就用合并后重切的方式：先把同层级的所有面合并成一个大面，再按代码重新切分。这样能保证没有缝隙和重叠，代价是边界会有细微变化。

5.2 多部件面：一个行政区为什么不是"一块"

行政区很少是一个规整的多边形。它可能包含海上的岛屿、被河流切开的片区、被其他区包围的飞地。这些都会让几何类型从Polygon变成MultiPolygon。

多部件带来的第一个坑是中心点计算。centroid在MultiPolygon上的行为是计算所有部件的几何中心，对于形状不规则的区域，这个点很容易落在区域外部。如果这个区域是环形的，中心点几乎必然落在中间的空洞里；如果是狭长形带飞地，中心点可能落在两个部件之间的空隙上。

这会导致什么？你把中心点当成"区域代表点"去做空间连接，结果这个点归属了隔壁区。

正确的做法是用representative_point()，它保证返回的面内一点：

# 不推荐：可能落在面外 gdf['pt_wrong'] = gdf.geometry.centroid # 推荐：保证点在面内 gdf['pt_ok'] = gdf.geometry.representative_point()

第二个坑是面积和面积权重。多部件面的面积是各部件之和，这个没问题。但如果你要做"部件级"分析，就得先explode()拆开。拆开之后每行是一个独立部件，此时要注意：部件的面积和原要素的代码不再是唯一对应关系，一个代码会对应多行，后续做聚合时要重新groupby回去。

parts = gdf.explode(index_parts=False).reset_index(drop=True) parts['part_area'] = parts.geometry.area part_summary = parts.groupby('code')['part_area'].agg(['sum', 'count'])

5.3 无效几何的修复与取舍

无效几何的典型形态：自相交（边界线自己穿过自己）、环方向错误（外环写成顺时针）、重复顶点、零面积多边形。这些在shapely里会被is_valid标出来，但不同的修复方式结果不同。

from shapely.validation import make_valid from shapely import set_precision # 方式一：make_valid，返回可能是GeometryCollection，需要提取面 gdf['geometry'] = gdf.geometry.apply(make_valid) # 方式二：buffer(0)，简单粗暴，但对自相交多边形效果不稳定 gdf['geometry'] = gdf.geometry.buffer(0) # 方式三：精度归一，处理近乎重合的顶点 gdf['geometry'] = gdf.geometry.apply(lambda g: set_precision(g, 0.000001))

三种方式的取舍：make_valid最稳妥，是首选；buffer(0)快但有风险，因为它会把一些本不应该修的形状也改掉；set_precision适合那种"坐标精度过高导致浮点比较失败"的场景，做完全局精度归一之后，很多共享边的匹配问题会自动解决。

我个人的顺序是：先set_precision统一精度，再对剩下的无效几何用make_valid。两步之后还有问题的，单独拎出来人工看，因为那说明数据源头就有问题。

6. 面转点与空间连接：业务数据挂到行政区上的完整链路

6.1 三种挂接方式怎么选

把业务数据（通常是带经纬度的点、或者带地址的文本）关联到行政区上，有三条路：

实际项目里，优先走代码关联，其次空间连接，最后名称兜底。三条路组合起来用，覆盖率能到接近百分之百。

空间连接的核心是判断"点落在哪个面里"。这里有个细节很容易被忽略：within和intersects的区别。

• within要求点严格在面内部，落在边界线上会返回False。
• intersects包括边界接触，但会导致一个点同时匹配到多个面（因为它正好在共享边上）。

一个点落在两个区的公共边界上，用intersects会得到两条匹配记录，你的统计就重复了；用within会得到零条记录，这个点就丢了。两种都不理想。

我的处理策略是分两步：先用within做主匹配，把所有匹配上的点固定下来；剩下没匹配上的（就是那些落在边界上或者坐标有偏移的），单独用最近邻回填，并且打一个match_type标记。这样既不重复，也不丢数据，而且事后能追溯哪些点是靠回填得来的，数据质量报告里可以单独说明。

6.2 完整链路代码

import geopandas as gpd import pandas as pd # 业务点数据：至少要有 lon, lat biz = pd.read_csv("biz_points.csv") biz_gdf = gpd.GeoDataFrame( biz, geometry=gpd.points_from_xy(biz['lon'], biz['lat']), crs="EPSG:4326" ) # 行政区面数据，只保留需要的字段，减少连接开销 adm = gpd.read_file("admin_county.shp", encoding='utf-8') adm = adm[['code', 'name', 'geometry']].to_crs(biz_gdf.crs) # 第一步：主匹配，严格点内 matched = gpd.sjoin( biz_gdf, adm, how='left', predicate='within' )[['biz_id', 'code', 'name', 'geometry']].copy() matched['match_type'] = matched['code'].notna().map({True: 'inside', False: 'unmatched'}) # 第二步：对未匹配点做最近邻回填 left = matched[matched['match_type'] == 'unmatched'].drop(columns=['code', 'name']).copy() if len(left) > 0: backfill = gpd.sjoin_nearest( left, adm, how='left', max_distance=0.05, # 单位是度，约5公里，按需调整 distance_col='dist_deg' )[['biz_id', 'code', 'name', 'dist_deg']] backfill['match_type'] = 'nearest' matched = matched.merge( backfill[['biz_id', 'code', 'name', 'dist_deg']], on='biz_id', how='left', suffixes=('', '_bf') ) matched['code'] = matched['code'].fillna(matched['code_bf']) matched['name'] = matched['name'].fillna(matched['name_bf']) matched.loc[matched['match_type'] == 'unmatched', 'match_type'] = 'nearest' matched = matched.drop(columns=['code_bf', 'name_bf']) # 覆盖率报告 print(matched['match_type'].value_counts(normalize=True))

max_distance这个参数需要反复试。设太小，边界附近的点回填不上；设太大，明显偏远的点会被硬塞进最近的区，产生错误归属。我一般从目标区域平均区县尺度的十分之一开始试，看未匹配率和回填距离分布，再调整。

6.3 坐标精度和坐标系的隐性错位

还有个坑：业务点的坐标系和面数据的坐标系可能一致，但精度不同。业务坐标如果是经过脱敏处理、偏移了几百米的，那你用严格点内匹配，会有相当比例的点落在错误的区里——尤其是那些靠近边界的点。

判断方法：做一次匹配，然后看match_type为nearest的比例。如果这个比例异常高（比如超过10%），大概率是坐标精度或者坐标系有问题，而不是数据本身写得差。

另一个隐性错位是坐标系同名不同实。同样是标称WGS84的数据，有的是原始采集的，有的是从其他坐标系转换过来的，转换参数不同会带来几十米的偏差。这种偏差在省级尺度看不出问题，在区县边界附近就会导致归属错误。

7. 区划调整与版本管理：数据更新时怎么不被"改名"坑死

7.1 区划调整的四种形态，影响完全不同

行政区划调整不是单一动作，至少有四种形态，每种对数据处理的影响方式不一样：

撤并：两个或多个单元合并成一个。原代码失效，新代码出现。业务数据里如果还存着老代码，直接挂接会全部失配。

拆分：一个单元分成多个。原代码可能保留给其中一个，也可能整体作废。

更名：代码不变，名称变了。这种情况最隐蔽——因为代码能对上，挂接不会报错，但展示出来的名字是旧的，用户会觉得数据没更新。

边界微调：代码和名称都不变，几何变了。这种最难发现，因为代码、名称全部对得上，只有几何比对才能看出来。如果你的业务统计依赖于精确的空间范围，边界微调会导致历史数据和当前数据不可比。

我在一个项目里吃过更名的亏：区划数据更新了一版，代码全部对得上，我就直接替换了。结果页面上展示的还是老名字，用户投诉"信息是错的"。查了半天才发现是名称字段变了但代码没变，而我的关联逻辑是"按代码关联取新名称"——理论上应该取到新名字的，问题出在我缓存了旧版本的数据没有刷新。

7.2 用代码做版本比对的完整方法

两份不同时间的行政区划数据，做一次全字段比对，能清楚看出所有变化：

old = gpd.read_file("admin_2022.shp", encoding='utf-8').set_index('code') new = gpd.read_file("admin_2023.shp", encoding='utf-8').set_index('code') old_codes, new_codes = set(old.index), set(new.index) added = new_codes - old_codes # 新增单元 removed = old_codes - new_codes # 消失单元 common = old_codes & new_codes # 代码未变的单元 # 名称变化 name_changed = [ c for c in common if old.loc[c, 'name'] != new.loc[c, 'name'] ] # 几何变化：用面积差 + 几何哈希双重判断 import hashlib def geo_hash(geom): return hashlib.md5(geom.wkb).hexdigest() geom_changed = [ c for c in common if geo_hash(old.loc[c, 'geometry']) != geo_hash(new.loc[c, 'geometry']) ] print(f"新增 {len(added)} 个，消失 {len(removed)} 个，" f"更名 {len(name_changed)} 个，几何变化 {len(geom_changed)} 个")

用wkb哈希判断几何变化有个前提：两份数据的坐标系必须一致，否则同样的几何会算出不同的哈希。所以比对之前先统一to_crs。

如果只想判断"形状是否变了"而不关心坐标精度带来的微小抖动，那就用面积差做主判断：面积差超过千分之一就算变化，这样能过滤掉重采样带来的噪声。

7.3 业务表该带哪个版本的区划

这是个设计问题，不是技术问题，但影响深远。

业务数据在产生的那一刻，它归属于当时生效的区划。如果你用最新的区划去回溯历史数据，会得到两种失真：一是历史上属于A区的数据被算到B区去了；二是历史上存在的区已经不存在了，数据无处安放。

所以我的建议是：业务表里存区划代码，同时存一个区划版本标识。做当期统计时用当期版本，做历史回溯时用历史版本对应的映射关系。这套机制在建表的时候就要定下来，后期再补成本极高。

如果没有版本标识怎么办？退一步的做法是建一张代码变更映射表，记录"老代码 → 新代码"的对应关系。对于撤并的情况，一个老代码可能对应多个新代码（按面积或人口比例分摊），这个分摊规则要显式写进表里，而不是拍脑袋。

8. 一次区域统计对不上账的完整排查链路

8.1 现象：汇总数字比预期多了百分之十五

说一个真实的排查过程。项目背景是给一批业务数据做区域汇总，口径是县级。跑完之后发现全国汇总的业务量比业务方给的总数多了大约15%。数据量对不上，说明中间有重复计数。

排查顺序是这样的，供参考：

第一步，确认业务源表本身有没有重复行。按主键去重之后重新计数，数字没变，排除数据源问题。

第二步，确认空间连接有没有产生重复。检查sjoin之后的记录数是否大于源表记录数。如果大于，说明一个点匹配到了多个面。我这里用的within，理论上不会重复，但为了排除，还是查了一下——结果确实没有重复。所以问题不在连接环节。

第三步，反向检查行政区数据本身。按代码统计要素数，发现有几个代码出现了两次。再看这几个代码对应的几何，是两个几乎完全相同的面，只是顶点顺序反了。这就是典型的重复要素——同一块区域被存了两遍。

重复要素怎么来的？大概率是数据合并的时候没去重，或者是边界微调之后新旧两个版本都留在了同一个文件里。

第四步，去掉重复后重跑。重复去掉之后，汇总数字和业务方给的总数对上了，误差在百分之一以内，属于正常范围。

整个排查过程花了大半天，但如果一开始就做要素唯一性检查，五分钟就能定位。

8.2 一份可复用的排查清单

从这次经历之后，我把区域统计的排查整理成了一个固定顺序，遇到任何"数字对不上"的情况，按这个顺序走：

这个清单的价值在于顺序不能乱。很多人一上来就怀疑坐标系，花两小时换投影、调参数，最后发现是数据里有重复行。从最简单、最可能的原因开始查，是排查的基本原则。

8.3 几个容易被忽略但很致命的小细节

浮点精度导致的匹配失败。两个面共享一条边，理论上坐标应该完全一致，但一个是116.123456789，一个是116.123456788，浮点比较就不相等。这种问题在大范围数据里极其普遍。解决办法是做精度归一，把所有坐标截断到统一的小数位数（比如6位，约0.1米精度），再比较。

索引重置遗漏导致的错位。用explode()或者sjoin之后如果忘了reset_index，后面做按位置赋值的时候会出现索引对不上的情况，赋值结果全乱但不报错。这个坑我踩过两次，现在养成习惯：凡是做了会改变行数的操作，下一行必写reset_index(drop=True)。

编码在不同工具间不一致。用QGIS导出、用Python读取、再用另一个工具打开，三处的中文可能都不一样。建议在项目里统一约定：存储用UTF-8，读取时显式指定编码，中间不做隐式转换。

多部件面的代表点选错。前面提过centroid的问题，但还有一个变体：用representative_point虽然保证点在面内，但它返回的是"最靠近中心的一个内部点"，对于狭长形状的区域，这个点可能落在边缘附近，拿去做半径测算或者距离计算会有偏差。如果业务上需要"区域的几何中心"，那就要先做explode，取最大部件的centroid，再判断这个点是否落在面内，不落在面内才回退到representative_point。

最后分享一个我在实际项目里总结的习惯：每次处理完行政区划数据，都存一份"处理后快照"和一份处理日志。日志里记清楚用的哪个版本、做了什么清洗、坐标系怎么转的、面积算法是什么。三个月后别人接手，或者你自己回头看，这份日志比代码本身更值钱。数据这东西，处理过程不记录，等于没处理过。