978 lines
35 KiB
Python
978 lines
35 KiB
Python
#!/usr/bin/env python3
|
||
# -*- coding: utf-8 -*-
|
||
"""
|
||
CM93 -> MySQL 无效 Area / Coverage 深度分类审计
|
||
Version: v4.0
|
||
|
||
背景：
|
||
- Dataset ID=1 为日本沿海局部 CM93 数据子集，不是完整 CM93。
|
||
- V3 已确认：
|
||
* Point 113346：全部有效
|
||
* Line 121773：全部有效
|
||
* Sounding 275700：全部有效
|
||
* Area 81080：52 个 Shapely invalid
|
||
* Coverage 1156：3 个 Shapely invalid
|
||
* feature ↔ geometry table mapping issue = 0
|
||
* raw_geotype=16 / META 的 4406 条为 EXPECTED_PARENT_NO_GEOMETRY
|
||
|
||
本脚本目标：
|
||
1. 重新找出 invalid Area / Coverage；
|
||
2. 对每一个 Shapely 报告的 invalidity point 做局部几何分析；
|
||
3. 区分：
|
||
- invalid 点本身位于 Cell 边界；
|
||
- geometry 只是“整体碰 Cell 边界”，但 invalid 点在内部；
|
||
- 重复顶点自接触；
|
||
- 线段内部真正交叉；
|
||
- shell / hole / ring 结构问题；
|
||
- Interior disconnected；
|
||
4. 找出多个 Feature 共享同一 invalid point 的簇；
|
||
5. 找出 Area 与 Coverage 在同一 Cell、同一 invalid point 的关联；
|
||
6. 输出当前数据库中 topology/edge/node/primitive 等相关 schema 线索，供下一阶段追 importer。
|
||
7. 全程只读，不修 geometry，不改 MySQL，不做 PBF。
|
||
|
||
重要：
|
||
“当前 Dataset 不是完整 CM93”只影响跨未导入海域的连续性判断。
|
||
一个已经生成的 Polygon 如果在自身内部发生真正 segment crossing，
|
||
仍属于需要调查的内在 geometry 问题。
|
||
但如果 invalidity point 正好位于 Cell 裁切边界、并被多个对象/coverage 共用，
|
||
更可能是 Cell clipping / shared topology / ring assembly 类问题，不能直接归为原始数据损坏。
|
||
|
||
依赖：
|
||
pip install pymysql shapely
|
||
|
||
运行：
|
||
python cm93_mysql_invalid_area_audit_v4.py ^
|
||
--host 127.0.0.1 ^
|
||
--port 3306 ^
|
||
--user root ^
|
||
--database cm93_chart ^
|
||
--dataset-id 1
|
||
|
||
输出目录会自动压缩成 zip，直接把 zip 发给 ChatGPT。
|
||
"""
|
||
|
||
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:
|
||
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
|
||
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 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 write_csv(path: Path, rows: Iterable[Dict[str, Any]], fieldnames: Optional[List[str]] = None):
|
||
rows = list(rows)
|
||
if fieldnames is None:
|
||
seen = set()
|
||
fieldnames = []
|
||
for row in rows:
|
||
for k in row.keys():
|
||
if k not in seen:
|
||
seen.add(k)
|
||
fieldnames.append(k)
|
||
|
||
path.parent.mkdir(parents=True, exist_ok=True)
|
||
with path.open("w", encoding="utf-8-sig", newline="") as f:
|
||
if not fieldnames:
|
||
return
|
||
w = csv.DictWriter(f, fieldnames=fieldnames, extrasaction="ignore")
|
||
w.writeheader()
|
||
for row in rows:
|
||
out = {}
|
||
for k in fieldnames:
|
||
v = row.get(k)
|
||
if isinstance(v, (tuple, list, dict)):
|
||
v = json.dumps(v, ensure_ascii=False)
|
||
out[k] = "" if v is None else v
|
||
w.writerow(out)
|
||
|
||
|
||
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}")
|
||
try:
|
||
cur.execute("START TRANSACTION")
|
||
except Exception:
|
||
pass
|
||
|
||
|
||
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], prefix: str = "") -> Optional[Tuple[float, float, float, float]]:
|
||
keys = [
|
||
prefix + "bbox_min_lon",
|
||
prefix + "bbox_min_lat",
|
||
prefix + "bbox_max_lon",
|
||
prefix + "bbox_max_lat",
|
||
]
|
||
vals = []
|
||
for k in keys:
|
||
try:
|
||
v = float(row[k])
|
||
if not math.isfinite(v):
|
||
return None
|
||
vals.append(v)
|
||
except Exception:
|
||
return None
|
||
return vals[0], vals[1], vals[2], vals[3]
|
||
|
||
|
||
def point_cell_edge_context(
|
||
p: Optional[Tuple[float, float]],
|
||
bbox: Optional[Tuple[float, float, float, float]],
|
||
tol: float,
|
||
) -> Tuple[bool, str, Optional[float]]:
|
||
if not p or not bbox:
|
||
return False, "", None
|
||
x, y = p
|
||
minx, miny, maxx, maxy = bbox
|
||
d = {
|
||
"WEST": abs(x - minx),
|
||
"EAST": abs(x - maxx),
|
||
"SOUTH": abs(y - miny),
|
||
"NORTH": abs(y - maxy),
|
||
}
|
||
edge, dist = min(d.items(), key=lambda kv: kv[1])
|
||
return dist <= tol, edge if dist <= tol else "", dist
|
||
|
||
|
||
def point_feature_bbox_edge_context(
|
||
p: Optional[Tuple[float, float]],
|
||
bbox: Optional[Tuple[float, float, float, float]],
|
||
tol: float,
|
||
) -> Tuple[bool, str, Optional[float]]:
|
||
return point_cell_edge_context(p, bbox, tol)
|
||
|
||
|
||
def ring_coordinate_lists(geom) -> List[Tuple[str, int, List[Tuple[float, float]]]]:
|
||
out = []
|
||
|
||
def add_polygon(poly, poly_index: int):
|
||
try:
|
||
out.append((f"poly{poly_index}.shell", 0, [(float(x), float(y)) for x, y, *rest in poly.exterior.coords]))
|
||
except Exception:
|
||
pass
|
||
for i, ring in enumerate(poly.interiors):
|
||
try:
|
||
out.append((f"poly{poly_index}.hole", i, [(float(x), float(y)) for x, y, *rest in ring.coords]))
|
||
except Exception:
|
||
pass
|
||
|
||
if geom.geom_type == "Polygon":
|
||
add_polygon(geom, 0)
|
||
elif geom.geom_type == "MultiPolygon":
|
||
for pi, poly in enumerate(geom.geoms):
|
||
add_polygon(poly, pi)
|
||
|
||
return out
|
||
|
||
|
||
def point_close(a: Tuple[float, float], b: Tuple[float, float], tol: float) -> bool:
|
||
return abs(a[0]-b[0]) <= tol and abs(a[1]-b[1]) <= tol
|
||
|
||
|
||
def analyze_invalid_point_on_rings(
|
||
geom,
|
||
invalid_point: Optional[Tuple[float, float]],
|
||
tol: float,
|
||
) -> Dict[str, Any]:
|
||
result = {
|
||
"polygon_count": 0,
|
||
"ring_count": 0,
|
||
"hole_count": 0,
|
||
"total_vertices": 0,
|
||
"invalid_point_vertex_occurrences": 0,
|
||
"invalid_point_distinct_ring_hits": 0,
|
||
"invalid_point_segment_hits": 0,
|
||
"invalid_point_segment_interior_hits": 0,
|
||
"invalid_point_endpoint_hits": 0,
|
||
"repeated_nonconsecutive_vertex_count": 0,
|
||
"zero_length_segment_count": 0,
|
||
"mechanism": "UNRESOLVED",
|
||
"ring_hit_labels": "",
|
||
}
|
||
|
||
if geom.geom_type == "Polygon":
|
||
result["polygon_count"] = 1
|
||
elif geom.geom_type == "MultiPolygon":
|
||
result["polygon_count"] = len(geom.geoms)
|
||
else:
|
||
return result
|
||
|
||
rings = ring_coordinate_lists(geom)
|
||
result["ring_count"] = len(rings)
|
||
result["hole_count"] = sum(1 for label, _, _ in rings if ".hole" in label)
|
||
|
||
ring_hits = set()
|
||
vertex_occurrences = 0
|
||
segment_hits = 0
|
||
segment_interior_hits = 0
|
||
endpoint_hits = 0
|
||
repeated_nonconsecutive = 0
|
||
zero_segments = 0
|
||
|
||
pgeom = Point(invalid_point) if invalid_point else None
|
||
|
||
for label, ring_index, coords in rings:
|
||
result["total_vertices"] += len(coords)
|
||
|
||
# 去掉 ring 最后一个闭合重复点后做“非连续重复顶点”统计。
|
||
core = coords[:-1] if len(coords) > 1 and point_close(coords[0], coords[-1], tol) else coords[:]
|
||
seen_positions: Dict[Tuple[int, int], List[int]] = defaultdict(list)
|
||
scale = 1.0 / max(tol, 1e-12)
|
||
for i, c in enumerate(core):
|
||
key = (round(c[0] * scale), round(c[1] * scale))
|
||
seen_positions[key].append(i)
|
||
for positions in seen_positions.values():
|
||
if len(positions) > 1:
|
||
# 相邻重复点不归到 nonconsecutive；它们通常会形成 zero segment。
|
||
has_nonconsecutive = False
|
||
for i in range(len(positions)):
|
||
for j in range(i+1, len(positions)):
|
||
if abs(positions[j] - positions[i]) > 1:
|
||
has_nonconsecutive = True
|
||
break
|
||
if has_nonconsecutive:
|
||
break
|
||
if has_nonconsecutive:
|
||
repeated_nonconsecutive += 1
|
||
|
||
if invalid_point:
|
||
for c in core:
|
||
if point_close(c, invalid_point, tol):
|
||
vertex_occurrences += 1
|
||
ring_hits.add(f"{label}[{ring_index}]")
|
||
|
||
for i in range(len(coords)-1):
|
||
a = coords[i]
|
||
b = coords[i+1]
|
||
|
||
if point_close(a, b, tol):
|
||
zero_segments += 1
|
||
|
||
if not invalid_point:
|
||
continue
|
||
|
||
try:
|
||
seg = LineString([a, b])
|
||
if seg.distance(pgeom) <= tol:
|
||
segment_hits += 1
|
||
ring_hits.add(f"{label}[{ring_index}]")
|
||
at_a = point_close(a, invalid_point, tol)
|
||
at_b = point_close(b, invalid_point, tol)
|
||
if at_a or at_b:
|
||
endpoint_hits += 1
|
||
else:
|
||
segment_interior_hits += 1
|
||
except Exception:
|
||
pass
|
||
|
||
result["invalid_point_vertex_occurrences"] = vertex_occurrences
|
||
result["invalid_point_distinct_ring_hits"] = len(ring_hits)
|
||
result["invalid_point_segment_hits"] = segment_hits
|
||
result["invalid_point_segment_interior_hits"] = segment_interior_hits
|
||
result["invalid_point_endpoint_hits"] = endpoint_hits
|
||
result["repeated_nonconsecutive_vertex_count"] = repeated_nonconsecutive
|
||
result["zero_length_segment_count"] = zero_segments
|
||
result["ring_hit_labels"] = "|".join(sorted(ring_hits))
|
||
|
||
# 机制分类：只描述几何表现，不直接断言根因。
|
||
if not invalid_point:
|
||
result["mechanism"] = "NO_INVALID_POINT_FROM_VALIDATOR"
|
||
elif segment_interior_hits >= 2 and vertex_occurrences == 0:
|
||
result["mechanism"] = "TRUE_SEGMENT_CROSSING"
|
||
elif vertex_occurrences >= 2:
|
||
result["mechanism"] = "REPEATED_VERTEX_SELF_TOUCH"
|
||
elif vertex_occurrences >= 1 and segment_interior_hits >= 1:
|
||
result["mechanism"] = "VERTEX_TOUCHES_NONADJACENT_SEGMENT"
|
||
elif len(ring_hits) >= 2:
|
||
result["mechanism"] = "MULTI_RING_TOUCH_OR_CROSS"
|
||
elif segment_hits >= 3:
|
||
result["mechanism"] = "MULTI_SEGMENT_TOUCH"
|
||
elif segment_hits >= 2:
|
||
result["mechanism"] = "TWO_SEGMENT_INTERACTION"
|
||
else:
|
||
result["mechanism"] = "UNRESOLVED_LOCAL_TOPOLOGY"
|
||
|
||
return result
|
||
|
||
|
||
def validator_reason_category(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_GEOMETRY"
|
||
|
||
|
||
def scan_invalid_geometries(
|
||
cur,
|
||
*,
|
||
table: str,
|
||
kind: str,
|
||
dataset_id: int,
|
||
wkb_fn: str,
|
||
chunk_size: int,
|
||
point_tol: float,
|
||
edge_tol: float,
|
||
) -> List[Dict[str, Any]]:
|
||
last_id = 0
|
||
invalids = []
|
||
id_col = "feature_id" if kind == "AREA" else "id"
|
||
|
||
while True:
|
||
if kind == "AREA":
|
||
sql = f"""
|
||
SELECT
|
||
g.feature_id AS row_id,
|
||
g.feature_id,
|
||
NULL AS coverage_id,
|
||
{wkb_fn}(g.geom) AS geom_wkb,
|
||
f.object_class,
|
||
f.raw_geotype,
|
||
f.geometry_type,
|
||
f.geometry_issue,
|
||
f.geometry_repaired,
|
||
f.source_record_index,
|
||
f.cell_id,
|
||
c.cell_name,
|
||
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_area_geometry 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
|
||
"""
|
||
else:
|
||
sql = f"""
|
||
SELECT
|
||
v.id AS row_id,
|
||
NULL AS feature_id,
|
||
v.id AS coverage_id,
|
||
{wkb_fn}(v.geom) AS geom_wkb,
|
||
'COVERAGE' AS object_class,
|
||
NULL AS raw_geotype,
|
||
'COVERAGE' AS geometry_type,
|
||
v.geometry_issue,
|
||
v.geometry_repaired,
|
||
v.source_record_index,
|
||
v.cell_id,
|
||
c.cell_name,
|
||
NULL AS bbox_min_lon, NULL AS bbox_min_lat,
|
||
NULL AS bbox_max_lon, NULL AS bbox_max_lat,
|
||
c.bbox_min_lon AS cell_bbox_min_lon,
|
||
c.bbox_min_lat AS cell_bbox_min_lat,
|
||
c.bbox_max_lon AS cell_bbox_max_lon,
|
||
c.bbox_max_lat AS cell_bbox_max_lat
|
||
FROM chart_coverage v
|
||
JOIN chart_cell c ON c.id=v.cell_id
|
||
WHERE v.dataset_id=%s
|
||
AND v.id>%s
|
||
ORDER BY v.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["row_id"])
|
||
wkb = r["geom_wkb"]
|
||
if isinstance(wkb, memoryview):
|
||
wkb = wkb.tobytes()
|
||
elif isinstance(wkb, bytearray):
|
||
wkb = bytes(wkb)
|
||
|
||
try:
|
||
geom = load_wkb(wkb)
|
||
except Exception:
|
||
continue
|
||
|
||
try:
|
||
valid = bool(geom.is_valid)
|
||
except Exception:
|
||
continue
|
||
if valid:
|
||
continue
|
||
|
||
reason = explain_validity(geom)
|
||
p = parse_invalid_point(reason)
|
||
cb = bbox_from_row(r, "cell_")
|
||
fb = bbox_from_row(r)
|
||
|
||
on_cell_edge, cell_edge, cell_edge_dist = point_cell_edge_context(p, cb, edge_tol)
|
||
on_feature_bbox_edge, feature_bbox_edge, feature_bbox_edge_dist = point_feature_bbox_edge_context(
|
||
p, fb, edge_tol
|
||
)
|
||
|
||
local = analyze_invalid_point_on_rings(geom, p, point_tol)
|
||
|
||
geom_bounds = tuple(float(v) for v in geom.bounds)
|
||
|
||
row = {
|
||
"kind": kind,
|
||
"feature_id": r.get("feature_id"),
|
||
"coverage_id": r.get("coverage_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"),
|
||
"raw_geotype": r.get("raw_geotype"),
|
||
"geometry_type": r.get("geometry_type"),
|
||
"db_geometry_issue": r.get("geometry_issue"),
|
||
"db_geometry_repaired": r.get("geometry_repaired"),
|
||
"validator_category": validator_reason_category(reason),
|
||
"validator_reason": reason,
|
||
"invalid_x": p[0] if p else None,
|
||
"invalid_y": p[1] if p else None,
|
||
"cell_bbox": cb,
|
||
"feature_bbox": fb,
|
||
"geometry_bounds": geom_bounds,
|
||
"invalid_point_on_cell_edge": int(on_cell_edge),
|
||
"cell_edge": cell_edge,
|
||
"cell_edge_distance_deg": cell_edge_dist,
|
||
"invalid_point_on_feature_bbox_edge": int(on_feature_bbox_edge),
|
||
"feature_bbox_edge": feature_bbox_edge,
|
||
"feature_bbox_edge_distance_deg": feature_bbox_edge_dist,
|
||
}
|
||
row.update(local)
|
||
invalids.append(row)
|
||
|
||
print(f" {kind:<8} scanned_to={last_id} invalid_found={len(invalids)}")
|
||
|
||
return invalids
|
||
|
||
|
||
def rounded_point_key(row: Dict[str, Any], decimals: int = 7) -> Optional[Tuple[float, float]]:
|
||
try:
|
||
x = float(row["invalid_x"])
|
||
y = float(row["invalid_y"])
|
||
if math.isfinite(x) and math.isfinite(y):
|
||
return round(x, decimals), round(y, decimals)
|
||
except Exception:
|
||
pass
|
||
return None
|
||
|
||
|
||
def build_shared_point_groups(all_rows: List[Dict[str, Any]]) -> List[Dict[str, Any]]:
|
||
groups = defaultdict(list)
|
||
for r in all_rows:
|
||
k = rounded_point_key(r)
|
||
if k is not None:
|
||
groups[k].append(r)
|
||
|
||
out = []
|
||
for p, rs in groups.items():
|
||
if len(rs) < 2:
|
||
continue
|
||
kinds = Counter(r["kind"] for r in rs)
|
||
cells = sorted({str(r.get("cell_name") or "") for r in rs})
|
||
objs = Counter(str(r.get("object_class") or "") for r in rs)
|
||
out.append({
|
||
"invalid_x": p[0],
|
||
"invalid_y": p[1],
|
||
"record_count": len(rs),
|
||
"area_count": kinds.get("AREA", 0),
|
||
"coverage_count": kinds.get("COVERAGE", 0),
|
||
"cells": "|".join(cells),
|
||
"object_classes": "|".join(f"{k}:{v}" for k, v in objs.most_common()),
|
||
"cell_edge_records": sum(int(r.get("invalid_point_on_cell_edge") or 0) for r in rs),
|
||
"feature_ids": "|".join(str(r["feature_id"]) for r in rs if r.get("feature_id") is not None),
|
||
"coverage_ids": "|".join(str(r["coverage_id"]) for r in rs if r.get("coverage_id") is not None),
|
||
})
|
||
|
||
out.sort(key=lambda r: (-r["record_count"], r["invalid_x"], r["invalid_y"]))
|
||
return out
|
||
|
||
|
||
def mark_area_coverage_links(
|
||
areas: List[Dict[str, Any]],
|
||
coverages: List[Dict[str, Any]],
|
||
decimals: int = 7,
|
||
):
|
||
cov_index = defaultdict(list)
|
||
for c in coverages:
|
||
p = rounded_point_key(c, decimals)
|
||
if p is not None:
|
||
cov_index[(str(c.get("cell_name")), p)].append(c)
|
||
|
||
for a in areas:
|
||
p = rounded_point_key(a, decimals)
|
||
matches = cov_index.get((str(a.get("cell_name")), p), []) if p is not None else []
|
||
a["matches_invalid_coverage_same_cell_point"] = int(bool(matches))
|
||
a["matching_coverage_ids"] = "|".join(
|
||
str(x.get("coverage_id")) for x in matches
|
||
)
|
||
|
||
|
||
def schema_topology_inventory(cur, database: str) -> List[Dict[str, Any]]:
|
||
keywords = (
|
||
"edge", "node", "vector", "primitive", "segment",
|
||
"geom", "geometry", "coverage", "ring", "topolog"
|
||
)
|
||
rs = fetch_all(cur, """
|
||
SELECT TABLE_NAME, COLUMN_NAME, DATA_TYPE, COLUMN_TYPE, IS_NULLABLE
|
||
FROM information_schema.COLUMNS
|
||
WHERE TABLE_SCHEMA=%s
|
||
ORDER BY TABLE_NAME, ORDINAL_POSITION
|
||
""", (database,))
|
||
out = []
|
||
for r in rs:
|
||
blob = (r["TABLE_NAME"] + " " + r["COLUMN_NAME"]).lower()
|
||
if any(k in blob for k in keywords):
|
||
out.append(r)
|
||
return out
|
||
|
||
|
||
def table_inventory(cur, database: str) -> List[Dict[str, Any]]:
|
||
return fetch_all(cur, """
|
||
SELECT TABLE_NAME, TABLE_ROWS, ENGINE, TABLE_COMMENT
|
||
FROM information_schema.TABLES
|
||
WHERE TABLE_SCHEMA=%s
|
||
ORDER BY TABLE_NAME
|
||
""", (database,))
|
||
|
||
|
||
def classify_for_next_step(r: Dict[str, Any]) -> str:
|
||
"""
|
||
这是“调查优先级/表现分类”，不是最终根因。
|
||
"""
|
||
vc = r.get("validator_category")
|
||
mechanism = r.get("mechanism")
|
||
on_edge = bool(r.get("invalid_point_on_cell_edge"))
|
||
coverage_link = bool(r.get("matches_invalid_coverage_same_cell_point"))
|
||
|
||
if vc == "INTERIOR_DISCONNECTED":
|
||
return "INTERNAL_RING_STRUCTURE_NEEDS_REVIEW"
|
||
|
||
if on_edge and coverage_link:
|
||
return "CELL_BOUNDARY_SHARED_WITH_COVERAGE"
|
||
if on_edge and mechanism in {
|
||
"REPEATED_VERTEX_SELF_TOUCH",
|
||
"VERTEX_TOUCHES_NONADJACENT_SEGMENT",
|
||
"MULTI_RING_TOUCH_OR_CROSS",
|
||
"MULTI_SEGMENT_TOUCH",
|
||
"TWO_SEGMENT_INTERACTION",
|
||
}:
|
||
return "CELL_BOUNDARY_TOPOLOGY_CANDIDATE"
|
||
if on_edge and mechanism == "TRUE_SEGMENT_CROSSING":
|
||
return "CELL_BOUNDARY_SEGMENT_CROSSING"
|
||
if not on_edge and mechanism == "TRUE_SEGMENT_CROSSING":
|
||
return "INTERNAL_TRUE_SEGMENT_CROSSING"
|
||
if not on_edge and mechanism == "REPEATED_VERTEX_SELF_TOUCH":
|
||
return "INTERNAL_REPEATED_VERTEX_SELF_TOUCH"
|
||
if not on_edge:
|
||
return "INTERNAL_INVALID_GEOMETRY_NEEDS_REVIEW"
|
||
return "BOUNDARY_INVALID_GEOMETRY_NEEDS_REVIEW"
|
||
|
||
|
||
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/Coverage deep audit v4")
|
||
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(
|
||
"--point-tolerance",
|
||
type=float,
|
||
default=1e-8,
|
||
help="invalid point 与 ring 顶点/segment 的比较容差（degree）"
|
||
)
|
||
p.add_argument(
|
||
"--cell-edge-tolerance",
|
||
type=float,
|
||
default=5e-6,
|
||
help="invalidity point 是否位于 Cell bbox 边界的容差（degree）"
|
||
)
|
||
p.add_argument("--output", default=None)
|
||
return p.parse_args()
|
||
|
||
|
||
def main():
|
||
args = parse_args()
|
||
if not args.password:
|
||
args.password = getpass.getpass("MySQL password: ")
|
||
|
||
stamp = datetime.now().strftime("%Y%m%d_%H%M%S")
|
||
outdir = Path(args.output or f"cm93_invalid_area_audit_v4_dataset_{args.dataset_id}_{stamp}")
|
||
outdir.mkdir(parents=True, exist_ok=True)
|
||
|
||
report = {
|
||
"version": "v4.0",
|
||
"generated_at": datetime.now().isoformat(timespec="seconds"),
|
||
"database": args.database,
|
||
"dataset_id": args.dataset_id,
|
||
"read_only": True,
|
||
"pbf": "OUT_OF_SCOPE",
|
||
"dataset_scope_note": (
|
||
"CM93JAPAN 只是局部子集。跨未导入海域连续性不作为 FAIL；"
|
||
"本轮只深挖已经生成的 invalid Polygon/Coverage 自身。"
|
||
),
|
||
"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] 深扫 invalid Area ...")
|
||
areas = scan_invalid_geometries(
|
||
cur,
|
||
table="chart_area_geometry",
|
||
kind="AREA",
|
||
dataset_id=args.dataset_id,
|
||
wkb_fn=wkb_fn,
|
||
chunk_size=args.chunk_size,
|
||
point_tol=args.point_tolerance,
|
||
edge_tol=args.cell_edge_tolerance,
|
||
)
|
||
|
||
print("[2/6] 深扫 invalid Coverage ...")
|
||
coverages = scan_invalid_geometries(
|
||
cur,
|
||
table="chart_coverage",
|
||
kind="COVERAGE",
|
||
dataset_id=args.dataset_id,
|
||
wkb_fn=wkb_fn,
|
||
chunk_size=args.chunk_size,
|
||
point_tol=args.point_tolerance,
|
||
edge_tol=args.cell_edge_tolerance,
|
||
)
|
||
|
||
print("[3/6] 建立 Area ↔ Coverage 同点关联 ...")
|
||
mark_area_coverage_links(areas, coverages)
|
||
|
||
for r in areas:
|
||
r["next_step_class"] = classify_for_next_step(r)
|
||
for r in coverages:
|
||
r["matches_invalid_coverage_same_cell_point"] = 0
|
||
r["matching_coverage_ids"] = ""
|
||
r["next_step_class"] = classify_for_next_step(r)
|
||
|
||
all_invalid = areas + coverages
|
||
shared_groups = build_shared_point_groups(all_invalid)
|
||
|
||
print("[4/6] 导出 schema/topology 线索 ...")
|
||
topology_schema = schema_topology_inventory(cur, args.database)
|
||
tables = table_inventory(cur, args.database)
|
||
|
||
print("[5/6] 汇总分类 ...")
|
||
validator_counts = Counter(r["validator_category"] for r in all_invalid)
|
||
mechanism_counts = Counter(r["mechanism"] for r in all_invalid)
|
||
next_counts = Counter(r["next_step_class"] for r in all_invalid)
|
||
area_next_counts = Counter(r["next_step_class"] for r in areas)
|
||
cov_next_counts = Counter(r["next_step_class"] for r in coverages)
|
||
|
||
summary = {
|
||
"invalid_area_count": len(areas),
|
||
"invalid_coverage_count": len(coverages),
|
||
"invalid_total": len(all_invalid),
|
||
"invalid_point_on_cell_edge_total": sum(int(r["invalid_point_on_cell_edge"]) for r in all_invalid),
|
||
"invalid_point_internal_total": sum(1 for r in all_invalid if not r["invalid_point_on_cell_edge"]),
|
||
"area_invalid_point_on_cell_edge": sum(int(r["invalid_point_on_cell_edge"]) for r in areas),
|
||
"area_invalid_point_internal": sum(1 for r in areas if not r["invalid_point_on_cell_edge"]),
|
||
"coverage_invalid_point_on_cell_edge": sum(int(r["invalid_point_on_cell_edge"]) for r in coverages),
|
||
"area_matches_invalid_coverage_same_cell_point": sum(
|
||
int(r.get("matches_invalid_coverage_same_cell_point") or 0)
|
||
for r in areas
|
||
),
|
||
"validator_categories": dict(validator_counts),
|
||
"mechanisms": dict(mechanism_counts),
|
||
"next_step_classes": dict(next_counts),
|
||
"area_next_step_classes": dict(area_next_counts),
|
||
"coverage_next_step_classes": dict(cov_next_counts),
|
||
"shared_invalid_point_group_count": len(shared_groups),
|
||
}
|
||
report["summary"] = summary
|
||
|
||
write_csv(outdir / "invalid_areas_deep_classified.csv", areas)
|
||
write_csv(outdir / "invalid_coverages_deep_classified.csv", coverages)
|
||
write_csv(outdir / "invalid_all_deep_classified.csv", all_invalid)
|
||
write_csv(outdir / "shared_invalid_point_groups.csv", shared_groups)
|
||
write_csv(outdir / "topology_schema_inventory.csv", topology_schema)
|
||
write_csv(outdir / "database_table_inventory.csv", tables)
|
||
|
||
write_csv(
|
||
outdir / "validator_category_counts.csv",
|
||
[{"category": k, "count": v} for k, v in validator_counts.most_common()],
|
||
["category", "count"]
|
||
)
|
||
write_csv(
|
||
outdir / "local_mechanism_counts.csv",
|
||
[{"mechanism": k, "count": v} for k, v in mechanism_counts.most_common()],
|
||
["mechanism", "count"]
|
||
)
|
||
write_csv(
|
||
outdir / "next_step_class_counts.csv",
|
||
[{"next_step_class": k, "count": v} for k, v in next_counts.most_common()],
|
||
["next_step_class", "count"]
|
||
)
|
||
|
||
# 重点样本：每一 next_step_class 最多 10 条
|
||
samples = []
|
||
sc = Counter()
|
||
for r in sorted(
|
||
all_invalid,
|
||
key=lambda x: (
|
||
x.get("next_step_class", ""),
|
||
str(x.get("cell_name", "")),
|
||
int(x.get("feature_id") or x.get("coverage_id") or 0)
|
||
)
|
||
):
|
||
c = r["next_step_class"]
|
||
if sc[c] < 10:
|
||
samples.append(r)
|
||
sc[c] += 1
|
||
write_csv(outdir / "deep_classification_samples.csv", samples)
|
||
|
||
(outdir / "invalid_area_audit_report_v4.json").write_text(
|
||
json.dumps(report, ensure_ascii=False, indent=2, default=str),
|
||
encoding="utf-8"
|
||
)
|
||
|
||
lines = []
|
||
lines.append("CM93 -> MySQL Invalid Area/Coverage 深度分类审计 v4")
|
||
lines.append("=" * 82)
|
||
lines.append(f"Generated: {report['generated_at']}")
|
||
lines.append(f"Database: {args.database}")
|
||
lines.append(f"Dataset ID: {args.dataset_id}")
|
||
lines.append(f"MySQL: {report['mysql_version']}")
|
||
lines.append("")
|
||
lines.append("DATASET SCOPE")
|
||
lines.append(" 当前 CM93JAPAN 是局部数据，不要求与未导入海域连续。")
|
||
lines.append(" 但本轮检查的是已生成 Polygon 本身的自交/自接触机制。")
|
||
lines.append("")
|
||
lines.append("COUNTS")
|
||
lines.append(f" invalid Area: {len(areas)}")
|
||
lines.append(f" invalid Coverage: {len(coverages)}")
|
||
lines.append(f" total: {len(all_invalid)}")
|
||
lines.append("")
|
||
lines.append("INVALID POINT LOCATION")
|
||
lines.append(
|
||
f" Area invalid point ON cell edge: "
|
||
f"{summary['area_invalid_point_on_cell_edge']}"
|
||
)
|
||
lines.append(
|
||
f" Area invalid point INTERNAL: "
|
||
f"{summary['area_invalid_point_internal']}"
|
||
)
|
||
lines.append(
|
||
f" Coverage invalid point ON cell edge:"
|
||
f" {summary['coverage_invalid_point_on_cell_edge']}"
|
||
)
|
||
lines.append("")
|
||
lines.append("AREA ↔ COVERAGE")
|
||
lines.append(
|
||
f" Area same cell + same invalid point as invalid Coverage: "
|
||
f"{summary['area_matches_invalid_coverage_same_cell_point']}"
|
||
)
|
||
lines.append("")
|
||
lines.append("VALIDATOR CATEGORY")
|
||
for k, v in validator_counts.most_common():
|
||
lines.append(f" {k:42s} {v}")
|
||
lines.append("")
|
||
lines.append("LOCAL MECHANISM")
|
||
for k, v in mechanism_counts.most_common():
|
||
lines.append(f" {k:42s} {v}")
|
||
lines.append("")
|
||
lines.append("NEXT-STEP CLASS (表现分类，不是最终根因)")
|
||
for k, v in next_counts.most_common():
|
||
lines.append(f" {k:42s} {v}")
|
||
lines.append("")
|
||
lines.append("SHARED INVALID POINT GROUPS")
|
||
lines.append(f" groups with >=2 invalid records: {len(shared_groups)}")
|
||
lines.append("")
|
||
lines.append("INTERPRETATION RULES")
|
||
lines.append(
|
||
" CELL_BOUNDARY_SHARED_WITH_COVERAGE / CELL_BOUNDARY_TOPOLOGY_CANDIDATE "
|
||
"优先调查 Cell clipping / shared topology / importer ring assembly。"
|
||
)
|
||
lines.append(
|
||
" INTERNAL_TRUE_SEGMENT_CROSSING 优先调查源拓扑与 importer edge/ring 顺序。"
|
||
)
|
||
lines.append(
|
||
" INTERNAL_REPEATED_VERTEX_SELF_TOUCH 需要判断 CM93 是否允许 touch-ring 表达，"
|
||
"以及 importer 是否应输出 MultiPolygon/拆环。"
|
||
)
|
||
lines.append(
|
||
" 本脚本不执行 make_valid、buffer(0)、闭环、补点、删点或任何数据库写入。"
|
||
)
|
||
lines.append("")
|
||
lines.append("PBF: OUT OF SCOPE / NOT STARTED")
|
||
|
||
(outdir / "invalid_area_audit_summary_v4.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()
|