#!/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()