GCJ-02 / BD-09 坐标系偏移诊断与修正
何时用
- 用户说"shp 偏了""对不上影像""坐标系有问题"，且数据来自中国国内（天地图/高德/自然资源/百度/Web 地图）。
- 矢量（标称 EPSG:4326 WGS84）叠加在 S2/Landsat（真 WGS84）上整体偏移约 100~700m。
- 注意：不要先怀疑影像 CRS——S2/Landsat 几乎总是真 WGS84；问题 90% 出在国产矢量被 GCJ-02 加密却标成 WGS84。
诊断步骤（关键：用影像反证，别靠猜）
- 把矢量（红色轮廓）叠加到真 WGS84 影像（如 S2 RGB 合成）上，目视看是否偏离村镇/道路。
- 对同一矢量分别做两种反算，生成三幅叠加对比：
- Original（标称 WGS84 原值）
-gcj02_to_wgs84（GCJ-02→WGS84）
-bd09_to_wgs84（BD-09→WGS84） - 哪幅扣住真实地物，就是哪种加密。通常 GCJ-02 命中（政府/高德/腾讯数据）；BD-09 仅百度源。
坐标变换核心代码（Python，逐顶点）
import math
def _delta(lat, lon):
a=6378245.0; ee=0.00669342162296594323
dlat=_tlat(lon-105.0, lat-35.0); dlng=_tlon(lon-105.0, lat-35.0)
rl=lat/180.0*math.pi; m=math.sin(rl); m=1-ee*m*m; sm=math.sqrt(m)
dlat=(dlat*180.0)/((a*(1-ee))/(m*sm)*math.pi)
dlng=(dlng*180.0)/(a/sm*math.cos(rl)*math.pi)
return dlat, dlng
def gcj02_to_wgs84(lat, lon):
if not (72.004<=lon<=137.8347 and 0.8293<=lat<=55.8271): return lat, lon
dlat,dlng=_delta(lat,lon); wlat,wlon=lat-dlat,lon-dlng
for _ in range(5):
dlat,dlng=_delta(wlat,wlon); wlat,wlon=lat-dlat,lon-dlng
return wlat,wlon
# _tlat/_tlon 为 GCJ-02 公开多项式（见下方附注），对每个多边形顶点 (x=lon,y=lat) 调用 gcj02_to_wgs84(y,x)
对 shapely Polygon/MultiPolygon，遍历 exterior.coords 与每个 interior，逐点反算后重建几何。CGCS2000 与 WGS84 差异 <2m，可忽略，无需处理。
> _tlat(lng,lat)、_tlon(lng,lat) 为 GCJ-02 标准多项式（网上公开实现一致），含 sin/cos 多项与 sqrt(abs)，直接套用即可。
修正流程
- 读原始 shp（标称 4326）→ 逐顶点
gcj02_to_wgs84→ 写新 shp（仍是set_crs('EPSG:4326')，但现在是真 WGS84）。保留原文件不覆盖。 - 校验：全部
geometry.is_valid（通常仍有效，变换平滑）；可算平均 centroid 偏移米数（山东实测 ~500m）作为修正量证据。 - 下游（掩膜/分类/面积统计）一律用修正后 shp。
★★★ 反向陷阱：_wgs84 文件不是权威，必须先判定再使用
> 实战事故（2026-09-17，某项目区）：项目里同时存在 村界.shp（原始）和
> 村界_wgs84.shp（由 fix_village_crs_liyang.py 做 GCJ-02→WGS84"纠偏"产出）。
> 纠偏脚本本身写错了：它把已是真 WGS84 的原始村界又当 GCJ-02 转了一次，
> 于是 _wgs84 反而偏了 495 m。
> 更严重的是：上一版已交付的报告全流程用了这个错文件（build_report_stats_fix9.py
> / make_village_stats_fix9.py 都 read_file('村界_wgs84.shp')），
> 导致每个村的数字都"搬"到隔壁村。直到一个月后两版报告对不上才暴露。
教训：文件名里的 _wgs84 / _fix / _corrected 不构成任何证据。见到这类文件，
先做第 4 节的三步判定，再决定用哪套。
4 两套矢量并存时，如何判定哪套正确（4 组独立判据）
前提：两套键（镇+村）完全一致、只是几何整体错开。不要用"看起来对"定案，
下面 4 条都不依赖任何一套矢量自身：
判据 ① 贴河率（中国村界最强判据）
绝大多数村界以河为界。用 S2 L2A 的 SCL 层 class=6（水体） 或 JRC 水体，
统计界线附近的水体像元占比。正确的那套会高出一个量级。
实测：距水 ≤30m 的边界像元 正确版 9.6% vs 错误版 2.4% （4 倍差）
判据 ② 界线切穿建成区（越少越好）
用 ESA WorldCover 2021 class=50（built-up），统计村界线上（含 ±10m/±20m 膨胀）的建成区像元。
实测：距界 ≤0m 正确 219.0 亩 vs 错误 326.4 亩
距界 ≤20m 正确 842.5 亩 vs 错误 1290.5 亩
界线切穿村庄居民点说明它错位了。
判据 ③ 已交付成果是被哪套裁的（追溯传播面）
把既有交付层（作物/地块 shp）分别与两套的并集求 clip，看落在谁里面比例更高。
实测：交付层 430808.9 亩 → 错误版并集内 99.92% ，正确版并集内 96.76%
→ 交付层是按错误版裁的 ⇒ 交付层也需要重裁
这条判据同时暴露下游污染范围，是排查必备。
判据 ④ 高倍目视（1200×1200 栅格窗口）
红/青双色叠图，放大到河段级别。正确的那套严丝合缝沿河道走，错误的那套斜切田地。
补充：排除"栅格错位"的干扰
若两套成果各有栅格，先做 IoU + 平移自检（np.roll dy/dx = −2…+2）：
零平移时 IoU 最优（实测 0.855，任一方向平移都下降）⇒ 无重投影错位，
差异只来自矢量裁剪，可以放心用矢量判据定案。
★ 判据 ⑤ 官方承包地块质心命中率（有地块 shp 时最快的判据，秒级定案）
农险项目里常同时拿到官方承包地块 shp（如 Parcels_*.shp，EPSG:4490/4553）。
它是逐块实测矢量，与影像无关 → 可当独立参照用，不需要打开影像目视。
vil = vil.to_crs(32651); par = par.to_crs(32651)
u = unary_union(vil[vil['zhen']==镇].geometry.values)
cen = par.geometry.representative_point() # 用质心，不用交面积（抗边界抖动）
print('质心命中率 %.2f%%' % (100*cen.within(u).mean()))
print('面积占比 %.2f%%' % (100*par.geometry.intersection(u).area.sum()/par.geometry.area.sum()))
判读：命中率与面积占比都 >95% ⇒ 村界与地块同一坐标系（真 WGS84），可用；
若出现 整体错开数百米 / 命中率骤降到 <50% ⇒ 村界是 GCJ-02，走上面 4 组判据定案。
为什么用质心而不是交面积：地块沿村界密集排列，边界抖动会同时吃进/吐出面积、
把真信号淹没；质心是离散的是/否判定，跨界地块只占极少数，信噪比高得多。
实测（某项目区，2026-09-20）：某镇 11,992 块官方地块（EPSG:4490）
对 <村界目录>\村界.shp 的某镇村界 → 质心命中 99.46% / 面积占比 99.60% ⇒ 一次判定通过，
不必再开影像目视。这一条应在拿到新村界后第一时间跑，成本 3 秒。
反向验证（拿到"谁对"的结论后必做）
用判为正确的那套重算任一份既有报告的县级汇总，与报告原文逐档对比。
若用错误的那套能复现出 差 0.0，则"报告是用错误的那套算的"这个结论闭合。
5 纠正后的传播面排查（必做）
grep -rl "村界_wgs84\|_wgs84\.shp" --include=*.py --include=*.json .
- 列出所有引用错文件的下游脚本 / 报告 / 图件（含
.mxd、出图任务表）→ 全部需要重跑 - 逐项目检查其他地区是否也有一份错文件（本例：只有某市受影响，其余各市只有一套村界，未受影响）
- 报告里的图件与数字是同源的：换村界必须同批重出，只订正数字会图文不一致
易错点
- 样本点若是在 ArcGIS 里按真 WGS84 影像打点的，坐标正确，无需修正；只有标称 WGS84 的国产矢量才需要修。
- 别对已经纠偏过的数据再纠偏一次——迭代次数写成 3~5 次是给单次反算收敛用的，
不是让你对同一份数据反复调用。动手前先判定原始文件是 GCJ-02 还是真 WGS84。 - 别用
os.remove删原文件；沙箱会拦。直接写_wgs84新文件。 - 一个 tile 内偏移量随经纬度略变（非严格平移），但量级一致，整体视觉修正即够用。
- 逐村质心偏移量应作为诊断输出打印出来（中位/极值/方位角）。若中位 400~550m 且
71/71 村全部错开，就是典型的 GCJ-02 量级；若只有个别村偏，那是数据质量问题，不是坐标系问题。
产出清单示例（山东项目）
- 修正后村界
merged_villages_wgs84.shp - 选景映射 JSON、逐村面积 CSV
- 诊断图
alignment_compare.png（original/GCJ-02/BD-09 三色叠加）
