行政区划shp数据这几个字看着平平无奇但凡做过区域类项目的人都清楚它几乎是最容易拿到、也最容易用歪的一份基础数据。下载下来解压一堆同名不同后缀的文件摊在文件夹里拖进桌面软件能正常显示颜色一配、图例一挂看起来就能用了。可一旦你要拿它去挂业务数据、做区域统计、做看板下钻、做面积测算前面偷的每一步懒都会在后面变成对不上账的麻烦。我这些年在区域分析、门店选址、物流分单、看板地图这几类活儿里反复跟这份数据打交道踩过的坑从编码乱码到坐标系误判从共享边不重合到区划调整导致的代码失效几乎把能踩的都踩了一遍。这篇就把这套经验摊开讲从拿到文件那一刻该做什么体检到层级代码怎么设计主键到投影怎么选、面积怎么算再到空间连接和版本回溯的完整链路。不管你是刚接触地理数据的开发还是被区域统计折磨过一轮的分析师下面这些内容应该都能直接拿去用。1. 拿到行政区划shp先别急着加载先想清楚你要用它回答什么问题1.1 三种用法对应三套完全不同的处理路径很多人把行政区划shp数据当成一张底图这个定位本身就限制了它的价值。它其实是一套带层级关系的空间主数据用法的差异会直接决定后面所有技术选择。制图展示型只关心形状好看不好看。这种场景下坐标系随便Web墨卡托就行拓扑有点小缝隙也看不出来字段乱码更是无所谓反正只画不标。这是成本最低的用法也是唯一的拖进去就能用的用法。统计挂接型才是主流需求。比如把订单表按区域汇总、把门店地址归到某个区县、把用户数按行政区聚合。这种场景的核心矛盾是主键唯一性和名称标准化——区划代码必须干净、不能重复、不能有前后空格名称必须能对上业务表里的写法。这时候坐标系反而次要因为你多半是在做属性关联而不是几何运算。分析计算型最苛刻。面积、人口密度、辐射范围、缓冲区覆盖、区域间距离这些都要算数值。一旦涉及数值坐标系选择就成了生死线投影错了算出来的面积能差出一个量级拓扑不干净叠加分析会产生一堆碎渣多边形。我的习惯是在打开文件之前先用一句话写清楚这个项目最终要输出什么。是一张能看的地图还是一张区域汇总表还是一组空间测算指标。写不出来就先别动数据因为你根本不知道自己该清洗到什么程度。1.2 把所有层级合并成一个图层是最常见的第一个错误新手最容易干的一件事是把省级、市级、县级三个shp合并成一个大图层然后一次性渲染觉得这样图上有省界也有县界很完整。问题在于行政区划是嵌套关系。一个省的面几何上就包含了它下属所有市的面。你把省和市叠在同一个图层里同一个空间位置会被重复绘制两到三次。视觉效果上你可能觉得边界线变粗了但实际后果是任何基于这个图层的统计都会重复计数。举个具体场景。你想统计全国有多少个区域单元合并后的图层会给你省数市数县数这个数字没有任何业务含义。你想算总面积三个层级的面叠在一起总面积会是实际面积的两到三倍。正确的做法是分层存放、按需调用。省级一份、市级一份、县级一份各自独立成层需要哪个层级就加载哪个。确实需要在一个视图里同时展示两级边界时用样式区分省界粗线、县界细线而不是把数据物理合并。1.3 先确认这份数据的时间截面行政区划不是静态的。撤县设区、合并、新设、更名每年都在发生。一份数据里如果混了不同年份的调整结果会出现非常诡异的状况某个代码在省级图层里存在在县级图层里却找不到对应的父级。我在一个项目里遇到过市级图层里有一个单元代码前缀在省级图层里查不到。排查了半小时才反应过来这是两个不同批次的数据拼在一起的中间正好赶上一轮调整。所以拿到数据第一件事翻一下有没有元数据、有没有说明文档、有没有标注数据截止时间。没有的话就自己记一个入库时间后面所有对不上的问题都从这里开始回溯。2. Shapefile四件套与DBF字段限制决定了你逃不掉的一轮体检2.1 .shp、.shx、.dbf、.prj 各管什么缺一个会怎样Shapefile从来不是一个文件它是至少四个文件组成的集合。这个设计来自很早的年代一直沿用至今理解它各自的职责能省掉很多莫名其妙的故障。文件后缀作用缺失后的典型症状.shp存储几何坐标本身文件根本打不开.shx几何索引记录每个要素在shp中的偏移能打开但读取极慢或直接报索引错误.dbf属性表存代码、名称等字段只剩图形没有任何属性无法关联.prj坐标系定义WKT文本软件按猜的方式套坐标系位置可能整体偏移.prj缺失是最阴险的一种。因为软件不会报错它会默默按某种默认坐标系渲染图上看起来形状是对的你根本发现不了。等到你把这份数据和另一份带正确坐标系的数据叠在一起才发现整体错位了几百公里。我见过的补救办法是从同一批数据的其他图层里复制一份.prj过来前提是你确认它们坐标系一致。如果连参照都没有只能靠已知地标的实际经纬度去反推成本很高。还有个小细节文件名千万别带中文和空格。某些老版本驱动在路径解析上会出问题尤其是路径里带空格的时候读出来是空表但不是报错排查起来很费时间。养成习惯解压后先重命名成英文短名。2.2 DBF的字段名10字符上限与中文编码DBF这个格式的字段名上限是10个字符。这意味着district_name会被截断成district_nadministrative_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, encodingutf-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, encodingencoding) 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, encodingutf-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] / 1e64.2 等积、等距、等角选错了就是白算投影不是随便挑一个能用的不同投影保留的特性完全不同。做区域统计你需要的是等积投影做距离测算你需要的是等距投影做导航和切片底图等角投影才合适。投影类型保留特性适合场景不适合场景等积圆锥投影面积比例准确全国/大区域面积统计、人口密度形状会变形不适合精细制图等距投影特定方向距离准确辐射范围、服务半径测算面积会有偏差等角投影如墨卡托族角度和形状不变瓦片底图、导航、形状展示高纬度面积严重放大绝对不能算面积通用横轴分带投影局部精度高单个中小区域内的精细测量跨带拼接麻烦大范围不适用Web墨卡托是面积计算的重灾区。它为了保持形状和角度把高纬度区域在纵向拉伸得越来越厉害。在它的坐标系下算面积纬度越高误差越离谱越靠近极地偏差越夸张。而偏偏很多在线地图服务和前端组件默认就用这个投影好多人直接在上面算面积结果当然是错的。跨越大范围区域做等积投影时参数怎么定常见做法是标准纬线取目标区域南北跨度的约六分之一和六分之五位置中央经线取区域中轴经度。这样能让整个区域内的面积变形控制在很小的范围内。4.3 面积计算的完整写法与误差校验import geopandas as gpd from pyproj import CRS # 假设目标区域跨纬度 20~50 度中轴经度 105 度 albers CRS.from_proj4( projaea lat_125 lat_245 lat_00 lon_0105 x_00 y_00 datumWGS84 unitsm 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_partsFalse).reset_index(dropTrue) 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, geometrygpd.points_from_xy(biz[lon], biz[lat]), crsEPSG:4326 ) # 行政区面数据只保留需要的字段减少连接开销 adm gpd.read_file(admin_county.shp, encodingutf-8) adm adm[[code, name, geometry]].to_crs(biz_gdf.crs) # 第一步主匹配严格点内 matched gpd.sjoin( biz_gdf, adm, howleft, predicatewithin )[[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, howleft, max_distance0.05, # 单位是度约5公里按需调整 distance_coldist_deg )[[biz_id, code, name, dist_deg]] backfill[match_type] nearest matched matched.merge( backfill[[biz_id, code, name, dist_deg]], onbiz_id, howleft, 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(normalizeTrue))max_distance这个参数需要反复试。设太小边界附近的点回填不上设太大明显偏远的点会被硬塞进最近的区产生错误归属。我一般从目标区域平均区县尺度的十分之一开始试看未匹配率和回填距离分布再调整。6.3 坐标精度和坐标系的隐性错位还有个坑业务点的坐标系和面数据的坐标系可能一致但精度不同。业务坐标如果是经过脱敏处理、偏移了几百米的那你用严格点内匹配会有相当比例的点落在错误的区里——尤其是那些靠近边界的点。判断方法做一次匹配然后看match_type为nearest的比例。如果这个比例异常高比如超过10%大概率是坐标精度或者坐标系有问题而不是数据本身写得差。另一个隐性错位是坐标系同名不同实。同样是标称WGS84的数据有的是原始采集的有的是从其他坐标系转换过来的转换参数不同会带来几十米的偏差。这种偏差在省级尺度看不出问题在区县边界附近就会导致归属错误。7. 区划调整与版本管理数据更新时怎么不被改名坑死7.1 区划调整的四种形态影响完全不同行政区划调整不是单一动作至少有四种形态每种对数据处理的影响方式不一样撤并两个或多个单元合并成一个。原代码失效新代码出现。业务数据里如果还存着老代码直接挂接会全部失配。拆分一个单元分成多个。原代码可能保留给其中一个也可能整体作废。更名代码不变名称变了。这种情况最隐蔽——因为代码能对上挂接不会报错但展示出来的名字是旧的用户会觉得数据没更新。边界微调代码和名称都不变几何变了。这种最难发现因为代码、名称全部对得上只有几何比对才能看出来。如果你的业务统计依赖于精确的空间范围边界微调会导致历史数据和当前数据不可比。我在一个项目里吃过更名的亏区划数据更新了一版代码全部对得上我就直接替换了。结果页面上展示的还是老名字用户投诉信息是错的。查了半天才发现是名称字段变了但代码没变而我的关联逻辑是按代码关联取新名称——理论上应该取到新名字的问题出在我缓存了旧版本的数据没有刷新。7.2 用代码做版本比对的完整方法两份不同时间的行政区划数据做一次全字段比对能清楚看出所有变化old gpd.read_file(admin_2022.shp, encodingutf-8).set_index(code) new gpd.read_file(admin_2023.shp, encodingutf-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 一份可复用的排查清单从这次经历之后我把区域统计的排查整理成了一个固定顺序遇到任何数字对不上的情况按这个顺序走顺序检查项判断方法常见原因1源数据是否重复按业务主键 count 对比去重后 count数据导出重复、多次追加2空间连接是否重复sjoin 后行数 源表行数用了 intersects、点在边界上3行政区要素是否重复按代码 groupby count数据合并未去重、多版本混存4层级是否混用检查是否同时存在省市县三种层级合并图层后直接统计5是否有孤儿记录连接后 code 为空的数量名称不一致、代码格式错6面积/总量是否合理汇总值与已知总量对比坐标系错、投影选错7几何是否重叠union 前后面积对比拓扑没处理这个清单的价值在于顺序不能乱。很多人一上来就怀疑坐标系花两小时换投影、调参数最后发现是数据里有重复行。从最简单、最可能的原因开始查是排查的基本原则。8.3 几个容易被忽略但很致命的小细节浮点精度导致的匹配失败。两个面共享一条边理论上坐标应该完全一致但一个是116.123456789一个是116.123456788浮点比较就不相等。这种问题在大范围数据里极其普遍。解决办法是做精度归一把所有坐标截断到统一的小数位数比如6位约0.1米精度再比较。索引重置遗漏导致的错位。用explode()或者sjoin之后如果忘了reset_index后面做按位置赋值的时候会出现索引对不上的情况赋值结果全乱但不报错。这个坑我踩过两次现在养成习惯凡是做了会改变行数的操作下一行必写reset_index(dropTrue)。编码在不同工具间不一致。用QGIS导出、用Python读取、再用另一个工具打开三处的中文可能都不一样。建议在项目里统一约定存储用UTF-8读取时显式指定编码中间不做隐式转换。多部件面的代表点选错。前面提过centroid的问题但还有一个变体用representative_point虽然保证点在面内但它返回的是最靠近中心的一个内部点对于狭长形状的区域这个点可能落在边缘附近拿去做半径测算或者距离计算会有偏差。如果业务上需要区域的几何中心那就要先做explode取最大部件的centroid再判断这个点是否落在面内不落在面内才回退到representative_point。最后分享一个我在实际项目里总结的习惯每次处理完行政区划数据都存一份处理后快照和一份处理日志。日志里记清楚用的哪个版本、做了什么清洗、坐标系怎么转的、面积算法是什么。三个月后别人接手或者你自己回头看这份日志比代码本身更值钱。数据这东西处理过程不记录等于没处理过。
行政区划shp数据这几个字看着平平无奇但凡做过区域类项目的人都清楚它几乎是最容易拿到、也最容易用歪的一份基础数据。下载下来解压一堆同名不同后缀的文件摊在文件夹里拖进桌面软件能正常显示颜色一配、图例一挂看起来就能用了。可一旦你要拿它去挂业务数据、做区域统计、做看板下钻、做面积测算前面偷的每一步懒都会在后面变成对不上账的麻烦。我这些年在区域分析、门店选址、物流分单、看板地图这几类活儿里反复跟这份数据打交道踩过的坑从编码乱码到坐标系误判从共享边不重合到区划调整导致的代码失效几乎把能踩的都踩了一遍。这篇就把这套经验摊开讲从拿到文件那一刻该做什么体检到层级代码怎么设计主键到投影怎么选、面积怎么算再到空间连接和版本回溯的完整链路。不管你是刚接触地理数据的开发还是被区域统计折磨过一轮的分析师下面这些内容应该都能直接拿去用。1. 拿到行政区划shp先别急着加载先想清楚你要用它回答什么问题1.1 三种用法对应三套完全不同的处理路径很多人把行政区划shp数据当成一张底图这个定位本身就限制了它的价值。它其实是一套带层级关系的空间主数据用法的差异会直接决定后面所有技术选择。制图展示型只关心形状好看不好看。这种场景下坐标系随便Web墨卡托就行拓扑有点小缝隙也看不出来字段乱码更是无所谓反正只画不标。这是成本最低的用法也是唯一的拖进去就能用的用法。统计挂接型才是主流需求。比如把订单表按区域汇总、把门店地址归到某个区县、把用户数按行政区聚合。这种场景的核心矛盾是主键唯一性和名称标准化——区划代码必须干净、不能重复、不能有前后空格名称必须能对上业务表里的写法。这时候坐标系反而次要因为你多半是在做属性关联而不是几何运算。分析计算型最苛刻。面积、人口密度、辐射范围、缓冲区覆盖、区域间距离这些都要算数值。一旦涉及数值坐标系选择就成了生死线投影错了算出来的面积能差出一个量级拓扑不干净叠加分析会产生一堆碎渣多边形。我的习惯是在打开文件之前先用一句话写清楚这个项目最终要输出什么。是一张能看的地图还是一张区域汇总表还是一组空间测算指标。写不出来就先别动数据因为你根本不知道自己该清洗到什么程度。1.2 把所有层级合并成一个图层是最常见的第一个错误新手最容易干的一件事是把省级、市级、县级三个shp合并成一个大图层然后一次性渲染觉得这样图上有省界也有县界很完整。问题在于行政区划是嵌套关系。一个省的面几何上就包含了它下属所有市的面。你把省和市叠在同一个图层里同一个空间位置会被重复绘制两到三次。视觉效果上你可能觉得边界线变粗了但实际后果是任何基于这个图层的统计都会重复计数。举个具体场景。你想统计全国有多少个区域单元合并后的图层会给你省数市数县数这个数字没有任何业务含义。你想算总面积三个层级的面叠在一起总面积会是实际面积的两到三倍。正确的做法是分层存放、按需调用。省级一份、市级一份、县级一份各自独立成层需要哪个层级就加载哪个。确实需要在一个视图里同时展示两级边界时用样式区分省界粗线、县界细线而不是把数据物理合并。1.3 先确认这份数据的时间截面行政区划不是静态的。撤县设区、合并、新设、更名每年都在发生。一份数据里如果混了不同年份的调整结果会出现非常诡异的状况某个代码在省级图层里存在在县级图层里却找不到对应的父级。我在一个项目里遇到过市级图层里有一个单元代码前缀在省级图层里查不到。排查了半小时才反应过来这是两个不同批次的数据拼在一起的中间正好赶上一轮调整。所以拿到数据第一件事翻一下有没有元数据、有没有说明文档、有没有标注数据截止时间。没有的话就自己记一个入库时间后面所有对不上的问题都从这里开始回溯。2. Shapefile四件套与DBF字段限制决定了你逃不掉的一轮体检2.1 .shp、.shx、.dbf、.prj 各管什么缺一个会怎样Shapefile从来不是一个文件它是至少四个文件组成的集合。这个设计来自很早的年代一直沿用至今理解它各自的职责能省掉很多莫名其妙的故障。文件后缀作用缺失后的典型症状.shp存储几何坐标本身文件根本打不开.shx几何索引记录每个要素在shp中的偏移能打开但读取极慢或直接报索引错误.dbf属性表存代码、名称等字段只剩图形没有任何属性无法关联.prj坐标系定义WKT文本软件按猜的方式套坐标系位置可能整体偏移.prj缺失是最阴险的一种。因为软件不会报错它会默默按某种默认坐标系渲染图上看起来形状是对的你根本发现不了。等到你把这份数据和另一份带正确坐标系的数据叠在一起才发现整体错位了几百公里。我见过的补救办法是从同一批数据的其他图层里复制一份.prj过来前提是你确认它们坐标系一致。如果连参照都没有只能靠已知地标的实际经纬度去反推成本很高。还有个小细节文件名千万别带中文和空格。某些老版本驱动在路径解析上会出问题尤其是路径里带空格的时候读出来是空表但不是报错排查起来很费时间。养成习惯解压后先重命名成英文短名。2.2 DBF的字段名10字符上限与中文编码DBF这个格式的字段名上限是10个字符。这意味着district_name会被截断成district_nadministrative_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, encodingutf-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, encodingencoding) 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, encodingutf-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] / 1e64.2 等积、等距、等角选错了就是白算投影不是随便挑一个能用的不同投影保留的特性完全不同。做区域统计你需要的是等积投影做距离测算你需要的是等距投影做导航和切片底图等角投影才合适。投影类型保留特性适合场景不适合场景等积圆锥投影面积比例准确全国/大区域面积统计、人口密度形状会变形不适合精细制图等距投影特定方向距离准确辐射范围、服务半径测算面积会有偏差等角投影如墨卡托族角度和形状不变瓦片底图、导航、形状展示高纬度面积严重放大绝对不能算面积通用横轴分带投影局部精度高单个中小区域内的精细测量跨带拼接麻烦大范围不适用Web墨卡托是面积计算的重灾区。它为了保持形状和角度把高纬度区域在纵向拉伸得越来越厉害。在它的坐标系下算面积纬度越高误差越离谱越靠近极地偏差越夸张。而偏偏很多在线地图服务和前端组件默认就用这个投影好多人直接在上面算面积结果当然是错的。跨越大范围区域做等积投影时参数怎么定常见做法是标准纬线取目标区域南北跨度的约六分之一和六分之五位置中央经线取区域中轴经度。这样能让整个区域内的面积变形控制在很小的范围内。4.3 面积计算的完整写法与误差校验import geopandas as gpd from pyproj import CRS # 假设目标区域跨纬度 20~50 度中轴经度 105 度 albers CRS.from_proj4( projaea lat_125 lat_245 lat_00 lon_0105 x_00 y_00 datumWGS84 unitsm 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_partsFalse).reset_index(dropTrue) 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, geometrygpd.points_from_xy(biz[lon], biz[lat]), crsEPSG:4326 ) # 行政区面数据只保留需要的字段减少连接开销 adm gpd.read_file(admin_county.shp, encodingutf-8) adm adm[[code, name, geometry]].to_crs(biz_gdf.crs) # 第一步主匹配严格点内 matched gpd.sjoin( biz_gdf, adm, howleft, predicatewithin )[[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, howleft, max_distance0.05, # 单位是度约5公里按需调整 distance_coldist_deg )[[biz_id, code, name, dist_deg]] backfill[match_type] nearest matched matched.merge( backfill[[biz_id, code, name, dist_deg]], onbiz_id, howleft, 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(normalizeTrue))max_distance这个参数需要反复试。设太小边界附近的点回填不上设太大明显偏远的点会被硬塞进最近的区产生错误归属。我一般从目标区域平均区县尺度的十分之一开始试看未匹配率和回填距离分布再调整。6.3 坐标精度和坐标系的隐性错位还有个坑业务点的坐标系和面数据的坐标系可能一致但精度不同。业务坐标如果是经过脱敏处理、偏移了几百米的那你用严格点内匹配会有相当比例的点落在错误的区里——尤其是那些靠近边界的点。判断方法做一次匹配然后看match_type为nearest的比例。如果这个比例异常高比如超过10%大概率是坐标精度或者坐标系有问题而不是数据本身写得差。另一个隐性错位是坐标系同名不同实。同样是标称WGS84的数据有的是原始采集的有的是从其他坐标系转换过来的转换参数不同会带来几十米的偏差。这种偏差在省级尺度看不出问题在区县边界附近就会导致归属错误。7. 区划调整与版本管理数据更新时怎么不被改名坑死7.1 区划调整的四种形态影响完全不同行政区划调整不是单一动作至少有四种形态每种对数据处理的影响方式不一样撤并两个或多个单元合并成一个。原代码失效新代码出现。业务数据里如果还存着老代码直接挂接会全部失配。拆分一个单元分成多个。原代码可能保留给其中一个也可能整体作废。更名代码不变名称变了。这种情况最隐蔽——因为代码能对上挂接不会报错但展示出来的名字是旧的用户会觉得数据没更新。边界微调代码和名称都不变几何变了。这种最难发现因为代码、名称全部对得上只有几何比对才能看出来。如果你的业务统计依赖于精确的空间范围边界微调会导致历史数据和当前数据不可比。我在一个项目里吃过更名的亏区划数据更新了一版代码全部对得上我就直接替换了。结果页面上展示的还是老名字用户投诉信息是错的。查了半天才发现是名称字段变了但代码没变而我的关联逻辑是按代码关联取新名称——理论上应该取到新名字的问题出在我缓存了旧版本的数据没有刷新。7.2 用代码做版本比对的完整方法两份不同时间的行政区划数据做一次全字段比对能清楚看出所有变化old gpd.read_file(admin_2022.shp, encodingutf-8).set_index(code) new gpd.read_file(admin_2023.shp, encodingutf-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 一份可复用的排查清单从这次经历之后我把区域统计的排查整理成了一个固定顺序遇到任何数字对不上的情况按这个顺序走顺序检查项判断方法常见原因1源数据是否重复按业务主键 count 对比去重后 count数据导出重复、多次追加2空间连接是否重复sjoin 后行数 源表行数用了 intersects、点在边界上3行政区要素是否重复按代码 groupby count数据合并未去重、多版本混存4层级是否混用检查是否同时存在省市县三种层级合并图层后直接统计5是否有孤儿记录连接后 code 为空的数量名称不一致、代码格式错6面积/总量是否合理汇总值与已知总量对比坐标系错、投影选错7几何是否重叠union 前后面积对比拓扑没处理这个清单的价值在于顺序不能乱。很多人一上来就怀疑坐标系花两小时换投影、调参数最后发现是数据里有重复行。从最简单、最可能的原因开始查是排查的基本原则。8.3 几个容易被忽略但很致命的小细节浮点精度导致的匹配失败。两个面共享一条边理论上坐标应该完全一致但一个是116.123456789一个是116.123456788浮点比较就不相等。这种问题在大范围数据里极其普遍。解决办法是做精度归一把所有坐标截断到统一的小数位数比如6位约0.1米精度再比较。索引重置遗漏导致的错位。用explode()或者sjoin之后如果忘了reset_index后面做按位置赋值的时候会出现索引对不上的情况赋值结果全乱但不报错。这个坑我踩过两次现在养成习惯凡是做了会改变行数的操作下一行必写reset_index(dropTrue)。编码在不同工具间不一致。用QGIS导出、用Python读取、再用另一个工具打开三处的中文可能都不一样。建议在项目里统一约定存储用UTF-8读取时显式指定编码中间不做隐式转换。多部件面的代表点选错。前面提过centroid的问题但还有一个变体用representative_point虽然保证点在面内但它返回的是最靠近中心的一个内部点对于狭长形状的区域这个点可能落在边缘附近拿去做半径测算或者距离计算会有偏差。如果业务上需要区域的几何中心那就要先做explode取最大部件的centroid再判断这个点是否落在面内不落在面内才回退到representative_point。最后分享一个我在实际项目里总结的习惯每次处理完行政区划数据都存一份处理后快照和一份处理日志。日志里记清楚用的哪个版本、做了什么清洗、坐标系怎么转的、面积算法是什么。三个月后别人接手或者你自己回头看这份日志比代码本身更值钱。数据这东西处理过程不记录等于没处理过。
