Files
CM95/cm93_mysql_normalization_canonical_v10.py
2026-09-02 09:50:12 +08:00

1634 lines
47 KiB
Python
Raw Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
CM93 -> MySQL OGC Normalization 可重复性 / Canonical WKB 验证
Version: v10.0
背景
====
V9 已确认:
- 55/55 invalid Area/Coverage 使用 odd-even face selection 可生成 valid Polygon/MultiPolygon;
- 55/55 与 make_valid polygonal 部分一致;
- 55/55 不引入原始 linework 之外的新边界;
- 55/55 被丢弃的 linework 与已确认 retraced connector 完全一致;
- 55/55 normalization 幂等;
- 82181 条原本有效的 Area/Coverage 全部 PASS-THROUGH,WKB hash 0 变化;
- 27 条 candidate bbox 缩小,其中 27/27 都由被去掉的零面积 retraced connector 解释。
V10 目的
=======
在真正设计 importer/schema 修改前,再验证一件实现层关键要求:
“相同的 CM93 面语义,在输入 ring 起点、ring 方向、MultiPolygon component 顺序、
polygonize face 顺序不同的情况下,最终是否能得到稳定的 canonical geometry / WKB hash。”
原因:
如果 normalized geometry 语义完全相同,但 WKB 因 ring 起点或 component 顺序变化而不同,
以后重新导入、增量同步、缓存或 hash 对账可能出现伪变化。
V10 只读、只在内存运行。
验证内容
========
对 55 个 invalid Area/Coverage:
1. 基线:
- original invalid geometry
- odd-even normalization candidate
- GEOS/Shapely normalize() canonicalization
- canonical WKB SHA-256
2. 输入等价变体:
A. 原始输入
B. 每个 ring 全部反向
C. 每个 ring 改变起始 vertex
D. ring 反向 + 改变起点
E. MultiPolygon component 反序(如适用)
以上变体不改变原始线网/odd-even 面语义。
3. 对每个变体重新执行 odd-even normalization;
然后 canonicalize;
必须满足:
- candidate valid
- candidate polygonal-only
- 与基线 candidate topology equals
- canonical WKB SHA-256 与基线一致
4. polygonize selected-face 顺序:
- 原顺序 union
- reverse 顺序 union
- 按面积排序 union
- 按 centroid 排序 union
canonical hash 必须一致。
5. 重复执行:
同一输入连续 normalization 3 次,canonical hash 必须一致。
6. Canonicalization 自身不得改变 topology:
- symdiff area = 0
- Hausdorff distance = 0(数值容差内)
7. 输出 55 条 golden vectors:
- kind
- feature_id / coverage_id
- object_class
- cell
- source_record_index
- normalized type
- normalized part count
- canonical WKB SHA-256
- candidate bbox
- candidate area
- discarded connector length
- feature bbox 是否未来需要更新
这些 golden vectors 可供后续 importer 实现回归测试。
8. 有效 geometry policy:
仍然只做 PASS-THROUGH,不 canonicalize 后写回;
统计 valid records,确认 policy 不会改它们。
重要:
- CM93JAPAN 是日本沿海局部数据,不是完整 CM93;
- 跨未导入海域边界不连续不作为错误;
- 本脚本不修改数据库;
- 不开始 PBF。
依赖
====
pip install pymysql shapely
当前验证环境已知 Shapely 2.0.7。
执行
====
python cm93_mysql_normalization_canonical_v10.py ^
--host 127.0.0.1 ^
--port 3306 ^
--user root ^
--database cm93_chart ^
--dataset-id 1
输出
====
- normalization_canonical_summary_v10.txt
- normalization_canonical_report_v10.json
- normalization_golden_vectors_v10.csv
- normalization_variant_results.csv
- normalization_face_order_results.csv
- normalization_canonical_failures.csv
- normalization_bbox_metadata_impact.csv
- normalization_canonical_candidates.geojson
- normalization_golden_vectors_v10.json
并自动生成 ZIP。
禁止
====
- 不 UPDATE / INSERT / DELETE
- 不改 importer
- 不写 normalized geometry
- 不修改 bbox
- 不 PBF
"""
from __future__ import annotations
import argparse
import csv
import getpass
import hashlib
import json
import math
import os
import re
import sys
import traceback
import zipfile
from collections import Counter
from datetime import datetime
from pathlib import Path
from typing import Any, Dict, Iterable, List, Optional, Sequence, Tuple
try:
import pymysql
from pymysql.cursors import DictCursor
except ImportError:
print("缺少 pymysql:pip install pymysql shapely")
sys.exit(2)
try:
import shapely
try:
from shapely import from_wkb as _from_wkb
def load_wkb(data: bytes):
return _from_wkb(data)
except ImportError:
from shapely import wkb as _wkb
def load_wkb(data: bytes):
return _wkb.loads(data)
try:
from shapely import to_wkb as _to_wkb
def dump_wkb(g) -> bytes:
return bytes(_to_wkb(g, hex=False, byte_order=1))
except ImportError:
from shapely import wkb as _wkb2
def dump_wkb(g) -> bytes:
return _wkb2.dumps(g, hex=False, big_endian=False)
try:
from shapely import normalize as _normalize
def canonicalize(g):
return _normalize(g)
except ImportError:
raise ImportError("V10 需要 Shapely 2.x 的 shapely.normalize()")
try:
from shapely import make_valid as _make_valid
def make_valid_geom(g):
return _make_valid(g)
except ImportError:
from shapely.validation import make_valid as _make_valid
def make_valid_geom(g):
return _make_valid(g)
from shapely.geometry import (
Point, Polygon, MultiPolygon, mapping
)
from shapely.ops import unary_union, polygonize_full
from shapely.validation import explain_validity
except ImportError as e:
print(f"Shapely 初始化失败:{e}")
print("请执行:pip install -U shapely")
sys.exit(2)
INVALID_POINT_RE = re.compile(
r"\[\s*([-+]?\d+(?:\.\d+)?(?:[eE][-+]?\d+)?)\s+"
r"([-+]?\d+(?:\.\d+)?(?:[eE][-+]?\d+)?)\s*\]"
)
def scalar(cur, sql: str, params: Sequence[Any] = ()) -> Any:
cur.execute(sql, params)
row = cur.fetchone()
return next(iter(row.values())) if row else None
def begin_read_only(cur, warnings: List[str]):
try:
cur.execute("SET SESSION TRANSACTION READ ONLY")
except Exception as e:
warnings.append(f"SET SESSION TRANSACTION READ ONLY failed: {e}")
try:
cur.execute("START TRANSACTION READ ONLY")
except Exception as e:
warnings.append(f"START TRANSACTION READ ONLY failed: {e}")
cur.execute("START TRANSACTION")
def to_bytes(v: Any) -> Optional[bytes]:
if v is None:
return None
if isinstance(v, bytes):
return v
if isinstance(v, bytearray):
return bytes(v)
if isinstance(v, memoryview):
return v.tobytes()
try:
return bytes(v)
except Exception:
return None
def write_csv(
path: Path,
rows: Iterable[Dict[str, Any]],
fieldnames: Optional[List[str]] = None,
):
rows = list(rows)
if fieldnames is None:
seen = set()
fieldnames = []
for row in rows:
for k in row.keys():
if k not in seen:
seen.add(k)
fieldnames.append(k)
path.parent.mkdir(parents=True, exist_ok=True)
with path.open("w", encoding="utf-8-sig", newline="") as f:
if not fieldnames:
return
writer = csv.DictWriter(
f,
fieldnames=fieldnames,
extrasaction="ignore"
)
writer.writeheader()
for row in rows:
out = {}
for k in fieldnames:
v = row.get(k)
if isinstance(v, (dict, list, tuple)):
v = json.dumps(v, ensure_ascii=False, default=str)
out[k] = "" if v is None else v
writer.writerow(out)
def probe_wkb_fn(cur) -> str:
for fn in ("ST_AsBinary", "ST_AsWKB", "AsBinary", "AsWKB"):
try:
cur.execute(
f"SELECT {fn}(geom) AS x "
f"FROM chart_area_geometry LIMIT 1"
)
cur.fetchone()
return fn
except Exception:
pass
raise RuntimeError("找不到 geometry -> WKB 函数")
def geometry_hash(g) -> Optional[str]:
if g is None:
return None
try:
return hashlib.sha256(dump_wkb(g)).hexdigest()
except Exception:
return None
def safe_symdiff_area(a, b) -> Optional[float]:
if a is None or b is None:
return None
try:
return float(a.symmetric_difference(b).area)
except Exception:
return None
def safe_hausdorff(a, b) -> Optional[float]:
if a is None or b is None:
return None
try:
return float(a.hausdorff_distance(b))
except Exception:
return None
def safe_length(g) -> float:
if g is None:
return 0.0
try:
return float(g.length)
except Exception:
return 0.0
def flatten_components(g) -> Dict[str, List[Any]]:
result = {
"polygons": [],
"lines": [],
"points": [],
"others": [],
}
def rec(x):
if x is None:
return
t = x.geom_type
if t == "Polygon":
result["polygons"].append(x)
elif t == "MultiPolygon":
for y in x.geoms:
rec(y)
elif t in ("LineString", "LinearRing"):
result["lines"].append(x)
elif t == "MultiLineString":
for y in x.geoms:
rec(y)
elif t == "Point":
result["points"].append(x)
elif t == "MultiPoint":
for y in x.geoms:
rec(y)
elif hasattr(x, "geoms"):
for y in x.geoms:
rec(y)
else:
result["others"].append(x)
rec(g)
return result
def polygon_part_count(g) -> int:
if g is None:
return 0
if g.geom_type == "Polygon":
return 1
if g.geom_type == "MultiPolygon":
return len(g.geoms)
return len(flatten_components(g)["polygons"])
# ---------------------------------------------------------------------------
# Odd-even normalization
# ---------------------------------------------------------------------------
def original_rings(geom) -> List[List[Tuple[float, float]]]:
rings: List[List[Tuple[float, float]]] = []
def add_polygon(poly):
rings.append([
(float(c[0]), float(c[1]))
for c in poly.exterior.coords
])
for ring in poly.interiors:
rings.append([
(float(c[0]), float(c[1]))
for c in ring.coords
])
if geom.geom_type == "Polygon":
add_polygon(geom)
elif geom.geom_type == "MultiPolygon":
for poly in geom.geoms:
add_polygon(poly)
return rings
def point_on_segment(
px: float, py: float,
ax: float, ay: float,
bx: float, by: float,
eps: float = 1e-12,
) -> bool:
vx = bx - ax
vy = by - ay
wx = px - ax
wy = py - ay
c = abs(vx * wy - vy * wx)
if c > eps * max(1.0, abs(vx), abs(vy)):
return False
d = wx * vx + wy * vy
if d < -eps:
return False
vv = vx * vx + vy * vy
if d - vv > eps:
return False
return True
def even_odd_inside(
point: Tuple[float, float],
rings: List[List[Tuple[float, float]]],
) -> bool:
x, y = point
crossings = 0
for ring in rings:
for i in range(len(ring) - 1):
x1, y1 = ring[i]
x2, y2 = ring[i + 1]
if point_on_segment(x, y, x1, y1, x2, y2):
return True
if (y1 > y) != (y2 > y):
xinters = x1 + (y - y1) * (x2 - x1) / (y2 - y1)
if xinters > x:
crossings += 1
return crossings % 2 == 1
def select_faces_odd_even(geom) -> Tuple[List[Any], Dict[str, Any]]:
rings = original_rings(geom)
noded = unary_union(geom.boundary)
polygons_geom, cuts, dangles, invalid_rings = polygonize_full(noded)
faces = flatten_components(polygons_geom)["polygons"]
selected = []
for face in faces:
rp = face.representative_point()
p = (float(rp.x), float(rp.y))
if even_odd_inside(p, rings):
selected.append(face)
return selected, {
"face_count": len(faces),
"selected_face_count": len(selected),
"cuts_length_deg": safe_length(cuts),
"dangles_length_deg": safe_length(dangles),
"invalid_rings_length_deg": safe_length(invalid_rings),
}
def candidate_from_selected_faces(
selected: List[Any],
mode: str = "original",
):
if not selected:
return None
if mode == "original":
items = list(selected)
elif mode == "reverse":
items = list(reversed(selected))
elif mode == "area":
items = sorted(selected, key=lambda g: (float(g.area), g.bounds))
elif mode == "centroid":
items = sorted(
selected,
key=lambda g: (
float(g.centroid.x),
float(g.centroid.y),
float(g.area),
)
)
else:
raise ValueError(mode)
return unary_union(items)
def normalize_odd_even(geom):
if geom.is_valid:
return geom
selected, _ = select_faces_odd_even(geom)
if not selected:
return None
return unary_union(selected)
# ---------------------------------------------------------------------------
# Equivalent input variants
# ---------------------------------------------------------------------------
def rotate_closed_ring(
coords: List[Tuple[float, float]],
shift: int
) -> List[Tuple[float, float]]:
if not coords:
return coords
core = list(coords)
if len(core) > 1 and core[0] == core[-1]:
core = core[:-1]
if not core:
return coords
shift = shift % len(core)
rotated = core[shift:] + core[:shift]
rotated.append(rotated[0])
return rotated
def reverse_closed_ring(
coords: List[Tuple[float, float]]
) -> List[Tuple[float, float]]:
core = list(coords)
if len(core) > 1 and core[0] == core[-1]:
core = core[:-1]
rev = list(reversed(core))
if rev:
rev.append(rev[0])
return rev
def polygon_variant(poly, *, reverse: bool, rotate: bool):
ext = [
(float(c[0]), float(c[1]))
for c in poly.exterior.coords
]
holes = [
[
(float(c[0]), float(c[1]))
for c in ring.coords
]
for ring in poly.interiors
]
def tx(coords, index):
out = coords
if reverse:
out = reverse_closed_ring(out)
if rotate:
core_len = max(1, len(out) - 1)
shift = max(1, core_len // 3 + index)
out = rotate_closed_ring(out, shift)
return out
ext2 = tx(ext, 0)
holes2 = [tx(h, i + 1) for i, h in enumerate(holes)]
return Polygon(ext2, holes2)
def make_input_variant(geom, mode: str):
if mode == "original":
return geom
reverse = mode in {
"reverse_rings",
"reverse_and_rotate",
"reverse_components_and_rings",
}
rotate = mode in {
"rotate_rings",
"reverse_and_rotate",
}
if geom.geom_type == "Polygon":
return polygon_variant(
geom,
reverse=reverse,
rotate=rotate,
)
if geom.geom_type == "MultiPolygon":
polys = [
polygon_variant(
p,
reverse=reverse,
rotate=rotate,
)
for p in geom.geoms
]
if mode in {
"reverse_components",
"reverse_components_and_rings",
}:
polys = list(reversed(polys))
return MultiPolygon(polys)
return geom
# ---------------------------------------------------------------------------
# DB loader
# ---------------------------------------------------------------------------
def load_polygon_records(cur, dataset_id: int, wkb_fn: str):
records = []
cur.execute(f"""
SELECT
'AREA' AS kind,
f.id AS feature_id,
NULL AS coverage_geom_id,
f.object_class,
f.cell_id,
c.cell_name,
f.source_record_index,
f.coverage_id,
f.raw_geotype,
f.bbox_min_lon,
f.bbox_min_lat,
f.bbox_max_lon,
f.bbox_max_lat,
c.bbox_min_lon AS cell_bbox_min_lon,
c.bbox_min_lat AS cell_bbox_min_lat,
c.bbox_max_lon AS cell_bbox_max_lon,
c.bbox_max_lat AS cell_bbox_max_lat,
{wkb_fn}(g.geom) AS geom_wkb
FROM chart_feature f
JOIN chart_area_geometry g
ON g.feature_id=f.id
JOIN chart_cell c
ON c.id=f.cell_id
WHERE f.dataset_id=%s
AND f.geometry_type='AREA'
ORDER BY f.id
""", (dataset_id,))
for r in cur:
try:
geom = load_wkb(to_bytes(r["geom_wkb"]))
x = dict(r)
x["_geom"] = geom
records.append(x)
except Exception:
continue
cur.execute(f"""
SELECT
'COVERAGE' AS kind,
NULL AS feature_id,
v.id AS coverage_geom_id,
'COVERAGE' AS object_class,
v.cell_id,
c.cell_name,
v.source_record_index,
NULL AS coverage_id,
NULL AS raw_geotype,
NULL AS bbox_min_lon,
NULL AS bbox_min_lat,
NULL AS bbox_max_lon,
NULL AS bbox_max_lat,
c.bbox_min_lon AS cell_bbox_min_lon,
c.bbox_min_lat AS cell_bbox_min_lat,
c.bbox_max_lon AS cell_bbox_max_lon,
c.bbox_max_lat AS cell_bbox_max_lat,
{wkb_fn}(v.geom) AS geom_wkb
FROM chart_coverage v
JOIN chart_cell c
ON c.id=v.cell_id
WHERE v.dataset_id=%s
ORDER BY v.id
""", (dataset_id,))
for r in cur:
try:
geom = load_wkb(to_bytes(r["geom_wkb"]))
x = dict(r)
x["_geom"] = geom
records.append(x)
except Exception:
continue
return records
# ---------------------------------------------------------------------------
# GeoJSON
# ---------------------------------------------------------------------------
def make_geojson(records, path: Path):
features = []
def add(g, props):
if g is None:
return
try:
features.append({
"type": "Feature",
"geometry": mapping(g),
"properties": props,
})
except Exception:
pass
for r in records:
if r["_geom"].is_valid:
continue
rid = (
r.get("feature_id")
if r["kind"] == "AREA"
else r.get("coverage_geom_id")
)
common = {
"kind": r["kind"],
"record_id": rid,
"object_class": r.get("object_class"),
"cell_name": r.get("cell_name"),
"canonical_hash": r.get("_canonical_hash"),
}
add(
r["_geom"],
{
**common,
"layer": "original_invalid",
}
)
add(
r.get("_canonical_candidate"),
{
**common,
"layer": "canonical_candidate",
}
)
path.write_text(
json.dumps(
{
"type": "FeatureCollection",
"features": features,
},
ensure_ascii=False,
separators=(",", ":"),
default=str,
),
encoding="utf-8"
)
def parse_args():
p = argparse.ArgumentParser(
description="CM93 normalization canonical validation v10"
)
p.add_argument("--host", default=os.getenv("MYSQL_HOST", "127.0.0.1"))
p.add_argument("--port", type=int, default=int(os.getenv("MYSQL_PORT", "3306")))
p.add_argument("--user", default=os.getenv("MYSQL_USER", "root"))
p.add_argument("--password", default=os.getenv("MYSQL_PASSWORD"))
p.add_argument("--database", default=os.getenv("MYSQL_DATABASE", "cm93_chart"))
p.add_argument("--dataset-id", type=int, default=1)
p.add_argument(
"--topology-area-tolerance",
type=float,
default=1e-12,
)
p.add_argument(
"--hausdorff-tolerance",
type=float,
default=1e-12,
)
p.add_argument(
"--bbox-tolerance",
type=float,
default=5e-6,
)
p.add_argument("--output", default=None)
return p.parse_args()
def main():
args = parse_args()
if not args.password:
args.password = getpass.getpass("MySQL password: ")
stamp = datetime.now().strftime("%Y%m%d_%H%M%S")
outdir = Path(
args.output
or f"cm93_normalization_canonical_v10_dataset_{args.dataset_id}_{stamp}"
)
outdir.mkdir(parents=True, exist_ok=True)
report = {
"version": "v10.0",
"generated_at": datetime.now().isoformat(timespec="seconds"),
"database": args.database,
"dataset_id": args.dataset_id,
"read_only": True,
"pbf": "OUT_OF_SCOPE",
"shapely_version": getattr(shapely, "__version__", "unknown"),
"warnings": [],
}
conn = None
try:
conn = pymysql.connect(
host=args.host,
port=args.port,
user=args.user,
password=args.password,
database=args.database,
charset="utf8mb4",
cursorclass=DictCursor,
autocommit=False,
connect_timeout=20,
read_timeout=600,
write_timeout=30,
)
with conn.cursor() as cur:
begin_read_only(cur, report["warnings"])
report["mysql_version"] = scalar(cur, "SELECT VERSION()")
wkb_fn = probe_wkb_fn(cur)
print("[1/6] 加载 Area + Coverage ...")
records = load_polygon_records(
cur,
args.dataset_id,
wkb_fn
)
valid = [r for r in records if r["_geom"].is_valid]
invalid = [r for r in records if not r["_geom"].is_valid]
print(
f" total={len(records):,} "
f"valid={len(valid):,} invalid={len(invalid):,}"
)
print("[2/6] Valid PASS-THROUGH policy guard ...")
valid_policy_changed = 0
valid_hash_changed = 0
for i, r in enumerate(valid, start=1):
# policy: valid geometry 不进入 canonical rewrite
out = r["_geom"]
if out is not r["_geom"]:
valid_policy_changed += 1
if geometry_hash(out) != geometry_hash(r["_geom"]):
valid_hash_changed += 1
if i % 20000 == 0:
print(f" valid checked={i:,}")
print("[3/6] 55 invalid canonical / input-variant 验证 ...")
golden_rows = []
variant_rows = []
face_order_rows = []
failure_rows = []
bbox_rows = []
input_modes = [
"original",
"reverse_rings",
"rotate_rings",
"reverse_and_rotate",
"reverse_components",
"reverse_components_and_rings",
]
for i, r in enumerate(invalid, start=1):
geom = r["_geom"]
rid = (
r.get("feature_id")
if r["kind"] == "AREA"
else r.get("coverage_geom_id")
)
selected, sel_info = select_faces_odd_even(geom)
baseline = candidate_from_selected_faces(
selected,
"original"
)
if baseline is None:
fail = {
"kind": r["kind"],
"record_id": rid,
"object_class": r.get("object_class"),
"cell_name": r.get("cell_name"),
"failure": "BASELINE_CANDIDATE_NONE",
}
failure_rows.append(fail)
continue
baseline_canonical = canonicalize(baseline)
baseline_hash = geometry_hash(baseline_canonical)
baseline_valid = bool(baseline_canonical.is_valid)
baseline_polygon_only = (
baseline_canonical.geom_type
in {"Polygon", "MultiPolygon"}
)
# make_valid line remainder length only as metadata
mv = make_valid_geom(geom)
mvcomps = flatten_components(mv)
remainder = geom_union(mvcomps["lines"])
remainder_length = safe_length(remainder)
# canonicalization itself topology invariant
canonical_sym = safe_symdiff_area(
baseline,
baseline_canonical
)
canonical_haus = safe_hausdorff(
baseline,
baseline_canonical
)
canonical_topology_ok = (
canonical_sym is not None
and canonical_sym <= args.topology_area_tolerance
and canonical_haus is not None
and canonical_haus <= args.hausdorff_tolerance
)
# face order invariance
face_hashes = {}
for mode in [
"original",
"reverse",
"area",
"centroid",
]:
cand = candidate_from_selected_faces(
selected,
mode
)
canon = canonicalize(cand)
h = geometry_hash(canon)
face_hashes[mode] = h
face_order_rows.append({
"kind": r["kind"],
"feature_id": r.get("feature_id"),
"coverage_geom_id": r.get("coverage_geom_id"),
"object_class": r.get("object_class"),
"cell_name": r.get("cell_name"),
"face_order_mode": mode,
"canonical_wkb_sha256": h,
"matches_baseline_hash":
int(h == baseline_hash),
"topology_symdiff_area_deg2":
safe_symdiff_area(
baseline_canonical,
canon
),
})
face_order_stable = (
len(set(face_hashes.values())) == 1
)
# repeated same input
repeat_hashes = []
for _ in range(3):
c = normalize_odd_even(geom)
cc = canonicalize(c)
repeat_hashes.append(geometry_hash(cc))
repeat_stable = (
len(set(repeat_hashes)) == 1
and repeat_hashes[0] == baseline_hash
)
# input-equivalent variants
variant_fail = 0
for mode in input_modes:
try:
vg = make_input_variant(geom, mode)
vc = normalize_odd_even(vg)
if vc is None:
variant_fail += 1
variant_rows.append({
"kind": r["kind"],
"feature_id": r.get("feature_id"),
"coverage_geom_id": r.get("coverage_geom_id"),
"object_class": r.get("object_class"),
"cell_name": r.get("cell_name"),
"variant_mode": mode,
"status": "FAIL_NO_CANDIDATE",
})
continue
vcanon = canonicalize(vc)
vh = geometry_hash(vcanon)
sym = safe_symdiff_area(
baseline_canonical,
vcanon
)
haus = safe_hausdorff(
baseline_canonical,
vcanon
)
topology_equal = (
bool(
baseline_canonical.equals(
vcanon
)
)
and sym is not None
and sym <= args.topology_area_tolerance
)
hash_equal = (vh == baseline_hash)
if not (
vcanon.is_valid
and topology_equal
and hash_equal
):
variant_fail += 1
variant_rows.append({
"kind": r["kind"],
"feature_id": r.get("feature_id"),
"coverage_geom_id": r.get("coverage_geom_id"),
"object_class": r.get("object_class"),
"cell_name": r.get("cell_name"),
"variant_mode": mode,
"variant_input_valid":
int(bool(vg.is_valid)),
"candidate_valid":
int(bool(vcanon.is_valid)),
"candidate_type":
vcanon.geom_type,
"topology_equal_baseline":
int(topology_equal),
"symdiff_area_deg2": sym,
"hausdorff_distance_deg": haus,
"canonical_wkb_sha256": vh,
"matches_baseline_hash":
int(hash_equal),
"status":
"PASS"
if (
vcanon.is_valid
and topology_equal
and hash_equal
)
else "FAIL",
})
except Exception as e:
variant_fail += 1
variant_rows.append({
"kind": r["kind"],
"feature_id": r.get("feature_id"),
"coverage_geom_id": r.get("coverage_geom_id"),
"object_class": r.get("object_class"),
"cell_name": r.get("cell_name"),
"variant_mode": mode,
"status": "EXCEPTION",
"detail":
f"{type(e).__name__}: {e}",
})
# bbox metadata impact
old_bounds = tuple(float(v) for v in geom.bounds)
new_bounds = tuple(
float(v) for v in baseline_canonical.bounds
)
bbox_changed = any(
abs(a - b) > args.bbox_tolerance
for a, b in zip(old_bounds, new_bounds)
)
feature_bbox = None
feature_bbox_matches_candidate = None
if r["kind"] == "AREA":
vals = [
r.get("bbox_min_lon"),
r.get("bbox_min_lat"),
r.get("bbox_max_lon"),
r.get("bbox_max_lat"),
]
if all(v is not None for v in vals):
feature_bbox = tuple(
float(v) for v in vals
)
feature_bbox_matches_candidate = all(
abs(a - b) <= args.bbox_tolerance
for a, b in zip(
feature_bbox,
new_bounds
)
)
bbox_rows.append({
"kind": r["kind"],
"feature_id": r.get("feature_id"),
"coverage_geom_id": r.get("coverage_geom_id"),
"object_class": r.get("object_class"),
"cell_name": r.get("cell_name"),
"original_geometry_bounds": old_bounds,
"candidate_bounds": new_bounds,
"bbox_changed_over_tolerance":
int(bbox_changed),
"current_feature_bbox": feature_bbox,
"current_feature_bbox_matches_candidate":
(
int(feature_bbox_matches_candidate)
if feature_bbox_matches_candidate
is not None
else None
),
"future_feature_bbox_update_needed":
(
int(not feature_bbox_matches_candidate)
if feature_bbox_matches_candidate
is not None
else None
),
"discarded_connector_length_deg":
remainder_length,
})
record_pass = all([
baseline_valid,
baseline_polygon_only,
canonical_topology_ok,
face_order_stable,
repeat_stable,
variant_fail == 0,
])
golden = {
"kind": r["kind"],
"feature_id": r.get("feature_id"),
"coverage_geom_id": r.get("coverage_geom_id"),
"object_class": r.get("object_class"),
"cell_id": r.get("cell_id"),
"cell_name": r.get("cell_name"),
"source_record_index":
r.get("source_record_index"),
"coverage_id": r.get("coverage_id"),
"raw_geotype": r.get("raw_geotype"),
"normalization_rule":
"ODD_EVEN_FACE_SELECTION",
"canonicalization_rule":
"GEOS_NORMALIZE",
"candidate_type":
baseline_canonical.geom_type,
"candidate_polygon_part_count":
polygon_part_count(
baseline_canonical
),
"candidate_valid":
int(baseline_valid),
"polygonize_face_count":
sel_info["face_count"],
"selected_face_count":
sel_info["selected_face_count"],
"candidate_area_deg2":
float(baseline_canonical.area),
"candidate_bbox":
tuple(
float(v)
for v in baseline_canonical.bounds
),
"discarded_connector_length_deg":
remainder_length,
"canonical_wkb_sha256":
baseline_hash,
"canonicalization_symdiff_area_deg2":
canonical_sym,
"canonicalization_hausdorff_deg":
canonical_haus,
"canonical_topology_unchanged":
int(canonical_topology_ok),
"face_order_hash_stable":
int(face_order_stable),
"repeat_hash_stable":
int(repeat_stable),
"input_variant_failure_count":
variant_fail,
"record_status":
"PASS"
if record_pass
else "FAIL",
}
golden_rows.append(golden)
if not record_pass:
failure_rows.append(golden)
r["_canonical_candidate"] = baseline_canonical
r["_canonical_hash"] = baseline_hash
print(
f" [{i:02d}/{len(invalid)}] "
f"{r['kind']} id={rid} "
f"{r.get('object_class')} -> "
f"{'PASS' if record_pass else 'FAIL'}"
)
print("[4/6] 汇总 canonical stability ...")
record_status = Counter(
r["record_status"]
for r in golden_rows
)
variant_status = Counter(
r.get("status", "UNKNOWN")
for r in variant_rows
)
bbox_update_area_count = sum(
int(
r.get(
"future_feature_bbox_update_needed"
) or 0
)
for r in bbox_rows
)
report["summary"] = {
"total_polygon_records": len(records),
"valid_polygon_records": len(valid),
"invalid_polygon_records": len(invalid),
"valid_policy_changed_count":
valid_policy_changed,
"valid_hash_changed_count":
valid_hash_changed,
"golden_vector_count":
len(golden_rows),
"golden_record_status_counts":
dict(record_status),
"input_variant_case_count":
len(variant_rows),
"input_variant_status_counts":
dict(variant_status),
"face_order_case_count":
len(face_order_rows),
"face_order_hash_mismatch_count":
sum(
not bool(
r["matches_baseline_hash"]
)
for r in face_order_rows
),
"canonical_topology_failure_count":
sum(
not bool(
r["canonical_topology_unchanged"]
)
for r in golden_rows
),
"repeat_hash_failure_count":
sum(
not bool(
r["repeat_hash_stable"]
)
for r in golden_rows
),
"input_variant_failure_record_count":
sum(
int(
r["input_variant_failure_count"]
) > 0
for r in golden_rows
),
"future_area_feature_bbox_update_needed_count":
bbox_update_area_count,
}
overall_pass = (
len(invalid) == 55
and len(golden_rows) == 55
and record_status.get("PASS", 0) == 55
and record_status.get("FAIL", 0) == 0
and valid_policy_changed == 0
and valid_hash_changed == 0
and report["summary"][
"face_order_hash_mismatch_count"
] == 0
and report["summary"][
"canonical_topology_failure_count"
] == 0
and report["summary"][
"repeat_hash_failure_count"
] == 0
and report["summary"][
"input_variant_failure_record_count"
] == 0
)
report["summary"]["overall_pass"] = overall_pass
print("[5/6] 写 golden vectors / 报告 ...")
write_csv(
outdir / "normalization_golden_vectors_v10.csv",
golden_rows
)
(
outdir / "normalization_golden_vectors_v10.json"
).write_text(
json.dumps(
golden_rows,
ensure_ascii=False,
indent=2,
default=str,
),
encoding="utf-8"
)
write_csv(
outdir / "normalization_variant_results.csv",
variant_rows
)
write_csv(
outdir / "normalization_face_order_results.csv",
face_order_rows
)
write_csv(
outdir / "normalization_canonical_failures.csv",
failure_rows
)
write_csv(
outdir / "normalization_bbox_metadata_impact.csv",
bbox_rows
)
make_geojson(
records,
outdir / "normalization_canonical_candidates.geojson"
)
(
outdir / "normalization_canonical_report_v10.json"
).write_text(
json.dumps(
report,
ensure_ascii=False,
indent=2,
default=str,
),
encoding="utf-8"
)
lines = []
lines.append(
"CM93 -> MySQL OGC Normalization Canonical Validation v10"
)
lines.append("=" * 96)
lines.append(
f"Generated: {report['generated_at']}"
)
lines.append(
f"Dataset ID: {args.dataset_id}"
)
lines.append(
f"MySQL: {report['mysql_version']}"
)
lines.append(
f"Shapely: {report['shapely_version']}"
)
lines.append("")
lines.append("SCOPE")
lines.append(
" CM93JAPAN 是局部 CM93 数据集;"
"跨未导入海域不连续不作为错误。"
)
lines.append(
" V10 只验证 odd-even candidate 的 canonical WKB "
"可重复性,不修改数据库。"
)
lines.append("")
lines.append("BASELINE")
lines.append(
f" polygon records: {len(records)}"
)
lines.append(
f" valid: {len(valid)}"
)
lines.append(
f" invalid: {len(invalid)}"
)
lines.append("")
lines.append("VALID PASS-THROUGH")
lines.append(
f" policy changed: {valid_policy_changed}"
)
lines.append(
f" WKB hash changed: {valid_hash_changed}"
)
lines.append("")
lines.append("GOLDEN VECTOR STATUS")
for k, v in record_status.most_common():
lines.append(
f" {k:24s} {v}"
)
lines.append("")
lines.append("INPUT VARIANTS")
lines.append(
f" cases: {len(variant_rows)}"
)
for k, v in variant_status.most_common():
lines.append(
f" {k:24s} {v}"
)
lines.append("")
lines.append("STABILITY")
lines.append(
f" face-order hash mismatch: "
f"{report['summary']['face_order_hash_mismatch_count']}"
)
lines.append(
f" canonical topology failures: "
f"{report['summary']['canonical_topology_failure_count']}"
)
lines.append(
f" repeat hash failures: "
f"{report['summary']['repeat_hash_failure_count']}"
)
lines.append(
f" input-variant failure records: "
f"{report['summary']['input_variant_failure_record_count']}"
)
lines.append("")
lines.append("METADATA IMPACT")
lines.append(
f" future Area feature bbox updates needed: "
f"{bbox_update_area_count}"
)
lines.append(
" bbox 更新仅是未来若真正保存 normalized geometry "
"时的元数据同步要求;V10 不修改任何 bbox。"
)
lines.append("")
lines.append("OVERALL")
lines.append(
f" CANONICAL VALIDATION PASS: "
f"{'YES' if overall_pass else 'NO'}"
)
lines.append("")
lines.append("NEXT USE")
lines.append(
" normalization_golden_vectors_v10.csv/json "
"可作为未来 importer normalization 实现的 55 条 golden regression vectors。"
)
lines.append(
" 后续实现应固定:ODD_EVEN_FACE_SELECTION + GEOS_NORMALIZE,"
"并要求输出 canonical_wkb_sha256 与 golden vector 一致。"
)
lines.append("")
lines.append(
"NO DATABASE WRITE / NO IMPORTER CHANGE / NO PBF"
)
(
outdir / "normalization_canonical_summary_v10.txt"
).write_text(
"\n".join(lines),
encoding="utf-8"
)
conn.rollback()
print("[6/6] 压缩 ...")
zpath = outdir.with_suffix(".zip")
with zipfile.ZipFile(
zpath,
"w",
compression=zipfile.ZIP_DEFLATED
) as z:
for p in sorted(outdir.rglob("*")):
if p.is_file():
z.write(
p,
p.relative_to(outdir.parent)
)
print("")
print("完成。")
print(f"Output: {outdir.resolve()}")
print(f"ZIP: {zpath.resolve()}")
print("把 ZIP 整包上传给 ChatGPT。")
except KeyboardInterrupt:
if conn:
try:
conn.rollback()
except Exception:
pass
print("\n用户中断。")
sys.exit(130)
except Exception:
if conn:
try:
conn.rollback()
except Exception:
pass
tb = traceback.format_exc()
print(tb)
try:
(
outdir / "fatal_error.txt"
).write_text(
tb,
encoding="utf-8"
)
except Exception:
pass
sys.exit(1)
finally:
if conn:
try:
conn.close()
except Exception:
pass
if __name__ == "__main__":
main()