#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ CM93 -> MySQL Invalid Polygon 线性残片与填充规则诊断 Version: v8.0 背景 ==== 前序 V7 已确认： - invalid total = 55（Area 52 + Coverage 3） - 41 条 make_valid 后得到 Polygon + LineString - 14 条得到纯 Polygon/MultiPolygon - 其中 5 条 make_valid 面选择与“polygonize 全部 faces”不同 - 31 条 polygonize 有 dangles，10 条有 cuts，恰好合计 41 V8 目标 ======= 只做只读、内存诊断，不修复数据库。 A. 验证 41 条线性残片的来源 -------------------------------- 对 make_valid 产生的每一条线性 remainder： 1. 将 remainder 拆成最小线段； 2. 对每个 remainder segment 的中点，反查原始 invalid Polygon 的所有 ring segment； 3. 统计该点被多少条原始 segment 覆盖； 4. 判断覆盖 segment 相对于 remainder segment 的方向： - 正向 - 反向 5. 分类： - OPPOSITE_DIRECTION_RETRACE 同一位置至少有两条原始 segment，且存在正反两个方向； 这是“往返桥接/重走线”的强证据。 - SAME_DIRECTION_OVERLAP - SINGLE_TRAVERSAL - MULTI_COVER_NO_OPPOSITE - UNMATCHED 6. 对每个对象给出 record-level 分类： - RETRACED_CONNECTOR_CONFIRMED - RETRACED_CONNECTOR_MIXED - NONRETRACED_LINEAR_REMAINDER - NO_LINEAR_REMAINDER B. 验证 polygon face 填充规则 -------------------------------- 对原始 invalid geometry.boundary： 1. unary_union + polygonize_full，得到所有 planar faces； 2. 对每个 face 取 representative point； 3. 对原始 ring 坐标计算： - Odd-even parity - Non-zero winding number 4. 检查该 face 是否被 Shapely make_valid 的 polygonal result 选中； 5. 每条记录判定： - MATCHES_BOTH_ODD_EVEN_AND_NONZERO - MATCHES_ODD_EVEN_ONLY - MATCHES_NONZERO_ONLY - MATCHES_NEITHER - NO_FACES 这一步尤其解释 V7 中 5 条 POLYGON_ONLY_FACE_SELECTION_DIFFERS： 不是坐标变了，而是“全部 polygonized faces 中哪些应该填充”存在语义选择。 C. 输出结论候选 -------------------------------- 仅用于诊断，不自动修改数据库： - OGC_AREA_WITH_RETRACED_CONNECTOR 面部分可有效表达，线性残片基本确认是原 ring 的正反重走桥接线。 - OGC_AREA_WITH_NONRETRACED_REMAINDER 存在线性残片，但无法解释为简单正反重走，需要继续调查。 - OGC_POLYGONAL_FILL_RULE_RESOLVED 无线性残片，且 make_valid face selection 能被 odd-even/nonzero 规则解释。 - OGC_POLYGONAL_FILL_RULE_UNRESOLVED face selection 无法由两种常见规则解释。 重要边界 ======== 当前 CM93JAPAN 是日本沿海局部子集，不是完整 CM93： - 跨未导入海域连续性不作为 FAIL； - Cell 边界上的共享拓扑保留 boundary context； - V8 不将“Cell 边界 self-touch”自动判为源数据损坏。 禁止 ==== - 不 UPDATE / INSERT / DELETE - 不写回 make_valid - 不 buffer(0) - 不修改 geometry_issue - 不开始 PBF 依赖 ==== pip install pymysql shapely 执行 ==== python cm93_mysql_fillrule_audit_v8.py ^ --host 127.0.0.1 ^ --port 3306 ^ --user root ^ --database cm93_chart ^ --dataset-id 1 输出 ==== - fillrule_audit_summary_v8.txt - fillrule_audit_report_v8.json - invalid_polygon_fillrule_classification.csv - linear_remainder_segment_evidence.csv - polygon_face_fillrule_evidence.csv - fillrule_by_object_class.csv - fillrule_by_boundary_context.csv - retraced_connector_summary.csv - fillrule_diagnostic.geojson 并自动打 ZIP。 """ from __future__ import annotations import argparse import csv import getpass import json import math import os import re import sys import traceback import zipfile from collections import Counter, defaultdict 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 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, LineString, mapping from shapely.ops import unary_union, polygonize_full from shapely.validation import explain_validity except ImportError: print("缺少 shapely：pip install pymysql 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 r in rows: for k in r.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 w = csv.DictWriter(f, fieldnames=fieldnames, extrasaction="ignore") w.writeheader() for r in rows: out = {} for k in fieldnames: v = r.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 w.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 FROM chart_area_geometry LIMIT 1") cur.fetchone() return fn except Exception: pass raise RuntimeError("找不到 geometry -> WKB 函数") def parse_invalid_point(reason: str) -> Optional[Tuple[float, float]]: if not reason: return None m = INVALID_POINT_RE.search(reason) if not m: return None try: return float(m.group(1)), float(m.group(2)) except Exception: return None def bbox_from_row(row: Dict[str, Any]) -> Optional[Tuple[float, float, float, float]]: try: vals = ( float(row["cell_bbox_min_lon"]), float(row["cell_bbox_min_lat"]), float(row["cell_bbox_max_lon"]), float(row["cell_bbox_max_lat"]), ) if all(math.isfinite(v) for v in vals): return vals except Exception: pass return None def invalid_point_cell_edge( p: Optional[Tuple[float, float]], bbox: Optional[Tuple[float, float, float, float]], tol: float ) -> Tuple[int, str, Optional[float]]: if p is None or bbox is None: return 0, "", None x, y = p minx, miny, maxx, maxy = bbox d = { "WEST": abs(x - minx), "EAST": abs(x - maxx), "SOUTH": abs(y - miny), "NORTH": abs(y - maxy), } edge, dist = min(d.items(), key=lambda kv: kv[1]) return int(dist <= tol), edge if dist <= tol else "", dist 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 geom_union(items: Sequence[Any]) -> Optional[Any]: if not items: return None try: return unary_union(list(items)) except Exception: return None 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 p in geom.geoms: add_polygon(p) return rings def original_segments(geom) -> List[Dict[str, Any]]: segs = [] for ri, ring in enumerate(original_rings(geom)): for si in range(len(ring) - 1): a = ring[si] b = ring[si + 1] if a == b: continue try: line = LineString([a, b]) except Exception: continue segs.append({ "ring_index": ri, "segment_index": si, "a": a, "b": b, "line": line, }) return segs def vector_direction(a: Tuple[float, float], b: Tuple[float, float]) -> Tuple[float, float]: return b[0] - a[0], b[1] - a[1] def dot(a: Tuple[float, float], b: Tuple[float, float]) -> float: return a[0] * b[0] + a[1] * b[1] def cross(a: Tuple[float, float], b: Tuple[float, float]) -> float: return a[0] * b[1] - a[1] * b[0] def classify_remainder_segment( ra: Tuple[float, float], rb: Tuple[float, float], source_segments: List[Dict[str, Any]], tol: float, ) -> Dict[str, Any]: rvec = vector_direction(ra, rb) rlen = math.hypot(rvec[0], rvec[1]) if rlen <= 0: return { "coverage_count": 0, "same_direction_count": 0, "opposite_direction_count": 0, "segment_class": "ZERO_LENGTH_REMAINDER_SEGMENT", } mid = Point((ra[0] + rb[0]) / 2.0, (ra[1] + rb[1]) / 2.0) same = 0 opposite = 0 covers = 0 source_hits = [] for s in source_segments: line = s["line"] # midpoint 必须落在线段上或非常接近 try: if line.distance(mid) > tol: continue except Exception: continue svec = vector_direction(s["a"], s["b"]) slen = math.hypot(svec[0], svec[1]) if slen <= 0: continue # 要求方向基本共线，避免“恰好在交点处”的非同线 segment 被误计。 sin_angle = abs(cross(rvec, svec)) / (rlen * slen) if sin_angle > 1e-7: continue covers += 1 dp = dot(rvec, svec) if dp >= 0: same += 1 direction = "SAME" else: opposite += 1 direction = "OPPOSITE" source_hits.append({ "ring_index": s["ring_index"], "segment_index": s["segment_index"], "direction": direction, }) if covers >= 2 and same > 0 and opposite > 0: cls = "OPPOSITE_DIRECTION_RETRACE" elif covers >= 2 and opposite == 0: cls = "SAME_DIRECTION_OVERLAP" elif covers == 1: cls = "SINGLE_TRAVERSAL" elif covers >= 2: cls = "MULTI_COVER_NO_OPPOSITE" else: cls = "UNMATCHED" return { "coverage_count": covers, "same_direction_count": same, "opposite_direction_count": opposite, "segment_class": cls, "source_hits": source_hits, } def record_remainder_class(segment_classes: Counter) -> str: total = sum(segment_classes.values()) if total == 0: return "NO_LINEAR_REMAINDER" retraced = segment_classes.get("OPPOSITE_DIRECTION_RETRACE", 0) unmatched = segment_classes.get("UNMATCHED", 0) single = segment_classes.get("SINGLE_TRAVERSAL", 0) same = segment_classes.get("SAME_DIRECTION_OVERLAP", 0) other = total - retraced - unmatched - single - same if retraced == total: return "RETRACED_CONNECTOR_CONFIRMED" if retraced > 0 and (single > 0 or same > 0 or unmatched > 0 or other > 0): return "RETRACED_CONNECTOR_MIXED" return "NONRETRACED_LINEAR_REMAINDER" # --------------------------------------------------------------------------- # Fill rules # --------------------------------------------------------------------------- def point_on_segment(px, py, ax, ay, bx, by, eps=1e-12) -> bool: vx, vy = bx - ax, by - ay wx, wy = px - ax, 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): # representative point 理论上不会落边界；若发生，按 inside 处理。 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 is_left(ax, ay, bx, by, px, py) -> float: return (bx - ax) * (py - ay) - (px - ax) * (by - ay) def winding_number(point: Tuple[float, float], rings: List[List[Tuple[float, float]]]) -> int: x, y = point wn = 0 for ring in rings: for i in range(len(ring) - 1): x1, y1 = ring[i] x2, y2 = ring[i + 1] if y1 <= y: if y2 > y and is_left(x1, y1, x2, y2, x, y) > 0: wn += 1 else: if y2 <= y and is_left(x1, y1, x2, y2, x, y) < 0: wn -= 1 return wn def fill_rule_record_class(face_rows: List[Dict[str, Any]]) -> str: if not face_rows: return "NO_FACES" odd_ok = all( bool(r["selected_by_make_valid"]) == bool(r["odd_even_inside"]) for r in face_rows ) nonzero_ok = all( bool(r["selected_by_make_valid"]) == bool(r["nonzero_inside"]) for r in face_rows ) if odd_ok and nonzero_ok: return "MATCHES_BOTH_ODD_EVEN_AND_NONZERO" if odd_ok: return "MATCHES_ODD_EVEN_ONLY" if nonzero_ok: return "MATCHES_NONZERO_ONLY" return "MATCHES_NEITHER" def final_semantic_class(remainder_class: str, fill_class: str) -> str: if remainder_class == "RETRACED_CONNECTOR_CONFIRMED": return "OGC_AREA_WITH_RETRACED_CONNECTOR" if remainder_class in { "RETRACED_CONNECTOR_MIXED", "NONRETRACED_LINEAR_REMAINDER", }: return "OGC_AREA_WITH_NONRETRACED_REMAINDER" if fill_class != "MATCHES_NEITHER": return "OGC_POLYGONAL_FILL_RULE_RESOLVED" return "OGC_POLYGONAL_FILL_RULE_UNRESOLVED" # --------------------------------------------------------------------------- # DB # --------------------------------------------------------------------------- def load_invalid_records(cur, dataset_id: int, wkb_fn: str) -> List[Dict[str, Any]]: 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, 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"])) if not geom.is_valid: 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, 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"])) if not geom.is_valid: x = dict(r) x["_geom"] = geom records.append(x) except Exception: continue return records def make_geojson(records: List[Dict[str, Any]], 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: 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"), "remainder_class": r.get("_remainder_class"), "fill_rule_class": r.get("_fill_rule_class"), "final_semantic_class": r.get("_final_semantic_class"), } add(r.get("_geom"), {**common, "layer": "original_invalid"}) add(r.get("_mv_polygon_union"), {**common, "layer": "make_valid_polygonal"}) add(r.get("_mv_line_union"), {**common, "layer": "linear_remainder"}) p = r.get("_invalid_point") if p: add(Point(p), {**common, "layer": "invalidity_point"}) path.write_text( json.dumps( {"type": "FeatureCollection", "features": features}, ensure_ascii=False, separators=(",", ":"), default=str, ), encoding="utf-8" ) def zip_dir(folder: Path) -> Path: zpath = folder.with_suffix(".zip") with zipfile.ZipFile(zpath, "w", compression=zipfile.ZIP_DEFLATED) as z: for p in sorted(folder.rglob("*")): if p.is_file(): z.write(p, p.relative_to(folder.parent)) return zpath def parse_args(): p = argparse.ArgumentParser( description="CM93 invalid Polygon fill rule + retrace audit v8" ) 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( "--segment-tolerance", type=float, default=2e-8, help="remainder segment midpoint 与原始 ring segment 匹配容差（degree）" ) p.add_argument( "--cell-edge-tolerance", type=float, default=5e-6, help="invalid point 是否位于 Cell bbox 边界的容差（degree）" ) 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_fillrule_audit_v8_dataset_{args.dataset_id}_{stamp}" ) outdir.mkdir(parents=True, exist_ok=True) report = { "version": "v8.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"), "scope_note": ( "CM93JAPAN 是局部数据集。V8 只验证 invalid Polygon 的线性残片" "是否为重走桥接线，以及 face selection 对 odd-even/nonzero 填充规则的符合程度。" ), "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) report["wkb_fn"] = wkb_fn print("[1/5] 加载 invalid Area + Coverage ...") records = load_invalid_records(cur, args.dataset_id, wkb_fn) print(f" invalid records={len(records)}") classification_rows = [] remainder_evidence = [] face_evidence = [] print("[2/5] 线性残片 retrace 诊断 ...") for idx, r in enumerate(records, start=1): geom = r["_geom"] reason = explain_validity(geom) invalid_point = parse_invalid_point(reason) cell_bbox = bbox_from_row(r) on_edge, edge_name, edge_dist = invalid_point_cell_edge( invalid_point, cell_bbox, args.cell_edge_tolerance ) mv = make_valid_geom(geom) comps = flatten_components(mv) mv_poly_union = geom_union(comps["polygons"]) mv_line_union = geom_union(comps["lines"]) src_segments = original_segments(geom) seg_class_counts = Counter() segment_count = 0 for li, line in enumerate(comps["lines"]): coords = [(float(c[0]), float(c[1])) for c in line.coords] for si in range(len(coords) - 1): ra, rb = coords[si], coords[si+1] if ra == rb: continue ev = classify_remainder_segment( ra, rb, src_segments, args.segment_tolerance ) segment_count += 1 seg_class_counts[ev["segment_class"]] += 1 remainder_evidence.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"), "line_component_index": li, "remainder_segment_index": si, "ax": ra[0], "ay": ra[1], "bx": rb[0], "by": rb[1], "segment_length_deg": LineString([ra, rb]).length, "coverage_count": ev["coverage_count"], "same_direction_count": ev["same_direction_count"], "opposite_direction_count": ev["opposite_direction_count"], "segment_class": ev["segment_class"], "source_hits": ev.get("source_hits", []), }) rem_class = record_remainder_class(seg_class_counts) # Face analysis boundary = unary_union(geom.boundary) polygons_geom, cuts, dangles, invalid_rings = polygonize_full(boundary) face_comps = flatten_components(polygons_geom)["polygons"] rings = original_rings(geom) this_faces = [] for fi, face in enumerate(face_comps): rp = face.representative_point() pxy = (float(rp.x), float(rp.y)) odd = even_odd_inside(pxy, rings) wn = winding_number(pxy, rings) nonzero = wn != 0 selected = False if mv_poly_union is not None: try: selected = bool(mv_poly_union.covers(rp)) except Exception: selected = False row = { "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_index": fi, "face_area_deg2": float(face.area), "representative_x": pxy[0], "representative_y": pxy[1], "selected_by_make_valid": int(selected), "odd_even_inside": int(odd), "winding_number": wn, "nonzero_inside": int(nonzero), "odd_even_matches_selection": int(selected == odd), "nonzero_matches_selection": int(selected == nonzero), } face_evidence.append(row) this_faces.append(row) fill_class = fill_rule_record_class(this_faces) final_class = final_semantic_class(rem_class, fill_class) rid = ( r.get("feature_id") if r["kind"] == "AREA" else r.get("coverage_geom_id") ) c = { "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"), "validity_reason": reason, "invalid_point_on_cell_edge": on_edge, "cell_edge": edge_name, "cell_edge_distance_deg": edge_dist, "make_valid_polygon_count": len(comps["polygons"]), "make_valid_line_count": len(comps["lines"]), "linear_remainder_segment_count": segment_count, "remainder_segment_classes": dict(seg_class_counts), "remainder_class": rem_class, "polygonize_face_count": len(this_faces), "faces_selected_by_make_valid": sum( int(x["selected_by_make_valid"]) for x in this_faces ), "faces_odd_even_inside": sum( int(x["odd_even_inside"]) for x in this_faces ), "faces_nonzero_inside": sum( int(x["nonzero_inside"]) for x in this_faces ), "fill_rule_class": fill_class, "final_semantic_class": final_class, } classification_rows.append(c) r["_invalid_point"] = invalid_point r["_mv_polygon_union"] = mv_poly_union r["_mv_line_union"] = mv_line_union r["_remainder_class"] = rem_class r["_fill_rule_class"] = fill_class r["_final_semantic_class"] = final_class print( f" [{idx:02d}/{len(records)}] {r['kind']} id={rid} " f"{r.get('object_class')} -> " f"{rem_class} / {fill_class}" ) print("[3/5] 汇总 ...") remainder_counts = Counter(r["remainder_class"] for r in classification_rows) fill_counts = Counter(r["fill_rule_class"] for r in classification_rows) semantic_counts = Counter(r["final_semantic_class"] for r in classification_rows) segment_counts = Counter(r["segment_class"] for r in remainder_evidence) edge_ctx_counts = Counter( ( "CELL_EDGE" if r["invalid_point_on_cell_edge"] else "INTERNAL", r["final_semantic_class"], ) for r in classification_rows ) object_counts = Counter( (r["object_class"], r["final_semantic_class"]) for r in classification_rows ) report["summary"] = { "invalid_total": len(classification_rows), "remainder_class_counts": dict(remainder_counts), "fill_rule_class_counts": dict(fill_counts), "final_semantic_class_counts": dict(semantic_counts), "remainder_segment_class_counts": dict(segment_counts), "linear_remainder_record_count": sum( r["remainder_class"] != "NO_LINEAR_REMAINDER" for r in classification_rows ), "retrace_confirmed_record_count": sum( r["remainder_class"] == "RETRACED_CONNECTOR_CONFIRMED" for r in classification_rows ), "fill_rule_resolved_count": sum( r["fill_rule_class"] != "MATCHES_NEITHER" for r in classification_rows ), "fill_rule_odd_even_compatible_count": sum( r["fill_rule_class"] in { "MATCHES_BOTH_ODD_EVEN_AND_NONZERO", "MATCHES_ODD_EVEN_ONLY", } for r in classification_rows ), "fill_rule_nonzero_compatible_count": sum( r["fill_rule_class"] in { "MATCHES_BOTH_ODD_EVEN_AND_NONZERO", "MATCHES_NONZERO_ONLY", } for r in classification_rows ), } print("[4/5] 写文件 ...") write_csv( outdir / "invalid_polygon_fillrule_classification.csv", classification_rows ) write_csv( outdir / "linear_remainder_segment_evidence.csv", remainder_evidence ) write_csv( outdir / "polygon_face_fillrule_evidence.csv", face_evidence ) write_csv( outdir / "retraced_connector_summary.csv", [ {"remainder_class": k, "count": v} for k, v in remainder_counts.most_common() ], ["remainder_class", "count"] ) write_csv( outdir / "fillrule_by_object_class.csv", [ { "object_class": obj, "final_semantic_class": cls, "count": n, } for (obj, cls), n in object_counts.most_common() ] ) write_csv( outdir / "fillrule_by_boundary_context.csv", [ { "boundary_context": ctx, "final_semantic_class": cls, "count": n, } for (ctx, cls), n in edge_ctx_counts.most_common() ] ) make_geojson( records, outdir / "fillrule_diagnostic.geojson" ) (outdir / "fillrule_audit_report_v8.json").write_text( json.dumps(report, ensure_ascii=False, indent=2, default=str), encoding="utf-8" ) lines = [] lines.append("CM93 -> MySQL Invalid Polygon 线性残片与填充规则诊断 v8") lines.append("=" * 92) 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 为局部数据集；跨未导入海域不连续不作为 FAIL。") lines.append(" 本轮只诊断 retraced connector 与 polygon fill rule，不写数据库。") lines.append("") lines.append("REMAINDER RECORD CLASS") for k, v in remainder_counts.most_common(): lines.append(f" {k:42s} {v}") lines.append("") lines.append("REMAINDER SEGMENT EVIDENCE") for k, v in segment_counts.most_common(): lines.append(f" {k:42s} {v}") lines.append("") lines.append("FILL RULE") for k, v in fill_counts.most_common(): lines.append(f" {k:42s} {v}") lines.append("") lines.append("FINAL SEMANTIC CLASS") for k, v in semantic_counts.most_common(): lines.append(f" {k:42s} {v}") lines.append("") lines.append("KEY COUNTS") lines.append( f" linear remainder records: " f"{report['summary']['linear_remainder_record_count']}" ) lines.append( f" retraced connector confirmed: " f"{report['summary']['retrace_confirmed_record_count']}" ) lines.append( f" fill rule resolved: " f"{report['summary']['fill_rule_resolved_count']}" ) lines.append( f" odd-even compatible records: " f"{report['summary']['fill_rule_odd_even_compatible_count']}" ) lines.append( f" non-zero compatible records: " f"{report['summary']['fill_rule_nonzero_compatible_count']}" ) lines.append("") lines.append("INTERPRETATION") lines.append( " RETRACED_CONNECTOR_CONFIRMED：make_valid 的线性残片在原始 ring " "中由相反方向的重走 segment 解释，强烈支持“桥接/connector”而非随机坏线。" ) lines.append( " MATCHES_ODD_EVEN_ONLY / MATCHES_NONZERO_ONLY：" "说明 polygonize 的全部 planar faces 不能全部填充，必须应用对应 fill rule。" ) lines.append( " MATCHES_NEITHER：不能仅靠常见 odd-even / non-zero 规则解释，需要继续调查。" ) lines.append("") lines.append("NO WRITE / NO REPAIR / NO PBF") (outdir / "fillrule_audit_summary_v8.txt").write_text( "\n".join(lines), encoding="utf-8" ) conn.rollback() print("[5/5] 压缩 ...") zpath = zip_dir(outdir) 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()