#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ CM93 -> MySQL Dataset ID=1 实体空间 Geometry 只读审计 Version: v3.0 目的： 1. 审计已经真正写入 MySQL 的 Point / Line / Area / Sounding / Coverage geometry。 2. 不再把 raw_geotype=16 的正常 META parent 当作 geometry 错误。 3. 检查单个 geometry 的“内在有效性”，不要求当前日本局部 Dataset 与未导入海域连续。 4. 不修改数据库，不修 geometry，不做 PBF。 特别边界规则： 当前 CM93JAPAN 只是 CM93 的日本沿海局部子集，并非完整 CM93。 因此： - 不做跨 Dataset 外边界连续性 FAIL 判定； - 不因 Line 在 Cell/Dataset 边缘终止而判坏； - 不因缺少未导入邻区对象而判坏； - 对靠近 Cell/Dataset 边界的异常只增加 boundary_context 提示，不自动归因于 importer。 依赖： pip install pymysql shapely 示例： python cm93_mysql_spatial_audit_v3.py ^ --host 127.0.0.1 ^ --port 3306 ^ --user root ^ --database cm93_chart ^ --dataset-id 1 密码缺省时交互输入。 输出： - spatial_audit_summary.txt - spatial_audit_report.json - feature_geometry_mapping_issues.csv - intrinsic_geometry_issues.csv - intrinsic_issue_category_counts.csv - invalid_geometry_reason_counts.csv - geometry_stats_by_kind.csv - issue_object_class_counts.csv - issue_cell_counts.csv - issue_samples.csv - sounding_feature_bbox_issues.csv - coverage_geometry_issues.csv 并自动生成同名 zip。 安全： - SQL 仅包含 SELECT / SHOW / SET TRANSACTION READ ONLY - 无 INSERT / UPDATE / DELETE / DDL """ from __future__ import annotations import argparse import csv import getpass import json import math import os import sys import traceback import zipfile from collections import Counter, defaultdict from dataclasses import dataclass, asdict from datetime import datetime from pathlib import Path from typing import Any, Dict, Iterable, Iterator, 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.validation import explain_validity except Exception: def explain_validity(g): return "invalid" except ImportError: print("缺少 shapely。请执行：pip install pymysql shapely") sys.exit(2) EXPECTED_DATASET1 = { "cells": 784, "features": 322465, "attributes": 405153, "sounding_points": 275700, } EXPECTED_GEOMETRY_ISSUE_INTERPRETATION = { "raw_geotype_16_expected_parent_no_geometry": 4406, } def qi(name: str) -> str: return "`" + name.replace("`", "``") + "`" def finite(v: Optional[float]) -> bool: return v is not None and math.isfinite(float(v)) def safe_float(v: Any) -> Optional[float]: if v is None: return None try: x = float(v) return x if math.isfinite(x) else None except Exception: return None def bbox_tuple(row: Dict[str, Any], prefix: str = "") -> Optional[Tuple[float, float, float, float]]: vals = [ safe_float(row.get(prefix + "bbox_min_lon")), safe_float(row.get(prefix + "bbox_min_lat")), safe_float(row.get(prefix + "bbox_max_lon")), safe_float(row.get(prefix + "bbox_max_lat")), ] if any(v is None for v in vals): return None return vals[0], vals[1], vals[2], vals[3] def bbox_close(a: Tuple[float, float, float, float], b: Tuple[float, float, float, float], tol: float) -> bool: return all(abs(float(x) - float(y)) <= tol for x, y in zip(a, b)) def touches_bbox_edge(inner: Optional[Tuple[float, float, float, float]], outer: Optional[Tuple[float, float, float, float]], tol: float) -> bool: if not inner or not outer: return False return ( abs(inner[0] - outer[0]) <= tol or abs(inner[1] - outer[1]) <= tol or abs(inner[2] - outer[2]) <= tol or abs(inner[3] - outer[3]) <= tol ) def bbox_in_lonlat_range(b: Optional[Tuple[float, float, float, float]]) -> bool: if not b: return False minx, miny, maxx, maxy = b return ( -180.0 <= minx <= 180.0 and -180.0 <= maxx <= 180.0 and -90.0 <= miny <= 90.0 and -90.0 <= maxy <= 90.0 and minx <= maxx and miny <= maxy ) def write_csv(path: Path, rows: Iterable[Dict[str, Any]], fieldnames: Optional[List[str]] = None): rows = list(rows) path.parent.mkdir(parents=True, exist_ok=True) 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) with path.open("w", encoding="utf-8-sig", newline="") as f: w = csv.DictWriter(f, fieldnames=fieldnames or [], extrasaction="ignore") if fieldnames: w.writeheader() for r in rows: w.writerow({k: "" if r.get(k) is None else r.get(k) for k in fieldnames}) def scalar(cur, sql: str, params: Sequence[Any] = ()) -> Any: cur.execute(sql, params) row = cur.fetchone() if not row: return None return next(iter(row.values())) def rows(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 不支持/失败：{e}") try: cur.execute("START TRANSACTION READ ONLY") except Exception as e: warnings.append(f"START TRANSACTION READ ONLY 不支持/失败：{e}") try: cur.execute("START TRANSACTION") except Exception: pass def probe_spatial_functions(cur) -> Dict[str, Optional[str]]: """ 探测把 MySQL/MariaDB geometry 转成标准 WKB 的函数和 SRID 函数。 不使用 ST_IsValid，实际 validity 统一交给 Shapely，避免 DB 方言差异。 """ result = {"wkb_fn": None, "srid_fn": None} for fn in ["ST_AsBinary", "ST_AsWKB", "AsBinary", "AsWKB"]: try: cur.execute( f"SELECT {fn}(geom) AS x " f"FROM chart_point_geometry LIMIT 1" ) cur.fetchone() result["wkb_fn"] = fn break except Exception: continue for fn in ["ST_SRID", "SRID"]: try: cur.execute( f"SELECT {fn}(geom) AS x " f"FROM chart_point_geometry LIMIT 1" ) cur.fetchone() result["srid_fn"] = fn break except Exception: continue return result def classify_invalid_reason(reason: str) -> str: s = (reason or "").lower() if "self-intersection" in s: if "ring self-intersection" in s: return "RING_SELF_INTERSECTION" return "SELF_INTERSECTION" if "too few points" in s: return "TOO_FEW_POINTS" 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 "disconnected interior" in s: return "DISCONNECTED_INTERIOR" if "invalid coordinate" in s or "nan" in s or "inf" in s: return "INVALID_COORDINATE" if "ring" in s and "closed" in s: return "RING_NOT_CLOSED" return "OTHER_INVALID_GEOMETRY" def actual_type_ok(kind: str, geom_type: str) -> bool: if kind == "POINT": return geom_type == "Point" if kind == "LINE": return geom_type in {"LineString", "MultiLineString"} if kind == "AREA": return geom_type in {"Polygon", "MultiPolygon"} if kind == "SOUNDING": return geom_type == "Point" if kind == "COVERAGE": return geom_type in {"Polygon", "MultiPolygon"} return True def intrinsic_check( *, kind: str, feature_id: Optional[int], object_class: str, cell_id: Optional[int], cell_name: str, source_record_index: Optional[int], wkb_data: Any, srid: Any, normalized_srid: Any, feature_bbox: Optional[Tuple[float, float, float, float]], cell_bbox: Optional[Tuple[float, float, float, float]], dataset_bbox: Optional[Tuple[float, float, float, float]], bbox_tolerance: float, boundary_tolerance: float, extra_id: Optional[int] = None, ) -> Tuple[Dict[str, Any], List[Dict[str, Any]]]: """ 返回： stat: 单 geometry 统计 issues: 0..N 条 issue """ base = { "kind": kind, "feature_id": feature_id, "extra_id": extra_id, "object_class": object_class, "cell_id": cell_id, "cell_name": cell_name, "source_record_index": source_record_index, "srid": srid, "normalized_srid": normalized_srid, } stat = dict(base) stat.update({ "parsed": 0, "empty": 0, "valid": None, "geom_type": None, "simple": None, "bounds": None, "touches_cell_boundary": 0, "touches_dataset_bbox": 0, }) issues: List[Dict[str, Any]] = [] if wkb_data is None: x = dict(base) x.update({ "issue_category": "WKB_MISSING", "issue_detail": "DB geometry exists but WKB conversion returned NULL", "invalid_reason": "", "boundary_context": "", }) issues.append(x) return stat, issues if isinstance(wkb_data, memoryview): wkb_data = wkb_data.tobytes() if isinstance(wkb_data, bytearray): wkb_data = bytes(wkb_data) try: geom = load_wkb(wkb_data) stat["parsed"] = 1 except Exception as e: x = dict(base) x.update({ "issue_category": "WKB_PARSE_ERROR", "issue_detail": f"{type(e).__name__}: {e}", "invalid_reason": "", "boundary_context": "", }) issues.append(x) return stat, issues geom_type = getattr(geom, "geom_type", None) stat["geom_type"] = geom_type try: is_empty = bool(geom.is_empty) except Exception: is_empty = False stat["empty"] = int(is_empty) if is_empty: x = dict(base) x.update({ "issue_category": "EMPTY_GEOMETRY", "issue_detail": f"{kind} geometry is EMPTY", "invalid_reason": "", "actual_geom_type": geom_type, "boundary_context": "", }) issues.append(x) return stat, issues try: b = tuple(float(v) for v in geom.bounds) if len(b) == 4 and all(math.isfinite(v) for v in b): stat["bounds"] = b else: b = None except Exception: b = None stat["touches_cell_boundary"] = int(touches_bbox_edge(b, cell_bbox, boundary_tolerance)) stat["touches_dataset_bbox"] = int(touches_bbox_edge(b, dataset_bbox, boundary_tolerance)) boundary_bits = [] if stat["touches_cell_boundary"]: boundary_bits.append("CELL_EDGE") if stat["touches_dataset_bbox"]: boundary_bits.append("DATASET_BBOX_EDGE") boundary_context = "|".join(boundary_bits) if not actual_type_ok(kind, geom_type or ""): x = dict(base) x.update({ "issue_category": "GEOMETRY_TYPE_MISMATCH", "issue_detail": f"expected={kind}, actual={geom_type}", "invalid_reason": "", "actual_geom_type": geom_type, "boundary_context": boundary_context, }) issues.append(x) # 坐标范围：CM93 normalized SRID 当前应为经纬度数据。 if b is None: x = dict(base) x.update({ "issue_category": "NON_FINITE_OR_UNREADABLE_BOUNDS", "issue_detail": "geometry bounds unavailable/non-finite", "invalid_reason": "", "actual_geom_type": geom_type, "boundary_context": boundary_context, }) issues.append(x) elif not bbox_in_lonlat_range(b): x = dict(base) x.update({ "issue_category": "OUT_OF_RANGE_COORDINATE", "issue_detail": f"bounds={b}", "invalid_reason": "", "actual_geom_type": geom_type, "boundary_context": boundary_context, }) issues.append(x) # SRID 一致性。无法探测 SRID 函数时 srid=None，不判 FAIL。 if srid is not None and normalized_srid is not None: try: if int(srid) != int(normalized_srid): x = dict(base) x.update({ "issue_category": "SRID_MISMATCH", "issue_detail": f"geom_srid={srid}, feature_normalized_srid={normalized_srid}", "invalid_reason": "", "actual_geom_type": geom_type, "boundary_context": boundary_context, }) issues.append(x) except Exception: pass # Feature bbox 对账。Sounding 在后面按 feature 聚合单独核对。 if kind in {"POINT", "LINE", "AREA"}: if feature_bbox is None: x = dict(base) x.update({ "issue_category": "FEATURE_BBOX_MISSING", "issue_detail": f"{kind} has geometry but feature bbox is NULL/incomplete", "invalid_reason": "", "actual_geom_type": geom_type, "geometry_bounds": b, "boundary_context": boundary_context, }) issues.append(x) elif b is not None and not bbox_close(b, feature_bbox, bbox_tolerance): x = dict(base) x.update({ "issue_category": "FEATURE_BBOX_MISMATCH", "issue_detail": f"feature_bbox={feature_bbox}, geometry_bounds={b}, tol={bbox_tolerance}", "invalid_reason": "", "actual_geom_type": geom_type, "geometry_bounds": b, "feature_bbox": feature_bbox, "boundary_context": boundary_context, }) issues.append(x) # Shapely validity： # - Area/Coverage 的 invalid 是核心检查； # - Point/Line 也记录 validity； # - Line 非 simple 不等同 invalid，只作为统计，不生成 issue。 try: valid = bool(geom.is_valid) stat["valid"] = int(valid) except Exception: valid = True stat["valid"] = None try: stat["simple"] = int(bool(geom.is_simple)) except Exception: stat["simple"] = None if not valid: try: reason = explain_validity(geom) except Exception as e: reason = f"validity explanation unavailable: {e}" reason_class = classify_invalid_reason(reason) x = dict(base) x.update({ "issue_category": reason_class, "issue_detail": reason, "invalid_reason": reason, "actual_geom_type": geom_type, "geometry_bounds": b, "boundary_context": boundary_context, }) issues.append(x) return stat, issues def get_dataset_bbox(cur, dataset_id: int) -> Optional[Tuple[float, float, float, float]]: cur.execute(""" SELECT MIN(bbox_min_lon) AS minx, MIN(bbox_min_lat) AS miny, MAX(bbox_max_lon) AS maxx, MAX(bbox_max_lat) AS maxy FROM chart_cell WHERE dataset_id=%s AND bbox_min_lon IS NOT NULL AND bbox_min_lat IS NOT NULL AND bbox_max_lon IS NOT NULL AND bbox_max_lat IS NOT NULL """, (dataset_id,)) r = cur.fetchone() or {} vals = [safe_float(r.get("minx")), safe_float(r.get("miny")), safe_float(r.get("maxx")), safe_float(r.get("maxy"))] if any(v is None for v in vals): return None return vals[0], vals[1], vals[2], vals[3] def audit_mapping(cur, dataset_id: int) -> Tuple[List[Dict[str, Any]], Dict[str, Any]]: """ 检查 feature.geometry_type 与实际存储表是否一致。 """ mapping_issues: List[Dict[str, Any]] = [] stats: Dict[str, Any] = {} # Feature type counts ft = rows(cur, """ SELECT geometry_type, COUNT(*) AS n FROM chart_feature WHERE dataset_id=%s GROUP BY geometry_type ORDER BY geometry_type """, (dataset_id,)) stats["feature_type_counts"] = ft # 实体表 counts（只统计 Dataset 1） for name, table in [ ("point_rows", "chart_point_geometry"), ("line_rows", "chart_line_geometry"), ("area_rows", "chart_area_geometry"), ]: stats[name] = int(scalar(cur, f""" SELECT COUNT(*) FROM {table} g JOIN chart_feature f ON f.id=g.feature_id WHERE f.dataset_id=%s """, (dataset_id,)) or 0) stats["sounding_rows"] = int(scalar(cur, """ SELECT COUNT(*) FROM chart_sounding_point s JOIN chart_feature f ON f.id=s.feature_id WHERE f.dataset_id=%s """, (dataset_id,)) or 0) # 缺 expected geometry checks = [ ("POINT_MISSING_GEOMETRY", "POINT", "chart_point_geometry"), ("LINE_MISSING_GEOMETRY", "LINE", "chart_line_geometry"), ("AREA_MISSING_GEOMETRY", "AREA", "chart_area_geometry"), ] for issue_cat, gtype, table in checks: rs = rows(cur, f""" SELECT f.id AS feature_id, f.object_class, f.cell_id, c.cell_name, f.source_record_index, f.raw_geotype, f.geometry_type FROM chart_feature f JOIN chart_cell c ON c.id=f.cell_id LEFT JOIN {table} g ON g.feature_id=f.id WHERE f.dataset_id=%s AND f.geometry_type=%s AND g.feature_id IS NULL ORDER BY f.id """, (dataset_id, gtype)) for r in rs: r["issue_category"] = issue_cat mapping_issues.append(r) # SOUNDING feature 没有任何 sounding point rs = rows(cur, """ SELECT f.id AS feature_id, f.object_class, f.cell_id, c.cell_name, f.source_record_index, f.raw_geotype, f.geometry_type FROM chart_feature f JOIN chart_cell c ON c.id=f.cell_id LEFT JOIN chart_sounding_point s ON s.feature_id=f.id WHERE f.dataset_id=%s AND f.geometry_type='SOUNDING' GROUP BY f.id, f.object_class, f.cell_id, c.cell_name, f.source_record_index, f.raw_geotype, f.geometry_type HAVING COUNT(s.id)=0 ORDER BY f.id """, (dataset_id,)) for r in rs: r["issue_category"] = "SOUNDING_MISSING_POINTS" mapping_issues.append(r) # geometry 表串类型 for issue_cat, table, expected in [ ("POINT_TABLE_WRONG_FEATURE_TYPE", "chart_point_geometry", "POINT"), ("LINE_TABLE_WRONG_FEATURE_TYPE", "chart_line_geometry", "LINE"), ("AREA_TABLE_WRONG_FEATURE_TYPE", "chart_area_geometry", "AREA"), ]: rs = rows(cur, f""" SELECT f.id AS feature_id, f.object_class, f.cell_id, c.cell_name, f.source_record_index, f.raw_geotype, f.geometry_type FROM {table} g JOIN chart_feature f ON f.id=g.feature_id JOIN chart_cell c ON c.id=f.cell_id WHERE f.dataset_id=%s AND f.geometry_type<>%s ORDER BY f.id """, (dataset_id, expected)) for r in rs: r["issue_category"] = issue_cat mapping_issues.append(r) # sounding point 串到非 SOUNDING feature rs = rows(cur, """ SELECT DISTINCT f.id AS feature_id, f.object_class, f.cell_id, c.cell_name, f.source_record_index, f.raw_geotype, f.geometry_type FROM chart_sounding_point s JOIN chart_feature f ON f.id=s.feature_id JOIN chart_cell c ON c.id=f.cell_id WHERE f.dataset_id=%s AND f.geometry_type<>'SOUNDING' ORDER BY f.id """, (dataset_id,)) for r in rs: r["issue_category"] = "SOUNDING_POINT_WRONG_FEATURE_TYPE" mapping_issues.append(r) # META/UNKNOWN 不应偷偷挂实际 geometry rs = rows(cur, """ SELECT f.id AS feature_id, f.object_class, f.cell_id, c.cell_name, f.source_record_index, f.raw_geotype, f.geometry_type, CASE WHEN pg.feature_id IS NOT NULL THEN 1 ELSE 0 END AS has_point, CASE WHEN lg.feature_id IS NOT NULL THEN 1 ELSE 0 END AS has_line, CASE WHEN ag.feature_id IS NOT NULL THEN 1 ELSE 0 END AS has_area, COUNT(sp.id) AS sounding_point_count FROM chart_feature f JOIN chart_cell c ON c.id=f.cell_id LEFT JOIN chart_point_geometry pg ON pg.feature_id=f.id LEFT JOIN chart_line_geometry lg ON lg.feature_id=f.id LEFT JOIN chart_area_geometry ag ON ag.feature_id=f.id LEFT JOIN chart_sounding_point sp ON sp.feature_id=f.id WHERE f.dataset_id=%s AND f.geometry_type IN ('META','UNKNOWN') GROUP BY f.id, f.object_class, f.cell_id, c.cell_name, f.source_record_index, f.raw_geotype, f.geometry_type, pg.feature_id, lg.feature_id, ag.feature_id HAVING has_point=1 OR has_line=1 OR has_area=1 OR sounding_point_count>0 ORDER BY f.id """, (dataset_id,)) for r in rs: r["issue_category"] = "META_OR_UNKNOWN_HAS_UNEXPECTED_GEOMETRY" mapping_issues.append(r) stats["mapping_issue_count"] = len(mapping_issues) stats["mapping_issue_categories"] = dict(Counter(r["issue_category"] for r in mapping_issues)) return mapping_issues, stats def scan_feature_geometry_table( cur, *, table: str, kind: str, dataset_id: int, wkb_fn: str, srid_fn: Optional[str], dataset_bbox: Optional[Tuple[float, float, float, float]], chunk_size: int, bbox_tolerance: float, boundary_tolerance: float, ) -> Tuple[Dict[str, Any], List[Dict[str, Any]]]: last_id = 0 stats = Counter() geom_types = Counter() invalid_reasons = Counter() issues: List[Dict[str, Any]] = [] srid_expr = f"{srid_fn}(g.geom)" if srid_fn else "NULL" while True: sql = f""" SELECT g.feature_id, {wkb_fn}(g.geom) AS geom_wkb, {srid_expr} AS geom_srid, f.object_class, f.cell_id, c.cell_name, f.source_record_index, f.normalized_srid, 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 FROM {table} g JOIN chart_feature f ON f.id=g.feature_id JOIN chart_cell c ON c.id=f.cell_id WHERE f.dataset_id=%s AND g.feature_id>%s ORDER BY g.feature_id LIMIT %s """ cur.execute(sql, (dataset_id, last_id, chunk_size)) batch = list(cur.fetchall()) if not batch: break for r in batch: last_id = int(r["feature_id"]) fb = bbox_tuple(r) cb = bbox_tuple(r, "cell_") stat, found = intrinsic_check( kind=kind, feature_id=r["feature_id"], object_class=r.get("object_class") or "", cell_id=r.get("cell_id"), cell_name=r.get("cell_name") or "", source_record_index=r.get("source_record_index"), wkb_data=r.get("geom_wkb"), srid=r.get("geom_srid"), normalized_srid=r.get("normalized_srid"), feature_bbox=fb, cell_bbox=cb, dataset_bbox=dataset_bbox, bbox_tolerance=bbox_tolerance, boundary_tolerance=boundary_tolerance, ) stats["rows"] += 1 stats["parsed"] += stat["parsed"] stats["empty"] += stat["empty"] if stat["valid"] == 1: stats["valid"] += 1 elif stat["valid"] == 0: stats["invalid"] += 1 if kind == "LINE" and stat["simple"] == 0: stats["non_simple_line_info"] += 1 if stat["touches_cell_boundary"]: stats["touch_cell_boundary_info"] += 1 if stat["touches_dataset_bbox"]: stats["touch_dataset_bbox_info"] += 1 geom_types[stat.get("geom_type") or "UNPARSED"] += 1 for x in found: issues.append(x) invalid_reasons[x["issue_category"]] += 1 print( f" {kind:<8} scanned={stats['rows']:,} " f"issues={len(issues):,} last_feature_id={last_id}" ) return { "kind": kind, "table": table, "rows": stats["rows"], "parsed": stats["parsed"], "empty": stats["empty"], "valid": stats["valid"], "invalid": stats["invalid"], "non_simple_line_info": stats["non_simple_line_info"], "touch_cell_boundary_info": stats["touch_cell_boundary_info"], "touch_dataset_bbox_info": stats["touch_dataset_bbox_info"], "actual_geom_types": dict(geom_types), "issue_count": len(issues), "issue_categories": dict(invalid_reasons), }, issues def scan_soundings( cur, *, dataset_id: int, wkb_fn: str, srid_fn: Optional[str], dataset_bbox: Optional[Tuple[float, float, float, float]], chunk_size: int, bbox_tolerance: float, boundary_tolerance: float, ) -> Tuple[Dict[str, Any], List[Dict[str, Any]], List[Dict[str, Any]]]: last_id = 0 stats = Counter() geom_types = Counter() issues: List[Dict[str, Any]] = [] # feature_id -> [minx,miny,maxx,maxy,count] feature_bounds: Dict[int, List[float]] = {} srid_expr = f"{srid_fn}(s.geom)" if srid_fn else "NULL" while True: sql = f""" SELECT s.id AS sounding_id, s.feature_id, s.depth, {wkb_fn}(s.geom) AS geom_wkb, {srid_expr} AS geom_srid, f.object_class, f.cell_id, c.cell_name, f.source_record_index, f.normalized_srid, 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 FROM chart_sounding_point s JOIN chart_feature f ON f.id=s.feature_id JOIN chart_cell c ON c.id=f.cell_id WHERE f.dataset_id=%s AND s.id>%s ORDER BY s.id LIMIT %s """ cur.execute(sql, (dataset_id, last_id, chunk_size)) batch = list(cur.fetchall()) if not batch: break for r in batch: last_id = int(r["sounding_id"]) fb = bbox_tuple(r) cb = bbox_tuple(r, "cell_") stat, found = intrinsic_check( kind="SOUNDING", feature_id=r["feature_id"], extra_id=r["sounding_id"], object_class=r.get("object_class") or "", cell_id=r.get("cell_id"), cell_name=r.get("cell_name") or "", source_record_index=r.get("source_record_index"), wkb_data=r.get("geom_wkb"), srid=r.get("geom_srid"), normalized_srid=r.get("normalized_srid"), feature_bbox=fb, cell_bbox=cb, dataset_bbox=dataset_bbox, bbox_tolerance=bbox_tolerance, boundary_tolerance=boundary_tolerance, ) stats["rows"] += 1 stats["parsed"] += stat["parsed"] stats["empty"] += stat["empty"] if stat["valid"] == 1: stats["valid"] += 1 elif stat["valid"] == 0: stats["invalid"] += 1 geom_types[stat.get("geom_type") or "UNPARSED"] += 1 # depth 必须 finite depth = safe_float(r.get("depth")) if depth is None: issues.append({ "kind": "SOUNDING", "feature_id": r["feature_id"], "extra_id": r["sounding_id"], "object_class": r.get("object_class") or "", "cell_id": r.get("cell_id"), "cell_name": r.get("cell_name") or "", "source_record_index": r.get("source_record_index"), "issue_category": "INVALID_SOUNDING_DEPTH", "issue_detail": f"depth={r.get('depth')}", "boundary_context": "", }) issues.extend(found) b = stat.get("bounds") if b and len(b) == 4: fid = int(r["feature_id"]) if fid not in feature_bounds: feature_bounds[fid] = [b[0], b[1], b[2], b[3], 1] else: a = feature_bounds[fid] a[0] = min(a[0], b[0]) a[1] = min(a[1], b[1]) a[2] = max(a[2], b[2]) a[3] = max(a[3], b[3]) a[4] += 1 print( f" SOUNDING scanned={stats['rows']:,} " f"issues={len(issues):,} last_sounding_id={last_id}" ) # 按 feature 聚合 sounding bbox 与 chart_feature bbox 对账。 bbox_issues: List[Dict[str, Any]] = [] feature_ids = list(feature_bounds.keys()) # 704 个 feature 左右，可直接分批 IN for pos in range(0, len(feature_ids), 500): ids = feature_ids[pos:pos+500] placeholders = ",".join(["%s"] * len(ids)) rs = rows(cur, f""" SELECT f.id AS feature_id, f.object_class, f.cell_id, c.cell_name, f.source_record_index, 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 FROM chart_feature f JOIN chart_cell c ON c.id=f.cell_id WHERE f.id IN ({placeholders}) """, ids) for r in rs: fid = int(r["feature_id"]) a = feature_bounds[fid] actual = (a[0], a[1], a[2], a[3]) fb = bbox_tuple(r) cb = bbox_tuple(r, "cell_") context = [] if touches_bbox_edge(actual, cb, boundary_tolerance): context.append("CELL_EDGE") if touches_bbox_edge(actual, dataset_bbox, boundary_tolerance): context.append("DATASET_BBOX_EDGE") if fb is None: bbox_issues.append({ "feature_id": fid, "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"), "issue_category": "SOUNDING_FEATURE_BBOX_MISSING", "point_count": int(a[4]), "actual_bbox": actual, "feature_bbox": None, "boundary_context": "|".join(context), }) elif not bbox_close(actual, fb, bbox_tolerance): bbox_issues.append({ "feature_id": fid, "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"), "issue_category": "SOUNDING_FEATURE_BBOX_MISMATCH", "point_count": int(a[4]), "actual_bbox": actual, "feature_bbox": fb, "boundary_context": "|".join(context), }) return { "kind": "SOUNDING", "table": "chart_sounding_point", "rows": stats["rows"], "parsed": stats["parsed"], "empty": stats["empty"], "valid": stats["valid"], "invalid": stats["invalid"], "actual_geom_types": dict(geom_types), "issue_count_excluding_feature_bbox": len(issues), "feature_bbox_issue_count": len(bbox_issues), }, issues, bbox_issues def scan_coverage( cur, *, dataset_id: int, wkb_fn: str, srid_fn: Optional[str], dataset_bbox: Optional[Tuple[float, float, float, float]], bbox_tolerance: float, boundary_tolerance: float, ) -> Tuple[Dict[str, Any], List[Dict[str, Any]]]: srid_expr = f"{srid_fn}(v.geom)" if srid_fn else "NULL" cur.execute(f""" SELECT v.id AS coverage_id, v.cell_id, c.cell_name, v.source_record_index, {wkb_fn}(v.geom) AS geom_wkb, {srid_expr} AS geom_srid, 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 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,)) rs = list(cur.fetchall()) stats = Counter() geom_types = Counter() issues: List[Dict[str, Any]] = [] for r in rs: cb = bbox_tuple(r, "cell_") stat, found = intrinsic_check( kind="COVERAGE", feature_id=None, extra_id=r["coverage_id"], object_class="COVERAGE", cell_id=r.get("cell_id"), cell_name=r.get("cell_name") or "", source_record_index=r.get("source_record_index"), wkb_data=r.get("geom_wkb"), srid=r.get("geom_srid"), normalized_srid=None, feature_bbox=None, cell_bbox=cb, dataset_bbox=dataset_bbox, bbox_tolerance=bbox_tolerance, boundary_tolerance=boundary_tolerance, ) stats["rows"] += 1 stats["parsed"] += stat["parsed"] stats["empty"] += stat["empty"] if stat["valid"] == 1: stats["valid"] += 1 elif stat["valid"] == 0: stats["invalid"] += 1 geom_types[stat.get("geom_type") or "UNPARSED"] += 1 issues.extend(found) return { "kind": "COVERAGE", "table": "chart_coverage", "rows": stats["rows"], "parsed": stats["parsed"], "empty": stats["empty"], "valid": stats["valid"], "invalid": stats["invalid"], "actual_geom_types": dict(geom_types), "issue_count": len(issues), }, issues def zip_directory(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 MySQL 实体空间 geometry 只读审计 v3" ) 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("--chunk-size", type=int, default=2000) p.add_argument( "--bbox-tolerance", type=float, default=1e-6, help="Feature bbox 与实际 geometry bounds 对账容差（degree），默认 1e-6" ) p.add_argument( "--boundary-tolerance", type=float, default=1e-5, help="仅用于标注 geometry 是否碰 Cell/Dataset bbox 边缘，不用于 PASS/FAIL" ) p.add_argument("--output", default=None) p.add_argument( "--skip-coverage", action="store_true", help="跳过 chart_coverage geometry 审计" ) 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_spatial_audit_v3_dataset_{args.dataset_id}_{stamp}" ) outdir.mkdir(parents=True, exist_ok=True) report: Dict[str, Any] = { "version": "v3.0", "generated_at": datetime.now().isoformat(timespec="seconds"), "dataset_id": args.dataset_id, "database": args.database, "read_only": True, "pbf": "OUT_OF_SCOPE", "dataset_scope_note": ( "当前 CM93JAPAN 只是日本沿海局部子集，并非完整 CM93。" "不把跨未导入海域/邻区的不连续性作为 geometry corruption。" ), "warnings": [], "server": {}, "dataset_bbox": None, "counts": {}, "mapping": {}, "geometry_stats": [], "issue_summary": {}, } 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"]) print("[1/8] 读取服务器与 Dataset 基线...") report["server"]["version"] = scalar(cur, "SELECT VERSION()") try: report["server"]["version_comment"] = scalar(cur, "SELECT @@version_comment") except Exception: report["server"]["version_comment"] = None funcs = probe_spatial_functions(cur) report["server"]["spatial_functions"] = funcs if not funcs["wkb_fn"]: raise RuntimeError( "无法找到可用的 geometry -> WKB 函数 " "(ST_AsBinary/ST_AsWKB/AsBinary/AsWKB)" ) print(f" DB: {report['server']['version']}") print(f" WKB function: {funcs['wkb_fn']}") print(f" SRID function: {funcs['srid_fn']}") dataset_bbox = get_dataset_bbox(cur, args.dataset_id) report["dataset_bbox"] = dataset_bbox print(f" Dataset bbox: {dataset_bbox}") # 基线 counts report["counts"]["cells"] = int(scalar(cur, "SELECT COUNT(*) FROM chart_cell WHERE dataset_id=%s", (args.dataset_id,)) or 0) report["counts"]["features"] = int(scalar(cur, "SELECT COUNT(*) FROM chart_feature WHERE dataset_id=%s", (args.dataset_id,)) or 0) report["counts"]["attributes"] = int(scalar(cur, """ SELECT COUNT(*) FROM chart_feature_attribute a JOIN chart_feature f ON f.id=a.feature_id WHERE f.dataset_id=%s """, (args.dataset_id,)) or 0) report["counts"]["sounding_points"] = int(scalar(cur, """ SELECT COUNT(*) FROM chart_sounding_point s JOIN chart_feature f ON f.id=s.feature_id WHERE f.dataset_id=%s """, (args.dataset_id,)) or 0) report["counts"]["raw_geotype16_meta"] = int(scalar(cur, """ SELECT COUNT(*) FROM chart_feature WHERE dataset_id=%s AND raw_geotype=16 AND geometry_type='META' """, (args.dataset_id,)) or 0) report["counts"]["raw_geotype16_geometry_issue"] = int(scalar(cur, """ SELECT COUNT(*) FROM chart_feature WHERE dataset_id=%s AND raw_geotype=16 AND geometry_type='META' AND geometry_issue IS NOT NULL """, (args.dataset_id,)) or 0) print("[2/8] 检查 Feature geometry_type ↔ 实际 geometry 表映射...") mapping_issues, mapping_stats = audit_mapping(cur, args.dataset_id) report["mapping"] = mapping_stats write_csv(outdir / "feature_geometry_mapping_issues.csv", mapping_issues) print(f" mapping issues: {len(mapping_issues):,}") all_intrinsic_issues: List[Dict[str, Any]] = [] geometry_stats: List[Dict[str, Any]] = [] print("[3/8] 扫描 Point geometry...") st, iss = scan_feature_geometry_table( cur, table="chart_point_geometry", kind="POINT", dataset_id=args.dataset_id, wkb_fn=funcs["wkb_fn"], srid_fn=funcs["srid_fn"], dataset_bbox=dataset_bbox, chunk_size=args.chunk_size, bbox_tolerance=args.bbox_tolerance, boundary_tolerance=args.boundary_tolerance, ) geometry_stats.append(st) all_intrinsic_issues.extend(iss) print("[4/8] 扫描 Line geometry...") st, iss = scan_feature_geometry_table( cur, table="chart_line_geometry", kind="LINE", dataset_id=args.dataset_id, wkb_fn=funcs["wkb_fn"], srid_fn=funcs["srid_fn"], dataset_bbox=dataset_bbox, chunk_size=args.chunk_size, bbox_tolerance=args.bbox_tolerance, boundary_tolerance=args.boundary_tolerance, ) geometry_stats.append(st) all_intrinsic_issues.extend(iss) print("[5/8] 扫描 Area geometry...") st, iss = scan_feature_geometry_table( cur, table="chart_area_geometry", kind="AREA", dataset_id=args.dataset_id, wkb_fn=funcs["wkb_fn"], srid_fn=funcs["srid_fn"], dataset_bbox=dataset_bbox, chunk_size=args.chunk_size, bbox_tolerance=args.bbox_tolerance, boundary_tolerance=args.boundary_tolerance, ) geometry_stats.append(st) all_intrinsic_issues.extend(iss) print("[6/8] 扫描 Sounding points，并按 Feature 对账 bbox...") st, iss, sounding_bbox_issues = scan_soundings( cur, dataset_id=args.dataset_id, wkb_fn=funcs["wkb_fn"], srid_fn=funcs["srid_fn"], dataset_bbox=dataset_bbox, chunk_size=max(args.chunk_size, 5000), bbox_tolerance=args.bbox_tolerance, boundary_tolerance=args.boundary_tolerance, ) geometry_stats.append(st) all_intrinsic_issues.extend(iss) write_csv(outdir / "sounding_feature_bbox_issues.csv", sounding_bbox_issues) coverage_issues: List[Dict[str, Any]] = [] if not args.skip_coverage: print("[7/8] 扫描 Coverage geometry...") st, coverage_issues = scan_coverage( cur, dataset_id=args.dataset_id, wkb_fn=funcs["wkb_fn"], srid_fn=funcs["srid_fn"], dataset_bbox=dataset_bbox, bbox_tolerance=args.bbox_tolerance, boundary_tolerance=args.boundary_tolerance, ) geometry_stats.append(st) write_csv(outdir / "coverage_geometry_issues.csv", coverage_issues) else: print("[7/8] Coverage skipped by option.") print("[8/8] 汇总输出...") report["geometry_stats"] = geometry_stats # 主 intrinsic issue 输出（Feature Point/Line/Area + sounding point） write_csv(outdir / "intrinsic_geometry_issues.csv", all_intrinsic_issues) write_csv(outdir / "geometry_stats_by_kind.csv", geometry_stats) cat_counts = Counter( x.get("issue_category", "UNKNOWN") for x in all_intrinsic_issues ) cat_rows = [ {"issue_category": k, "count": v} for k, v in cat_counts.most_common() ] write_csv( outdir / "intrinsic_issue_category_counts.csv", cat_rows, ["issue_category", "count"] ) invalid_reason_counts = Counter( x.get("invalid_reason") for x in all_intrinsic_issues if x.get("invalid_reason") ) write_csv( outdir / "invalid_geometry_reason_counts.csv", [ {"reason": k, "count": v} for k, v in invalid_reason_counts.most_common() ], ["reason", "count"] ) obj_counts = Counter( (x.get("object_class") or "", x.get("issue_category") or "") for x in all_intrinsic_issues ) write_csv( outdir / "issue_object_class_counts.csv", [ {"object_class": k[0], "issue_category": k[1], "count": v} for k, v in obj_counts.most_common() ], ["object_class", "issue_category", "count"] ) cell_counts = Counter( (x.get("cell_name") or "", x.get("issue_category") or "") for x in all_intrinsic_issues ) write_csv( outdir / "issue_cell_counts.csv", [ {"cell_name": k[0], "issue_category": k[1], "count": v} for k, v in cell_counts.most_common() ], ["cell_name", "issue_category", "count"] ) # 每类最多 10 条样本 samples = [] sample_seen = Counter() for x in all_intrinsic_issues: cat = x.get("issue_category", "UNKNOWN") if sample_seen[cat] < 10: samples.append(x) sample_seen[cat] += 1 write_csv(outdir / "issue_samples.csv", samples) report["issue_summary"] = { "mapping_issue_count": len(mapping_issues), "intrinsic_issue_count": len(all_intrinsic_issues), "sounding_feature_bbox_issue_count": len(sounding_bbox_issues), "coverage_issue_count": len(coverage_issues), "intrinsic_by_category": dict(cat_counts), "issues_with_cell_boundary_context": sum( 1 for x in all_intrinsic_issues if "CELL_EDGE" in (x.get("boundary_context") or "") ), "issues_with_dataset_bbox_context": sum( 1 for x in all_intrinsic_issues if "DATASET_BBOX_EDGE" in (x.get("boundary_context") or "") ), } # 基线匹配 acceptance = {} for k, exp in EXPECTED_DATASET1.items(): actual = report["counts"].get(k) acceptance[k] = { "actual": actual, "expected": exp, "pass": actual == exp, } acceptance["raw_geotype16_expected_parent_count"] = { "actual": report["counts"]["raw_geotype16_meta"], "expected": 4406, "pass": report["counts"]["raw_geotype16_meta"] == 4406, "interpretation": "EXPECTED_PARENT_NO_GEOMETRY, not intrinsic geometry corruption", } acceptance["mapping"] = { "issue_count": len(mapping_issues), "pass": len(mapping_issues) == 0, } # 注意：intrinsic issue 非 0 时不自动称“Importer FAIL”。 # 需要结合 category、对象类型和 boundary_context 分析。 acceptance["intrinsic_geometry"] = { "issue_count": len(all_intrinsic_issues), "pass_if_zero": len(all_intrinsic_issues) == 0, "note": ( "非零时需要人工分类；局部 Dataset 边界不连续本身不属于 FAIL。" ), } report["acceptance"] = acceptance # JSON (outdir / "spatial_audit_report.json").write_text( json.dumps(report, ensure_ascii=False, indent=2, default=str), encoding="utf-8" ) # TXT summary lines = [] lines.append("CM93 -> MySQL 实体空间 Geometry 只读审计 v3") lines.append("=" * 78) lines.append(f"Generated: {report['generated_at']}") lines.append(f"Database: {args.database}") lines.append(f"Dataset ID: {args.dataset_id}") lines.append(f"DB version: {report['server']['version']}") lines.append(f"WKB fn: {funcs['wkb_fn']}") lines.append(f"SRID fn: {funcs['srid_fn']}") lines.append(f"Dataset bbox: {dataset_bbox}") lines.append("") lines.append("DATASET SCOPE") lines.append( " 当前 CM93JAPAN 是日本沿海局部子集，不是完整 CM93。" ) lines.append( " 跨未导入海域/邻区的几何或对象连续性不作为本轮 FAIL 条件。" ) lines.append( " boundary_context 仅用于后续解释，不会自动把记录判成 importer bug。" ) lines.append("") lines.append("BASELINE COUNTS") for k in ["cells", "features", "attributes", "sounding_points", "raw_geotype16_meta", "raw_geotype16_geometry_issue"]: lines.append(f" {k:36s} {report['counts'].get(k)}") lines.append("") lines.append("FEATURE ↔ GEOMETRY TABLE MAPPING") lines.append(f" mapping issues: {len(mapping_issues)}") for k, v in sorted(mapping_stats.get("mapping_issue_categories", {}).items()): lines.append(f" {k:40s} {v}") lines.append("") lines.append("INTRINSIC GEOMETRY STATS") for s in geometry_stats: lines.append( f" {s['kind']:<9} rows={s.get('rows',0):>9} " f"parsed={s.get('parsed',0):>9} " f"empty={s.get('empty',0):>6} " f"valid={s.get('valid',0):>9} " f"invalid={s.get('invalid',0):>6} " f"issues={s.get('issue_count', s.get('issue_count_excluding_feature_bbox',0)):>6}" ) lines.append( f" actual types={s.get('actual_geom_types')}" ) if s["kind"] == "LINE": lines.append( f" non_simple_line_info={s.get('non_simple_line_info',0)} " f"(informational; not automatically invalid)" ) lines.append("") lines.append("INTRINSIC ISSUE CATEGORIES") if cat_counts: for k, v in cat_counts.most_common(): lines.append(f" {k:36s} {v}") else: lines.append(" none") lines.append("") lines.append("SOUNDING FEATURE BBOX") lines.append(f" bbox issues: {len(sounding_bbox_issues)}") lines.append("") lines.append("COVERAGE") lines.append(f" coverage issues: {len(coverage_issues)}") lines.append("") lines.append("BOUNDARY CONTEXT (INFORMATIONAL)") lines.append( f" intrinsic issues touching cell bbox edge: " f"{report['issue_summary']['issues_with_cell_boundary_context']}" ) lines.append( f" intrinsic issues touching dataset bbox edge: " f"{report['issue_summary']['issues_with_dataset_bbox_context']}" ) lines.append("") lines.append("IMPORTANT INTERPRETATION") lines.append( " raw_geotype=16 / META 的 4406 条正常 parent-no-geometry " "不进入本轮 intrinsic geometry issue。" ) lines.append( " Line 在 Cell/Dataset 边缘结束、缺少未导入邻区对象，不属于本轮 geometry corruption。" ) lines.append( " Polygon 自相交、WKB 无法解析、坐标越界、类型串表等属于单 geometry 内在问题，" "需要单独分析。" ) lines.append("") lines.append("PBF: OUT OF SCOPE / NOT STARTED") (outdir / "spatial_audit_summary.txt").write_text( "\n".join(lines), encoding="utf-8" ) conn.rollback() zpath = zip_directory(outdir) print("") print("审计完成。") print(f"输出目录: {outdir.resolve()}") print(f"ZIP: {zpath.resolve()}") print("") 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()