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