#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ CM93 -> MySQL Invalid Polygon OGC Validity 表达诊断 Version: v7.0 背景 ==== 前序审计已经确认： - Dataset ID=1 是日本沿海局部 CM93 子集，不是完整 CM93； - raw_geotype=16 / META 的 4406 条不是 geometry corruption； - Point / Line / Sounding 已通过； - Area 共有 52 个 Shapely invalid； - Coverage 共有 3 个 Shapely invalid； - V6 从 raw_feature_blob + VECTOR 原始 edge chain 按 OpenCPN 规则重建： * 52/52 raw chain 仍然 invalid； * 52/52 与 DB Polygon 几乎完全重合； * 0 edge index 越界； * 0 edge junction 断裂； * 0 Feature 重复引用同一 edge； * importer edge 顺序/方向/坐标转换造成 geometry 损坏的假设已基本排除。 V7 目的 ======= 只做“OGC 有效面表达诊断”，不做任何数据库修复。 对 52 invalid Area + 3 invalid Coverage： 1. 调用 Shapely make_valid()，但只在内存中。 2. 分解 make_valid 输出： - Polygon / MultiPolygon - LineString / MultiLineString - Point / MultiPoint 3. 对原始 invalid geometry.boundary 做： unary_union(boundary) -> polygonize_full() 4. 统计： - 可 polygonize 的 face 数； - cuts / dangles / invalid rings； - make_valid polygon parts 数； - make_valid linear remainder 数和总长度； - make_valid point remainder 数； - 原始 invalid polygon area； - make_valid polygonal area； - polygonize faces union area； - make_valid polygonal 与 polygonize union 的 symmetric difference； - 原始 boundary 与 make_valid polygon boundary 的 Hausdorff distance； - boundary symmetric-difference length； 5. 精确标记 invalid point 是否在 Cell bbox 边界。 6. 给出“诊断分类”，注意这不是自动修复建议： - POLYGON_ONLY_CONSISTENT make_valid 只产生 polygonal geometry，且与 polygonize faces 基本一致。 - POLYGON_ONLY_FACE_SELECTION_DIFFERS make_valid 只产生面，但与 polygonize 全部 faces 的结果明显不同， 说明不能简单假设“所有 face 都应该填充”。 - POLYGON_PLUS_LINEAR_REMAINDER make_valid 同时产生面 + 线，说明原始 linework 中有不能直接归入有效面内部的线段。 - POLYGON_PLUS_POINT_REMAINDER - POLYGON_PLUS_MIXED_REMAINDER - NO_POLYGONAL_RESULT - MAKE_VALID_FAILED 7. 单独按： - Area / Coverage - Cell 边界 / Cell 内部 - Object Class - validity reason 统计。 8. 生成 GeoJSON： - original invalid - invalidity point - make_valid polygonal result - make_valid linear remainder - polygonize faces 9. 输出是否存在“可以统一无损转成 Polygon/MultiPolygon”的证据。 注意：即使所有记录都能变成有效面，也不能在本任务写回 DB。 非常重要 ======== 当前 Dataset 只覆盖日本沿海一部分，因此： - 数据集外海域不连续不算错误； - Cell 边界上的 self-touch / shared topology 要保留边界上下文； - V7 不把 Cell 边界异常自动判成“坏数据”； - V7 只回答“OGC Polygon/MultiPolygon 能否干净表达”。 禁止 ==== - 不 UPDATE / INSERT / DELETE - 不 ST_MakeValid 写回 - 不 buffer(0) - 不替换 chart_area_geometry - 不替换 chart_coverage - 不修改 geometry_issue - 不重新导入 - 不开始 PBF 依赖 ==== pip install pymysql shapely 执行 ==== python cm93_mysql_ogc_validity_audit_v7.py ^ --host 127.0.0.1 ^ --port 3306 ^ --user root ^ --database cm93_chart ^ --dataset-id 1 输出 ==== - ogc_validity_summary_v7.txt - ogc_validity_report_v7.json - invalid_polygon_ogc_diagnostic.csv - make_valid_component_detail.csv - polygonize_face_detail.csv - ogc_diagnostic_by_class.csv - ogc_diagnostic_by_object_class.csv - ogc_diagnostic_by_cell.csv - ogc_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, 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 fetch_all(cur, sql: str, params: Sequence[Any] = ()) -> List[Dict[str, Any]]: cur.execute(sql, params) return list(cur.fetchall()) 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, tuple, list)): 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 ds = { "WEST": abs(x - minx), "EAST": abs(x - maxx), "SOUTH": abs(y - miny), "NORTH": abs(y - maxy), } edge, dist = min(ds.items(), key=lambda kv: kv[1]) return int(dist <= tol), edge if dist <= tol else "", dist def validity_reason_class(reason: str) -> str: s = (reason or "").lower() if "ring self-intersection" in s: return "RING_SELF_INTERSECTION" if "self-intersection" in s: return "SELF_INTERSECTION" if "interior is disconnected" in s: return "INTERIOR_DISCONNECTED" if "hole lies outside shell" in s: return "HOLE_OUTSIDE_SHELL" if "nested shells" in s: return "NESTED_SHELLS" if "nested holes" in s: return "NESTED_HOLES" if "too few points" in s: return "TOO_FEW_POINTS" return "OTHER_INVALID" 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 safe_area(g) -> Optional[float]: if g is None: return None try: return float(g.area) except Exception: return None def safe_length(g) -> Optional[float]: if g is None: return None try: return float(g.length) 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_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_symdiff_length(a, b) -> Optional[float]: if a is None or b is None: return None try: return float(a.symmetric_difference(b).length) except Exception: return None def count_collection_parts(g) -> int: if g is None: return 0 try: if hasattr(g, "geoms"): return len(g.geoms) return 1 except Exception: return 0 def diagnostic_class( *, polygon_count: int, line_count: int, point_count: int, other_count: int, polygonize_agrees: Optional[bool], ) -> str: if polygon_count == 0: return "NO_POLYGONAL_RESULT" if line_count == 0 and point_count == 0 and other_count == 0: if polygonize_agrees is True: return "POLYGON_ONLY_CONSISTENT" if polygonize_agrees is False: return "POLYGON_ONLY_FACE_SELECTION_DIFFERS" return "POLYGON_ONLY_UNRESOLVED_COMPARISON" if line_count > 0 and point_count == 0 and other_count == 0: return "POLYGON_PLUS_LINEAR_REMAINDER" if line_count == 0 and point_count > 0 and other_count == 0: return "POLYGON_PLUS_POINT_REMAINDER" return "POLYGON_PLUS_MIXED_REMAINDER" 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, f.geometry_issue, f.geometry_repaired, 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, v.geometry_issue, v.geometry_repaired, 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"), "diagnostic_class": r.get("_diagnostic_class"), "validity_reason": r.get("_validity_reason"), "cell_edge": r.get("_invalid_point_on_cell_edge"), } add(r.get("_geom"), { **common, "layer": "original_invalid", }) p = r.get("_invalid_point") if p is not None: add(Point(p), { **common, "layer": "invalidity_point", }) add(r.get("_mv_polygon_union"), { **common, "layer": "make_valid_polygonal", }) add(r.get("_mv_line_union"), { **common, "layer": "make_valid_linear_remainder", }) add(r.get("_polygonize_union"), { **common, "layer": "polygonize_faces_union", }) 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 OGC validity audit v7" ) 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( "--cell-edge-tolerance", type=float, default=5e-6, help="invalid point 距 Cell bbox 边界判定容差（degree）" ) p.add_argument( "--face-symdiff-abs-tolerance", type=float, default=1e-12, help="make_valid polygon 与 polygonize union 面积差绝对容差（degree^2）" ) p.add_argument( "--face-symdiff-relative-tolerance", type=float, default=1e-6, help="make_valid polygon 与 polygonize union 面积差相对容差" ) 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_ogc_validity_audit_v7_dataset_{args.dataset_id}_{stamp}" ) outdir.mkdir(parents=True, exist_ok=True) report = { "version": "v7.0", "generated_at": datetime.now().isoformat(timespec="seconds"), "database": args.database, "dataset_id": args.dataset_id, "read_only": True, "pbf": "OUT_OF_SCOPE", "scope_note": ( "CM93JAPAN 是日本沿海局部数据集。" "V7 只诊断 invalid Polygon 如何映射到 OGC valid polygonal geometry，" "不做任何数据库写回。" ), "warnings": [], "shapely_version": getattr(shapely, "__version__", "unknown"), } 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" total={len(records)} " f"Area={sum(r['kind']=='AREA' for r in records)} " f"Coverage={sum(r['kind']=='COVERAGE' for r in records)}" ) print("[2/5] make_valid + polygonize_full 诊断 ...") diag_rows = [] component_rows = [] face_rows = [] for i, r in enumerate(records, start=1): geom = r["_geom"] reason = explain_validity(geom) reason_class = validity_reason_class(reason) 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 ) base = { "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"), "validity_reason": reason, "validity_reason_class": reason_class, "invalid_x": invalid_point[0] if invalid_point else None, "invalid_y": invalid_point[1] if invalid_point else None, "invalid_point_on_cell_edge": on_edge, "cell_edge": edge_name, "cell_edge_distance_deg": edge_dist, } try: mv = make_valid_geom(geom) except Exception as e: diag_rows.append({ **base, "make_valid_type": None, "diagnostic_class": "MAKE_VALID_FAILED", "detail": f"{type(e).__name__}: {e}", }) continue comps = flatten_components(mv) mv_poly_union = geom_union(comps["polygons"]) mv_line_union = geom_union(comps["lines"]) mv_point_union = geom_union(comps["points"]) for idx, p in enumerate(comps["polygons"]): component_rows.append({ **base, "component_dimension": "POLYGON", "component_index": idx, "geom_type": p.geom_type, "area_deg2": safe_area(p), "length_deg": safe_length(p.boundary), }) for idx, l in enumerate(comps["lines"]): component_rows.append({ **base, "component_dimension": "LINE", "component_index": idx, "geom_type": l.geom_type, "area_deg2": 0.0, "length_deg": safe_length(l), }) for idx, p in enumerate(comps["points"]): component_rows.append({ **base, "component_dimension": "POINT", "component_index": idx, "geom_type": p.geom_type, "area_deg2": 0.0, "length_deg": 0.0, }) for idx, o in enumerate(comps["others"]): component_rows.append({ **base, "component_dimension": "OTHER", "component_index": idx, "geom_type": getattr(o, "geom_type", type(o).__name__), "area_deg2": safe_area(o), "length_deg": safe_length(o), }) # polygonize_full on noded linework polygonize_error = "" polygonize_union = None polygons_geom = cuts = dangles = invalid_rings = None poly_faces: List[Any] = [] try: noded_boundary = unary_union(geom.boundary) polygons_geom, cuts, dangles, invalid_rings = polygonize_full( noded_boundary ) pcomps = flatten_components(polygons_geom) poly_faces = pcomps["polygons"] polygonize_union = geom_union(poly_faces) except Exception as e: polygonize_error = f"{type(e).__name__}: {e}" for idx, face in enumerate(poly_faces): face_rows.append({ **base, "face_index": idx, "area_deg2": safe_area(face), "boundary_length_deg": safe_length(face.boundary), }) orig_area = safe_area(geom) mv_poly_area = safe_area(mv_poly_union) polyz_area = safe_area(polygonize_union) mv_vs_polyz_symarea = safe_symdiff_area( mv_poly_union, polygonize_union ) if ( mv_vs_polyz_symarea is not None and mv_poly_area is not None and polyz_area is not None ): scale = max(abs(mv_poly_area), abs(polyz_area), 1e-30) allowed = max( args.face_symdiff_abs_tolerance, args.face_symdiff_relative_tolerance * scale ) polygonize_agrees = mv_vs_polyz_symarea <= allowed else: allowed = None polygonize_agrees = None cls = diagnostic_class( polygon_count=len(comps["polygons"]), line_count=len(comps["lines"]), point_count=len(comps["points"]), other_count=len(comps["others"]), polygonize_agrees=polygonize_agrees, ) original_boundary = None mv_boundary = None try: original_boundary = unary_union(geom.boundary) except Exception: pass try: if mv_poly_union is not None: mv_boundary = unary_union(mv_poly_union.boundary) except Exception: pass boundary_haus = safe_hausdorff( original_boundary, mv_boundary ) boundary_sym_len = safe_symdiff_length( original_boundary, mv_boundary ) line_remainder_len = sum( safe_length(x) or 0.0 for x in comps["lines"] ) point_remainder_count = len(comps["points"]) diag = { **base, "original_geom_type": geom.geom_type, "original_area_deg2": orig_area, "make_valid_type": mv.geom_type, "make_valid_is_valid": int(bool(mv.is_valid)), "make_valid_polygon_count": len(comps["polygons"]), "make_valid_line_count": len(comps["lines"]), "make_valid_point_count": len(comps["points"]), "make_valid_other_count": len(comps["others"]), "make_valid_polygon_area_deg2": mv_poly_area, "make_valid_linear_remainder_length_deg": line_remainder_len, "make_valid_point_remainder_count": point_remainder_count, "polygonize_face_count": len(poly_faces), "polygonize_union_area_deg2": polyz_area, "polygonize_cuts_length_deg": safe_length(cuts), "polygonize_dangles_length_deg": safe_length(dangles), "polygonize_invalid_rings_length_deg": safe_length(invalid_rings), "polygonize_error": polygonize_error, "make_valid_vs_polygonize_symdiff_area_deg2": mv_vs_polyz_symarea, "face_agreement_allowed_area_deg2": allowed, "make_valid_polygonize_agree": ( int(polygonize_agrees) if polygonize_agrees is not None else None ), "original_vs_make_valid_boundary_hausdorff_deg": boundary_haus, "original_vs_make_valid_boundary_symdiff_length_deg": boundary_sym_len, "diagnostic_class": cls, } diag_rows.append(diag) r["_validity_reason"] = reason r["_invalid_point"] = invalid_point r["_invalid_point_on_cell_edge"] = on_edge r["_mv_polygon_union"] = mv_poly_union r["_mv_line_union"] = mv_line_union r["_mv_point_union"] = mv_point_union r["_polygonize_union"] = polygonize_union r["_diagnostic_class"] = cls rid = ( r.get("feature_id") if r["kind"] == "AREA" else r.get("coverage_geom_id") ) print( f" [{i:02d}/{len(records)}] " f"{r['kind']} id={rid} {r.get('object_class')} " f"{reason_class} -> {cls} " f"poly={len(comps['polygons'])} " f"line={len(comps['lines'])}" ) print("[3/5] 汇总分类 ...") class_counts = Counter( r.get("diagnostic_class", "UNKNOWN") for r in diag_rows ) kind_class_counts = Counter( (r.get("kind"), r.get("diagnostic_class")) for r in diag_rows ) edge_class_counts = Counter( ( "CELL_EDGE" if r.get("invalid_point_on_cell_edge") else "INTERNAL", r.get("diagnostic_class"), ) for r in diag_rows ) object_counts = Counter( (r.get("object_class"), r.get("diagnostic_class")) for r in diag_rows ) cell_counts = Counter( (r.get("cell_name"), r.get("diagnostic_class")) for r in diag_rows ) report["summary"] = { "invalid_total": len(diag_rows), "invalid_area_count": sum( r.get("kind") == "AREA" for r in diag_rows ), "invalid_coverage_count": sum( r.get("kind") == "COVERAGE" for r in diag_rows ), "diagnostic_class_counts": dict(class_counts), "cell_edge_count": sum( int(r.get("invalid_point_on_cell_edge") or 0) for r in diag_rows ), "internal_count": sum( 1 for r in diag_rows if not r.get("invalid_point_on_cell_edge") ), "make_valid_with_linear_remainder": sum( int(r.get("make_valid_line_count") or 0) > 0 for r in diag_rows ), "make_valid_polygon_only": sum( int(r.get("make_valid_polygon_count") or 0) > 0 and int(r.get("make_valid_line_count") or 0) == 0 and int(r.get("make_valid_point_count") or 0) == 0 and int(r.get("make_valid_other_count") or 0) == 0 for r in diag_rows ), "make_valid_vs_polygonize_agree": sum( int(r.get("make_valid_polygonize_agree") or 0) for r in diag_rows ), "records_with_polygonize_dangles": sum( (r.get("polygonize_dangles_length_deg") or 0.0) > 0 for r in diag_rows ), "records_with_polygonize_cuts": sum( (r.get("polygonize_cuts_length_deg") or 0.0) > 0 for r in diag_rows ), } by_class_rows = [] for (kind, cls), n in kind_class_counts.most_common(): by_class_rows.append({ "group": "KIND", "key": kind, "diagnostic_class": cls, "count": n, }) for (ctx, cls), n in edge_class_counts.most_common(): by_class_rows.append({ "group": "BOUNDARY_CONTEXT", "key": ctx, "diagnostic_class": cls, "count": n, }) obj_rows = [ { "object_class": obj, "diagnostic_class": cls, "count": n, } for (obj, cls), n in object_counts.most_common() ] cell_rows = [ { "cell_name": cell, "diagnostic_class": cls, "count": n, } for (cell, cls), n in cell_counts.most_common() ] print("[4/5] 写输出文件 ...") write_csv( outdir / "invalid_polygon_ogc_diagnostic.csv", diag_rows ) write_csv( outdir / "make_valid_component_detail.csv", component_rows ) write_csv( outdir / "polygonize_face_detail.csv", face_rows ) write_csv( outdir / "ogc_diagnostic_by_class.csv", by_class_rows ) write_csv( outdir / "ogc_diagnostic_by_object_class.csv", obj_rows ) write_csv( outdir / "ogc_diagnostic_by_cell.csv", cell_rows ) make_geojson( [r for r in records if "_diagnostic_class" in r], outdir / "ogc_diagnostic.geojson" ) (outdir / "ogc_validity_report_v7.json").write_text( json.dumps(report, ensure_ascii=False, indent=2, default=str), encoding="utf-8" ) lines = [] lines.append("CM93 -> MySQL Invalid Polygon OGC Validity 表达诊断 v7") lines.append("=" * 90) 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 是局部数据集；Cell/Dataset 外边界不连续不作为错误。" ) lines.append( " 本轮仅在内存中运行 make_valid / polygonize_full 作为诊断，" "绝不写回数据库。" ) lines.append("") lines.append("BASELINE") lines.append( f" invalid total: {report['summary']['invalid_total']}" ) lines.append( f" Area: {report['summary']['invalid_area_count']}" ) lines.append( f" Coverage: {report['summary']['invalid_coverage_count']}" ) lines.append( f" invalid point Cell edge: {report['summary']['cell_edge_count']}" ) lines.append( f" invalid point internal: {report['summary']['internal_count']}" ) lines.append("") lines.append("DIAGNOSTIC CLASS") for k, v in class_counts.most_common(): lines.append(f" {k:46s} {v}") lines.append("") lines.append("MAKE_VALID COMPOSITION") lines.append( f" polygon-only records: " f"{report['summary']['make_valid_polygon_only']}" ) lines.append( f" records with linear remainder: " f"{report['summary']['make_valid_with_linear_remainder']}" ) lines.append("") lines.append("POLYGONIZE") lines.append( f" make_valid polygon agrees with polygonize faces: " f"{report['summary']['make_valid_vs_polygonize_agree']}" ) lines.append( f" records with polygonize dangles: " f"{report['summary']['records_with_polygonize_dangles']}" ) lines.append( f" records with polygonize cuts: " f"{report['summary']['records_with_polygonize_cuts']}" ) lines.append("") lines.append("HOW TO READ") lines.append( " POLYGON_ONLY_CONSISTENT：存在比较干净的 OGC Polygon/MultiPolygon 表达候选，" "但本轮仍不允许写回。" ) lines.append( " POLYGON_PLUS_LINEAR_REMAINDER：make_valid 为了得到有效面，" "还必须留下线性几何；这说明不能粗暴把结果只取 Polygon 后当成无损修复。" ) lines.append( " POLYGON_ONLY_FACE_SELECTION_DIFFERS：make_valid 与把所有 noded faces " "全部 polygonize 的结果不同，说明 fill/face 选择存在语义问题。" ) lines.append( " 当前目标只是确定表达性质，不决定 importer 是否要规范化。" ) lines.append("") lines.append("NO WRITE / NO REPAIR / NO PBF") (outdir / "ogc_validity_summary_v7.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()