#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ CM93 -> MySQL invalid Area 原始 Edge Chain 重建与 Importer 对账 Version: v6.0 背景 ==== Dataset ID=1 为日本沿海局部 CM93 数据，并非完整 CM93。 V3/V4/V5 已确认： - 4406 条 raw_geotype=16 / META 是正常 parent-no-geometry，不是真 geometry 损坏； - Point / Line / Sounding 均通过； - Area 81080 中仅 52 个 Shapely invalid； - Coverage 1156 中仅 3 个 invalid； - 52 个 invalid Area 中 50 个的 invalid point 还被至少一个“有效 Area”边界共享； - chart_feature.raw_feature_blob 已逐 Feature 保存原始（已解码）Feature record； - chart_raw_record 中保存每个 Cell 的 HEADER / VECTOR 等整表 blob。 本脚本目的 ========== 只针对 invalid Area： 1. 从 chart_feature.raw_feature_blob 解码 CM93 Area 的 edge references： - object_type - geotype - obj_desc_bytes - n_elements - 每个 geometry word - edge_index = word & 0x1fff - segment_usage = word >> 13 - usage bit 4 表示该 edge 反向使用 2. 从 chart_raw_record(record_kind='VECTOR') 解析该 Cell 的全部原始 edge： - ushort npoints - npoints * (ushort x, ushort y) 3. 按 OpenCPN CM93 BuildGeom 的规则重建 ring： - 根据 segment_usage 选择正向/反向； - 每段拼接时跳过前一段共享的最后点； - 当前 edge 终点回到 ring 起点时结束一个 ring。 4. 检查： - Feature edge index 是否越界； - 相邻 edge 是否连续； - ring 是否闭合； - 同一 edge 是否重复使用； - ring 是否有非连续重复顶点； - raw ring 在原始 0..65535 坐标空间是否自交/无效； - raw invalidity 与数据库 WGS84 Polygon invalidity 是否拓扑一致。 5. 使用 HEADER 中 easting/northing transform + chart_coverage.wgsox/wgsoy， 将 raw CM93 点转换到 WGS84，重建参考 Polygon。 6. 将“参考重建 Polygon”与数据库 chart_area_geometry.geom 对账： - geometry type - validity - validity reason - ring 数 - Hausdorff distance - symmetric difference area（仅辅助，不作为面积物理量） 7. 对 52 条给出表现分类： A. RAW_EDGE_CHAIN_ITSELF_INVALID 原始 edge chain 按 OpenCPN 规则组出来已经 invalid； 更接近 CM93 拓扑表达 / 数据本身 / OGC Polygon 表达差异。 B. RAW_CHAIN_VALID_DB_INVALID 原始 edge chain 重建有效，但 DB Polygon 无效； 强烈指向 importer geometry assembly / conversion bug。 C. EDGE_CHAIN_DISCONTINUITY Feature 引用的相邻 edge 不连续或 ring 不闭合； 需查 edge ordering / source data。 D. RECONSTRUCTION_MISMATCH 参考重建与 DB geometry 明显不同； 指向 importer transform / ring assembly / coverage offset 差异。 E. RECONSTRUCTION_MATCHES_DB 参考重建和 DB topology/shape 基本一致。 8. 输出最值得进一步人工看的孤立 Feature（V5 中无共享邻居的 2 条）详细 edge chain。 9. 全程只读；不 make_valid、不 buffer(0)、不修改 DB、不做 PBF。 参考 ==== 算法按 OpenCPN 当前 cm93.cpp 的公开实现复核： - read_vector_record_table() - read_feature_record_table() - cm93chart::BuildGeom() 注意： 本脚本是“对账/取证工具”，不是修复器。 即使 raw chain 本身 invalid，也不能直接断言 CM93 原始数据“坏了”。 CM93 的共享拓扑表达转换为 OGC Polygon 后，可能出现 self-touch 等 validity 差异。 依赖 ==== pip install pymysql shapely 执行 ==== python cm93_mysql_raw_edge_audit_v6.py ^ --host 127.0.0.1 ^ --port 3306 ^ --user root ^ --database cm93_chart ^ --dataset-id 1 输出 ==== 自动生成目录和 ZIP： - raw_edge_audit_summary_v6.txt - raw_edge_audit_report_v6.json - invalid_area_raw_edge_classification.csv - invalid_area_edge_references.csv - invalid_area_ring_summary.csv - invalid_area_edge_junction_issues.csv - raw_vs_db_geometry_compare.csv - isolated_invalid_area_edge_detail.csv - affected_cell_vector_table_audit.csv - invalid_area_raw_reconstruction.geojson 安全 ==== 只执行 SELECT / read-only transaction。 """ from __future__ import annotations import argparse import csv import getpass import json import math import os import struct 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: 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) from shapely.geometry import ( Point, LineString, LinearRing, Polygon, MultiPolygon, mapping ) from shapely.validation import explain_validity except ImportError: print("缺少 shapely：pip install pymysql shapely") sys.exit(2) CM93_SEMIMAJOR_AXIS_METERS = 6378388.0 DEGREE = math.pi / 180.0 def qi(name: str) -> str: return "`" + name.replace("`", "``") + "`" 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) elif isinstance(v, (bytes, bytearray, memoryview)): v = bytes(v).hex() 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 函数") # --------------------------------------------------------------------------- # CM93 raw structures # --------------------------------------------------------------------------- def u16le(buf: bytes, offset: int) -> int: if offset + 2 > len(buf): raise ValueError(f"ushort 越界 offset={offset}, len={len(buf)}") return struct.unpack_from(" Dict[str, Any]: """ OpenCPN read_feature_record_table() 对 Area 的前部格式： uint8 object_type uint8 geom_prim uint16 obj_desc_bytes uint16 n_elements uint16 edge_ref[n_elements] ... related/attrs follow according to geotype flags edge_ref: edge_index = edge_ref & 0x1fff segment_usage = edge_ref >> 13 reversed when segment_usage & 4 """ if raw is None or len(raw) < 6: raise ValueError("raw_feature_blob 太短") object_type = raw[0] geotype = raw[1] obj_desc_bytes = u16le(raw, 2) primitive = geotype & 0x0F if primitive != 4: raise ValueError( f"不是 Area primitive: geotype={geotype}, low_nibble={primitive}" ) n_elements = u16le(raw, 4) need = 6 + n_elements * 2 if len(raw) < need: raise ValueError( f"raw_feature_blob edge refs 不完整: " f"n_elements={n_elements}, need={need}, len={len(raw)}" ) refs = [] pos = 6 for i in range(n_elements): word = u16le(raw, pos) pos += 2 usage = word >> 13 edge_index = word & 0x1FFF refs.append({ "sequence": i, "raw_word": word, "raw_word_hex": f"0x{word:04x}", "edge_index": edge_index, "segment_usage": usage, "reverse": int(bool(usage & 4)), "usage_bit0": int(bool(usage & 1)), "usage_bit1": int(bool(usage & 2)), "usage_bit2_reverse": int(bool(usage & 4)), }) return { "object_type": object_type, "geotype": geotype, "obj_desc_bytes": obj_desc_bytes, "primitive": primitive, "has_children_related": int(bool(geotype & 0x10)), "has_related_count_word": int(bool(geotype & 0x20)), "has_bit_0x40": int(bool(geotype & 0x40)), "has_attributes": int(bool(geotype & 0x80)), "n_elements": n_elements, "geometry_refs_end_offset": pos, "raw_length": len(raw), "refs": refs, } def parse_vector_table(raw: bytes) -> List[List[Tuple[int, int]]]: """ OpenCPN read_vector_record_table(): for each edge: uint16 npoints repeat npoints: uint16 x uint16 y """ if raw is None: raise ValueError("VECTOR raw_blob is NULL") edges: List[List[Tuple[int, int]]] = [] pos = 0 while pos < len(raw): if pos + 2 > len(raw): raise ValueError( f"VECTOR 尾部不足 2 bytes: pos={pos}, len={len(raw)}" ) npoints = u16le(raw, pos) pos += 2 need = npoints * 4 if pos + need > len(raw): raise ValueError( f"VECTOR edge#{len(edges)} 数据越界: " f"npoints={npoints}, pos={pos}, len={len(raw)}" ) pts = [] for _ in range(npoints): x = u16le(raw, pos) y = u16le(raw, pos + 2) pos += 4 pts.append((x, y)) edges.append(pts) if pos != len(raw): raise ValueError(f"VECTOR 未精确消费: pos={pos}, len={len(raw)}") return edges def parse_header_metadata(raw_metadata_json: Any) -> Dict[str, Any]: if raw_metadata_json is None: raise ValueError("HEADER raw_metadata_json is NULL") if isinstance(raw_metadata_json, (bytes, bytearray, memoryview)): s = bytes(raw_metadata_json).decode("utf-8", errors="replace") obj = json.loads(s) elif isinstance(raw_metadata_json, str): obj = json.loads(raw_metadata_json) elif isinstance(raw_metadata_json, dict): obj = raw_metadata_json else: obj = json.loads(str(raw_metadata_json)) h = obj.get("header", obj) required = [ "easting_min", "easting_max", "northing_min", "northing_max", "usn_vector_records" ] missing = [k for k in required if k not in h] if missing: raise ValueError(f"HEADER metadata 缺字段: {missing}") return h def transform_cm93_point_to_wgs84( p: Tuple[int, int], header: Dict[str, Any], wgsox: float = 0.0, wgsoy: float = 0.0, ) -> Tuple[float, float]: """ 对齐 OpenCPN： transform_x_rate = delta_x / 65535 transform_y_rate = delta_y / 65535 valx = x * rate + origin valy = y * rate + origin valx -= wgsox valy -= wgsoy lat = inverse mercator lon = valx / (DEGREE * CM93_semimajor_axis_meters) 返回 (lon, lat) """ x, y = p e_min = float(header["easting_min"]) e_max = float(header["easting_max"]) n_min = float(header["northing_min"]) n_max = float(header["northing_max"]) dx = e_max - e_min if dx < 0: dx += CM93_SEMIMAJOR_AXIS_METERS * 2.0 * math.pi x_rate = dx / 65535.0 y_rate = (n_max - n_min) / 65535.0 valx = x * x_rate + e_min valy = y * y_rate + n_min valx -= float(wgsox or 0.0) valy -= float(wgsoy or 0.0) lat = ( 2.0 * math.atan(math.exp(valy / CM93_SEMIMAJOR_AXIS_METERS)) - math.pi / 2.0 ) / DEGREE lon = valx / (DEGREE * CM93_SEMIMAJOR_AXIS_METERS) return lon, lat def orient_edge(points: List[Tuple[int, int]], reverse: bool) -> List[Tuple[int, int]]: return list(reversed(points)) if reverse else list(points) def build_rings_like_opencpn( refs: List[Dict[str, Any]], edges: List[List[Tuple[int, int]]], ) -> Tuple[List[List[Tuple[int, int]]], List[Dict[str, Any]], List[Dict[str, Any]]]: """ 按 OpenCPN BuildGeom Area 的 ring 判定方式重建。 OpenCPN 对每个 segment 把全部点写入，然后 ip-- 跳掉该段最后点； 当 cur_end_point == start_point 时，把 closure point 留下并结束 ring。 这里用更直观的列表实现，但保留相同语义。 """ rings: List[List[Tuple[int, int]]] = [] edge_rows: List[Dict[str, Any]] = [] junction_issues: List[Dict[str, Any]] = [] current: List[Tuple[int, int]] = [] ring_start: Optional[Tuple[int, int]] = None previous_end: Optional[Tuple[int, int]] = None ring_index = 0 for ref in refs: ei = int(ref["edge_index"]) if ei < 0 or ei >= len(edges): edge_rows.append({ **ref, "edge_exists": 0, "edge_point_count": None, "ring_index": ring_index, }) junction_issues.append({ "ring_index": ring_index, "sequence": ref["sequence"], "issue_category": "EDGE_INDEX_OUT_OF_RANGE", "detail": f"edge_index={ei}, vector_edge_count={len(edges)}", }) continue raw_pts = edges[ei] pts = orient_edge(raw_pts, bool(ref["reverse"])) if not pts: edge_rows.append({ **ref, "edge_exists": 1, "edge_point_count": 0, "ring_index": ring_index, }) junction_issues.append({ "ring_index": ring_index, "sequence": ref["sequence"], "issue_category": "EMPTY_EDGE", "detail": f"edge_index={ei}", }) continue start = pts[0] end = pts[-1] if ring_start is None: ring_start = start previous_end = None continuous_from_previous = ( previous_end is None or previous_end == start ) if previous_end is not None and previous_end != start: junction_issues.append({ "ring_index": ring_index, "sequence": ref["sequence"], "issue_category": "EDGE_JUNCTION_DISCONTINUITY", "detail": ( f"previous_end={previous_end}, current_start={start}, " f"edge_index={ei}" ), "previous_end_x": previous_end[0], "previous_end_y": previous_end[1], "current_start_x": start[0], "current_start_y": start[1], }) # OpenCPN 在 segment 之间跳过共享 end point。 # 若连续，当前 ring 已包含 previous_end；所以追加本 edge 除首点外的点。 # 第一段则全部加入。 if not current: current.extend(pts) else: if current[-1] == pts[0]: current.extend(pts[1:]) else: # 为保持 source 顺序证据，不自动补连；直接追加完整 edge。 current.extend(pts) closed_here = (end == ring_start) edge_rows.append({ **ref, "edge_exists": 1, "edge_point_count": len(raw_pts), "oriented_start_x": start[0], "oriented_start_y": start[1], "oriented_end_x": end[0], "oriented_end_y": end[1], "continuous_from_previous": int(continuous_from_previous), "ring_index": ring_index, "closes_ring": int(closed_here), }) previous_end = end if closed_here: if current and current[0] != current[-1]: current.append(current[0]) rings.append(current) current = [] ring_start = None previous_end = None ring_index += 1 if current: junction_issues.append({ "ring_index": ring_index, "sequence": None, "issue_category": "RING_NOT_CLOSED_AT_END_OF_REFS", "detail": ( f"unclosed ring vertices={len(current)}, " f"start={ring_start}, end={current[-1] if current else None}" ), }) rings.append(current) return rings, edge_rows, junction_issues # --------------------------------------------------------------------------- # Geometry construction and comparison # --------------------------------------------------------------------------- def repeated_nonconsecutive_vertices( ring: List[Tuple[int, int]] ) -> List[Tuple[int, int]]: if not ring: return [] core = ring[:-1] if len(ring) > 1 and ring[0] == ring[-1] else ring positions = defaultdict(list) for i, p in enumerate(core): positions[p].append(i) out = [] n = len(core) for p, pos in positions.items(): if len(pos) < 2: continue nonconsecutive = False for i in range(len(pos)): for j in range(i + 1, len(pos)): d = abs(pos[j] - pos[i]) if d > 1 and d < n - 1: nonconsecutive = True break if nonconsecutive: break if nonconsecutive: out.append(p) return out def signed_area(coords: List[Tuple[float, float]]) -> float: if len(coords) < 3: return 0.0 s = 0.0 for i in range(len(coords) - 1): x1, y1 = coords[i] x2, y2 = coords[i + 1] s += x1 * y2 - x2 * y1 return s / 2.0 def ring_to_polygon_candidate( rings: List[List[Tuple[float, float]]] ): """ 为“对账”构造 Polygon/MultiPolygon。 不尝试 make_valid。 CM93/OpenCPN BuildGeom 只记录 contour 顺序，不在这里做复杂 shell/hole 修复。 优先策略： - 1 ring -> Polygon(shell) - 多 ring -> 先把第一个 ring 作为 shell，其余作为 holes； 若失败，再以每个 ring 独立 Polygon 组成 MultiPolygon，仅用于比较。 返回 (geom, construction_mode) """ clean = [] for ring in rings: if len(ring) >= 4: rr = list(ring) if rr[0] != rr[-1]: rr.append(rr[0]) clean.append(rr) if not clean: return None, "NO_USABLE_RING" if len(clean) == 1: try: return Polygon(clean[0]), "SINGLE_RING_POLYGON" except Exception: return None, "SINGLE_RING_CONSTRUCTION_FAILED" try: return Polygon(clean[0], clean[1:]), "FIRST_SHELL_REST_HOLES" except Exception: pass try: polys = [Polygon(r) for r in clean] return MultiPolygon(polys), "MULTIPOLYGON_PER_RING_FALLBACK" except Exception: return None, "MULTIRING_CONSTRUCTION_FAILED" def safe_validity(g) -> Tuple[Optional[bool], str]: if g is None: return None, "NULL_GEOMETRY" try: valid = bool(g.is_valid) return valid, "Valid Geometry" if valid else explain_validity(g) except Exception as e: return None, f"{type(e).__name__}: {e}" def safe_hausdorff(a, b) -> Optional[float]: try: return float(a.hausdorff_distance(b)) except Exception: return None def safe_symdiff_area(a, b) -> Optional[float]: try: return float(a.symmetric_difference(b).area) except Exception: return None def topology_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" in s: return "NESTED_RING" if "valid geometry" in s: return "VALID" return "OTHER" # --------------------------------------------------------------------------- # DB loaders # --------------------------------------------------------------------------- def load_invalid_areas(cur, dataset_id: int, wkb_fn: str) -> List[Dict[str, Any]]: cur.execute(f""" SELECT f.id AS feature_id, f.dataset_id, f.cell_id, f.coverage_id, f.source_record_index, f.object_class, f.geometry_type, f.raw_geotype, f.normalized_srid, f.raw_feature_blob, f.geometry_issue, f.geometry_repaired, c.cell_name, {wkb_fn}(g.geom) AS geom_wkb, cv.wgsox AS coverage_wgsox, cv.wgsoy AS coverage_wgsoy, cv.source_record_index AS coverage_source_record_index 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 LEFT JOIN chart_coverage cv ON cv.id=f.coverage_id WHERE f.dataset_id=%s AND f.geometry_type='AREA' ORDER BY f.id """, (dataset_id,)) out = [] for r in cur: try: db_geom = load_wkb(to_bytes(r["geom_wkb"])) if not db_geom.is_valid: x = dict(r) x["_db_geom"] = db_geom x["_db_reason"] = explain_validity(db_geom) out.append(x) except Exception: continue return out def load_cell_raw_tables( cur, cell_id: int ) -> Dict[str, Dict[str, Any]]: rs = fetch_all(cur, """ SELECT id, record_kind, source_record_index, raw_blob, raw_metadata_json FROM chart_raw_record WHERE cell_id=%s ORDER BY id """, (cell_id,)) d = {} for r in rs: kind = str(r["record_kind"]).upper() # 当前 import 设计中每种 table 通常一个 chunk。 # 如未来出现 chunk_count>1，先保留全部并在调用端报错。 if kind in d: # 用 list 暂存重复。 existing = d[kind] if not isinstance(existing, dict) or "_chunks" not in existing: d[kind] = {"_chunks": [existing, r]} else: existing["_chunks"].append(r) else: d[kind] = r return d def normalize_chunk_record(rec: Dict[str, Any], kind: str) -> Dict[str, Any]: if "_chunks" not in rec: return rec chunks = rec["_chunks"] parsed = [] for c in chunks: meta = c.get("raw_metadata_json") try: if isinstance(meta, str): mo = json.loads(meta) elif isinstance(meta, dict): mo = meta else: mo = json.loads(bytes(meta).decode("utf-8")) except Exception: mo = {} parsed.append((int(mo.get("chunk_index", 0)), c)) parsed.sort(key=lambda t: t[0]) raw = b"".join(to_bytes(c.get("raw_blob")) or b"" for _, c in parsed) first = dict(parsed[0][1]) first["raw_blob"] = raw first["_joined_chunks"] = len(parsed) return first # --------------------------------------------------------------------------- # GeoJSON # --------------------------------------------------------------------------- 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: fid = r["feature_id"] add(r.get("_db_geom"), { "layer": "db_invalid_area", "feature_id": fid, "object_class": r.get("object_class"), "cell_name": r.get("cell_name"), "db_reason": r.get("_db_reason"), "final_class": r.get("_final_class"), }) add(r.get("_raw_wgs_geom"), { "layer": "raw_reconstructed_area", "feature_id": fid, "object_class": r.get("object_class"), "cell_name": r.get("cell_name"), "raw_reason": r.get("_raw_wgs_reason"), "construction_mode": r.get("_construction_mode"), }) path.write_text( json.dumps( {"type": "FeatureCollection", "features": features}, ensure_ascii=False, separators=(",", ":"), default=str, ), encoding="utf-8" ) def classify_record( *, junction_issues: List[Dict[str, Any]], raw_valid: Optional[bool], db_valid: bool, hausdorff: Optional[float], symdiff: Optional[float], compare_tolerance: float, ) -> Tuple[str, List[str]]: flags = [] if junction_issues: flags.extend(sorted(set(x["issue_category"] for x in junction_issues))) if any( x["issue_category"] in { "EDGE_INDEX_OUT_OF_RANGE", "EDGE_JUNCTION_DISCONTINUITY", "RING_NOT_CLOSED_AT_END_OF_REFS", "EMPTY_EDGE", } for x in junction_issues ): return "EDGE_CHAIN_DISCONTINUITY", flags matches = ( hausdorff is not None and hausdorff <= compare_tolerance ) if raw_valid is True and db_valid is False: if matches: # 极少见：几何几乎一致但 validity 不同，可能构造策略/holes 差异。 return "RAW_CHAIN_VALID_DB_INVALID_NEAR_MATCH", flags return "RAW_CHAIN_VALID_DB_INVALID", flags if raw_valid is False: if matches: return "RAW_EDGE_CHAIN_ITSELF_INVALID_AND_MATCHES_DB", flags return "RAW_EDGE_CHAIN_ITSELF_INVALID_RECONSTRUCTION_DIFFERS", flags if matches: return "RECONSTRUCTION_MATCHES_DB", flags return "RECONSTRUCTION_MISMATCH", flags 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 Area raw edge chain audit v6" ) 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( "--compare-tolerance", type=float, default=2e-6, help="参考重建与 DB geometry Hausdorff distance 容差（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_raw_edge_audit_v6_dataset_{args.dataset_id}_{stamp}" ) outdir.mkdir(parents=True, exist_ok=True) report = { "version": "v6.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 是局部 CM93 数据集；本轮只重建已导入 invalid Area " "自身的 raw edge chain，不要求跨未导入海域连续。" ), "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/6] 找出 52 个 invalid Area ...") invalids = load_invalid_areas(cur, args.dataset_id, wkb_fn) print(f" invalid Area = {len(invalids)}") affected_cells = sorted(set(int(r["cell_id"]) for r in invalids)) print(f" affected Cells = {len(affected_cells)}") print("[2/6] 解析受影响 Cell 的 HEADER / VECTOR ...") cell_cache = {} cell_audit = [] for cell_id in affected_cells: rawtables = load_cell_raw_tables(cur, cell_id) if "HEADER" not in rawtables or "VECTOR" not in rawtables: cell_cache[cell_id] = { "error": "HEADER_OR_VECTOR_MISSING", "rawtables": rawtables, } cell_audit.append({ "cell_id": cell_id, "status": "FAIL", "issue": "HEADER_OR_VECTOR_MISSING", }) continue hrec = normalize_chunk_record(rawtables["HEADER"], "HEADER") vrec = normalize_chunk_record(rawtables["VECTOR"], "VECTOR") try: header = parse_header_metadata(hrec.get("raw_metadata_json")) vector_blob = to_bytes(vrec.get("raw_blob")) edges = parse_vector_table(vector_blob) expected_edges = int(header["usn_vector_records"]) exact = len(edges) == expected_edges cell_cache[cell_id] = { "header": header, "edges": edges, "rawtables": rawtables, } cell_audit.append({ "cell_id": cell_id, "status": "PASS" if exact else "CHECK", "parsed_vector_edges": len(edges), "header_usn_vector_records": expected_edges, "vector_count_match": int(exact), "vector_blob_bytes": len(vector_blob or b""), }) except Exception as e: cell_cache[cell_id] = { "error": f"{type(e).__name__}: {e}", "rawtables": rawtables, } cell_audit.append({ "cell_id": cell_id, "status": "FAIL", "issue": f"{type(e).__name__}: {e}", }) print("[3/6] 解码 invalid Feature raw_feature_blob 并重建 edge rings ...") classification_rows = [] edge_ref_rows = [] ring_rows = [] junction_rows = [] compare_rows = [] isolated_detail = [] # V5 已确认的两条无共享邻接对象 feature，用于特别详细导出。 # 不硬编码成“结论”，只是当前已知调查重点；若 ID 不存在也不影响主流程。 isolated_ids = {24057, 227278} for idx, r in enumerate(invalids, start=1): fid = int(r["feature_id"]) cell_id = int(r["cell_id"]) db_geom = r["_db_geom"] db_valid, db_reason = safe_validity(db_geom) cache = cell_cache.get(cell_id, {}) base = { "feature_id": fid, "object_class": r.get("object_class"), "cell_id": cell_id, "cell_name": r.get("cell_name"), "source_record_index": r.get("source_record_index"), "coverage_id": r.get("coverage_id"), "raw_geotype_db": r.get("raw_geotype"), "db_valid": int(bool(db_valid)) if db_valid is not None else None, "db_validity_reason": db_reason, } if "error" in cache: classification_rows.append({ **base, "final_class": "CELL_RAW_TABLE_UNAVAILABLE", "detail": cache["error"], }) continue raw_feature = to_bytes(r.get("raw_feature_blob")) if not raw_feature: classification_rows.append({ **base, "final_class": "RAW_FEATURE_BLOB_MISSING", }) continue try: parsed = parse_area_feature_blob(raw_feature) except Exception as e: classification_rows.append({ **base, "final_class": "RAW_FEATURE_PARSE_ERROR", "detail": f"{type(e).__name__}: {e}", "raw_feature_hex": raw_feature.hex(), }) continue refs = parsed["refs"] edges = cache["edges"] header = cache["header"] rings_raw, erows, jissues = build_rings_like_opencpn( refs, edges ) # Edge refs export for er in erows: rr = { **base, "feature_geotype_parsed": parsed["geotype"], "feature_object_type": parsed["object_type"], "feature_n_elements": parsed["n_elements"], **er, } edge_ref_rows.append(rr) if fid in isolated_ids: isolated_detail.append(rr) for ji in jissues: junction_rows.append({**base, **ji}) # raw coordinate ring stats repeated_total = 0 for ri, ring in enumerate(rings_raw): reps = repeated_nonconsecutive_vertices(ring) repeated_total += len(reps) raw_line = None try: raw_line = LinearRing(ring) if len(ring) >= 4 else None except Exception: raw_line = None ring_rows.append({ **base, "ring_index": ri, "raw_vertex_count": len(ring), "raw_closed": int(bool(ring) and ring[0] == ring[-1]), "repeated_nonconsecutive_vertex_count": len(reps), "repeated_vertices": reps[:50], "raw_ring_is_simple": ( int(bool(raw_line.is_simple)) if raw_line is not None else None ), "signed_area_raw_units": ( signed_area([(float(x), float(y)) for x, y in ring]) if ring else None ), }) # WGS84 transform wgsox = float(r.get("coverage_wgsox") or 0.0) wgsoy = float(r.get("coverage_wgsoy") or 0.0) rings_wgs = [] for ring in rings_raw: w = [ transform_cm93_point_to_wgs84( p, header, wgsox, wgsoy ) for p in ring ] rings_wgs.append(w) raw_wgs_geom, construction_mode = ring_to_polygon_candidate( rings_wgs ) raw_valid, raw_reason = safe_validity(raw_wgs_geom) haus = safe_hausdorff(raw_wgs_geom, db_geom) if raw_wgs_geom else None symdiff = safe_symdiff_area(raw_wgs_geom, db_geom) if raw_wgs_geom else None final_class, flags = classify_record( junction_issues=jissues, raw_valid=raw_valid, db_valid=bool(db_valid), hausdorff=haus, symdiff=symdiff, compare_tolerance=args.compare_tolerance, ) raw_reason_class = topology_reason_class(raw_reason) db_reason_class = topology_reason_class(db_reason) row = { **base, "parsed_object_type": parsed["object_type"], "parsed_geotype": parsed["geotype"], "parsed_obj_desc_bytes": parsed["obj_desc_bytes"], "parsed_n_elements": parsed["n_elements"], "parsed_has_attributes": parsed["has_attributes"], "vector_edge_count_cell": len(edges), "ring_count_reconstructed": len(rings_raw), "edge_junction_issue_count": len(jissues), "edge_junction_issue_flags": flags, "repeated_nonconsecutive_vertex_total": repeated_total, "coverage_wgsox": wgsox, "coverage_wgsoy": wgsoy, "construction_mode": construction_mode, "raw_reconstructed_valid": ( int(bool(raw_valid)) if raw_valid is not None else None ), "raw_validity_reason": raw_reason, "raw_validity_reason_class": raw_reason_class, "db_validity_reason_class": db_reason_class, "raw_db_same_validity_reason_class": int( raw_reason_class == db_reason_class ), "hausdorff_distance_deg": haus, "symmetric_difference_area_deg2": symdiff, "within_compare_tolerance": ( int(haus is not None and haus <= args.compare_tolerance) ), "final_class": final_class, } classification_rows.append(row) compare_rows.append({ **base, "construction_mode": construction_mode, "raw_geom_type": ( raw_wgs_geom.geom_type if raw_wgs_geom is not None else None ), "db_geom_type": db_geom.geom_type, "raw_valid": ( int(bool(raw_valid)) if raw_valid is not None else None ), "db_valid": int(bool(db_valid)), "raw_reason": raw_reason, "db_reason": db_reason, "hausdorff_distance_deg": haus, "symmetric_difference_area_deg2": symdiff, "within_compare_tolerance": ( int(haus is not None and haus <= args.compare_tolerance) ), }) # 挂到 GeoJSON 使用。 r["_raw_wgs_geom"] = raw_wgs_geom r["_raw_wgs_reason"] = raw_reason r["_construction_mode"] = construction_mode r["_final_class"] = final_class print( f" [{idx:02d}/{len(invalids)}] " f"feature={fid} {r.get('object_class')} " f"edges={len(refs)} rings={len(rings_raw)} " f"class={final_class}" ) print("[4/6] 汇总 raw vs DB 证据 ...") class_counts = Counter( r.get("final_class", "UNKNOWN") for r in classification_rows ) raw_reason_counts = Counter( r.get("raw_validity_reason_class", "UNKNOWN") for r in classification_rows if r.get("raw_validity_reason_class") ) db_reason_counts = Counter( r.get("db_validity_reason_class", "UNKNOWN") for r in classification_rows if r.get("db_validity_reason_class") ) vector_bad_cells = sum( 1 for r in cell_audit if r.get("vector_count_match") == 0 or r.get("status") == "FAIL" ) report["summary"] = { "invalid_area_count": len(invalids), "affected_cell_count": len(affected_cells), "cell_vector_audit_bad_count": vector_bad_cells, "classification_counts": dict(class_counts), "raw_validity_reason_counts": dict(raw_reason_counts), "db_validity_reason_counts": dict(db_reason_counts), "edge_junction_issue_count": len(junction_rows), "features_with_edge_junction_issue": len(set( int(r["feature_id"]) for r in junction_rows if r.get("feature_id") is not None )), "reference_reconstruction_within_tolerance": sum( int(r.get("within_compare_tolerance") or 0) for r in classification_rows ), "same_validity_reason_class_count": sum( int(r.get("raw_db_same_validity_reason_class") or 0) for r in classification_rows ), } print("[5/6] 写 CSV / JSON / GeoJSON ...") write_csv( outdir / "invalid_area_raw_edge_classification.csv", classification_rows ) write_csv( outdir / "invalid_area_edge_references.csv", edge_ref_rows ) write_csv( outdir / "invalid_area_ring_summary.csv", ring_rows ) write_csv( outdir / "invalid_area_edge_junction_issues.csv", junction_rows ) write_csv( outdir / "raw_vs_db_geometry_compare.csv", compare_rows ) write_csv( outdir / "isolated_invalid_area_edge_detail.csv", isolated_detail ) write_csv( outdir / "affected_cell_vector_table_audit.csv", cell_audit ) make_geojson( [ r for r in invalids if "_raw_wgs_geom" in r ], outdir / "invalid_area_raw_reconstruction.geojson" ) (outdir / "raw_edge_audit_report_v6.json").write_text( json.dumps(report, ensure_ascii=False, indent=2, default=str), encoding="utf-8" ) lines = [] lines.append("CM93 -> MySQL invalid Area 原始 Edge Chain 重建审计 v6") lines.append("=" * 88) lines.append(f"Generated: {report['generated_at']}") lines.append(f"Dataset ID: {args.dataset_id}") lines.append(f"MySQL: {report['mysql_version']}") lines.append("") lines.append("SCOPE") lines.append( " 当前 CM93JAPAN 是局部数据集；不要求与未导入海域跨边界连续。" ) lines.append( " 本轮直接从 DB 保存的 raw_feature_blob + VECTOR raw_blob " "重建 invalid Area 的 CM93 edge chain。" ) lines.append("") lines.append("BASELINE") lines.append(f" invalid Area: {len(invalids)}") lines.append(f" affected Cells: {len(affected_cells)}") lines.append(f" bad VECTOR-table audits: {vector_bad_cells}") lines.append("") lines.append("FINAL CLASSIFICATION") for k, v in class_counts.most_common(): lines.append(f" {k:56s} {v}") lines.append("") lines.append("RAW RECONSTRUCTED VALIDITY") for k, v in raw_reason_counts.most_common(): lines.append(f" {k:40s} {v}") lines.append("") lines.append("DB VALIDITY") for k, v in db_reason_counts.most_common(): lines.append(f" {k:40s} {v}") lines.append("") lines.append("EDGE CHAIN") lines.append( f" junction issue rows: {len(junction_rows)}" ) lines.append( f" features with junction issue:" f" {report['summary']['features_with_edge_junction_issue']}" ) lines.append("") lines.append("RAW ↔ DB") lines.append( f" within Hausdorff tolerance: " f"{report['summary']['reference_reconstruction_within_tolerance']}" f" / {len(classification_rows)}" ) lines.append( f" same validity reason class: " f"{report['summary']['same_validity_reason_class_count']}" f" / {len(classification_rows)}" ) lines.append("") lines.append("INTERPRETATION") lines.append( " RAW_CHAIN_VALID_DB_INVALID：优先判为 importer assembly/conversion 候选。" ) lines.append( " RAW_EDGE_CHAIN_ITSELF_INVALID_AND_MATCHES_DB：" "说明 DB 基本忠实反映了 CM93 edge chain 组装结果；" "后续重点研究 CM93 shared topology / OGC Polygon validity 语义，" "而不是先修改坐标。" ) lines.append( " EDGE_CHAIN_DISCONTINUITY：检查 edge index / usage / ordering；" "不得自动闭环或补边。" ) lines.append( " RECONSTRUCTION_MISMATCH：检查 coverage WGS offset、ring shell/hole " "解释以及 importer 的坐标转换。" ) lines.append("") lines.append("NO WRITE / NO REPAIR / NO PBF") (outdir / "raw_edge_audit_summary_v6.txt").write_text( "\n".join(lines), encoding="utf-8" ) conn.rollback() print("[6/6] 压缩输出 ...") 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()