1042 lines
37 KiB
Python
1042 lines
37 KiB
Python
#!/usr/bin/env python3
|
||
# -*- coding: utf-8 -*-
|
||
"""
|
||
CM93 -> MySQL invalid Area/Coverage 源记录 + 邻接拓扑追踪审计
|
||
Version: v5.0
|
||
|
||
用途：
|
||
- 延续 V4，只针对已确认的 52 个 invalid Area + 3 个 invalid Coverage。
|
||
- 不修复 geometry，不写数据库，不做 PBF。
|
||
- 当前 CM93JAPAN 是局部子集；跨未导入海域连续性不作为 FAIL。
|
||
|
||
V5 重点回答：
|
||
1. invalid point 是否被同 Cell 的其他 Area（包括有效 Area）共同使用？
|
||
2. invalid point 是邻接对象的 vertex、segment interior，还是仅落在内部？
|
||
3. 同一点是否形成多个 object class 的共享拓扑簇？
|
||
4. chart_feature.coverage_id 与 invalid Coverage 是否有关联？
|
||
5. source_record_index 能否在 chart_raw_record 中找到对应/邻近原始记录？
|
||
6. invalid Feature 的 attributes 是什么？
|
||
7. 当前 DB 是否保留足够的 raw/topology 中间数据，可以继续从 MySQL 追源？
|
||
8. 生成 GeoJSON，供人工用 QGIS/geojson.io 等查看 invalid polygon、invalid point、共享邻居。
|
||
|
||
依赖：
|
||
pip install pymysql shapely
|
||
|
||
运行：
|
||
python cm93_mysql_invalid_trace_v5.py ^
|
||
--host 127.0.0.1 ^
|
||
--port 3306 ^
|
||
--user root ^
|
||
--database cm93_chart ^
|
||
--dataset-id 1
|
||
|
||
输出自动打 zip，直接上传整个 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:
|
||
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, box, mapping
|
||
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, default=str)
|
||
elif isinstance(v, (bytes, bytearray, memoryview)):
|
||
b = bytes(v)
|
||
v = b.hex()
|
||
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}")
|
||
cur.execute("START TRANSACTION")
|
||
|
||
|
||
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 to_bytes(v: Any) -> Optional[bytes]:
|
||
if v is None:
|
||
return None
|
||
if isinstance(v, memoryview):
|
||
return v.tobytes()
|
||
if isinstance(v, bytearray):
|
||
return bytes(v)
|
||
if isinstance(v, bytes):
|
||
return v
|
||
return bytes(v)
|
||
|
||
|
||
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 tuple(vals) # type: ignore
|
||
|
||
|
||
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 ring_lists(geom) -> List[Tuple[str, List[Tuple[float, float]]]]:
|
||
out = []
|
||
|
||
def add_poly(poly, pi: int):
|
||
try:
|
||
out.append((f"p{pi}.shell", [(float(c[0]), float(c[1])) for c in poly.exterior.coords]))
|
||
except Exception:
|
||
pass
|
||
for hi, ring in enumerate(poly.interiors):
|
||
try:
|
||
out.append((f"p{pi}.hole{hi}", [(float(c[0]), float(c[1])) for c in ring.coords]))
|
||
except Exception:
|
||
pass
|
||
|
||
if geom.geom_type == "Polygon":
|
||
add_poly(geom, 0)
|
||
elif geom.geom_type == "MultiPolygon":
|
||
for pi, poly in enumerate(geom.geoms):
|
||
add_poly(poly, pi)
|
||
return out
|
||
|
||
|
||
def point_relation_to_area_boundary(
|
||
geom,
|
||
pxy: Tuple[float, float],
|
||
tol: float
|
||
) -> Dict[str, Any]:
|
||
"""
|
||
判断一个目标点与某 Area boundary/ring 的关系。
|
||
"""
|
||
p = Point(pxy)
|
||
vertex_occ = 0
|
||
segment_hits = 0
|
||
segment_interior_hits = 0
|
||
endpoint_hits = 0
|
||
ring_labels = set()
|
||
|
||
for label, coords in ring_lists(geom):
|
||
core = coords[:-1] if len(coords) > 1 and point_close(coords[0], coords[-1], tol) else coords
|
||
for c in core:
|
||
if point_close(c, pxy, tol):
|
||
vertex_occ += 1
|
||
ring_labels.add(label)
|
||
|
||
for i in range(len(coords) - 1):
|
||
a, b = coords[i], coords[i+1]
|
||
try:
|
||
seg = LineString([a, b])
|
||
if seg.distance(p) <= tol:
|
||
segment_hits += 1
|
||
ring_labels.add(label)
|
||
at_a = point_close(a, pxy, tol)
|
||
at_b = point_close(b, pxy, tol)
|
||
if at_a or at_b:
|
||
endpoint_hits += 1
|
||
else:
|
||
segment_interior_hits += 1
|
||
except Exception:
|
||
pass
|
||
|
||
try:
|
||
boundary_distance = float(geom.boundary.distance(p))
|
||
except Exception:
|
||
boundary_distance = None
|
||
|
||
try:
|
||
contains = bool(geom.contains(p))
|
||
except Exception:
|
||
contains = False
|
||
try:
|
||
covers = bool(geom.covers(p))
|
||
except Exception:
|
||
covers = False
|
||
|
||
if vertex_occ > 0:
|
||
relation = "VERTEX"
|
||
elif segment_interior_hits > 0:
|
||
relation = "SEGMENT_INTERIOR"
|
||
elif segment_hits > 0:
|
||
relation = "BOUNDARY"
|
||
elif boundary_distance is not None and boundary_distance <= tol:
|
||
relation = "BOUNDARY_NEAR"
|
||
elif contains:
|
||
relation = "INTERIOR"
|
||
elif covers:
|
||
relation = "COVERS_OTHER"
|
||
else:
|
||
relation = "NONE"
|
||
|
||
return {
|
||
"relation": relation,
|
||
"vertex_occurrences": vertex_occ,
|
||
"segment_hits": segment_hits,
|
||
"segment_interior_hits": segment_interior_hits,
|
||
"endpoint_hits": endpoint_hits,
|
||
"ring_labels": "|".join(sorted(ring_labels)),
|
||
"boundary_distance_deg": boundary_distance,
|
||
"contains_point": int(contains),
|
||
"covers_point": int(covers),
|
||
}
|
||
|
||
|
||
def get_columns(cur, db: str, table: str) -> List[Dict[str, Any]]:
|
||
return fetch_all(cur, """
|
||
SELECT COLUMN_NAME, DATA_TYPE, COLUMN_TYPE, IS_NULLABLE, COLUMN_KEY
|
||
FROM information_schema.COLUMNS
|
||
WHERE TABLE_SCHEMA=%s AND TABLE_NAME=%s
|
||
ORDER BY ORDINAL_POSITION
|
||
""", (db, table))
|
||
|
||
|
||
def col_names(cols: List[Dict[str, Any]]) -> List[str]:
|
||
return [str(x["COLUMN_NAME"]) for x in cols]
|
||
|
||
|
||
def find_col(cols: Sequence[str], candidates: Sequence[str]) -> Optional[str]:
|
||
lower = {c.lower(): c for c in cols}
|
||
for cand in candidates:
|
||
if cand.lower() in lower:
|
||
return lower[cand.lower()]
|
||
for cand in candidates:
|
||
for c in cols:
|
||
if cand.lower() in c.lower():
|
||
return c
|
||
return None
|
||
|
||
|
||
def detect_raw_record_link(cur, db: str) -> Dict[str, Any]:
|
||
cols_meta = get_columns(cur, db, "chart_raw_record")
|
||
cols = col_names(cols_meta)
|
||
return {
|
||
"columns": cols_meta,
|
||
"cell_col": find_col(cols, ["cell_id", "cellid"]),
|
||
"record_index_col": find_col(cols, [
|
||
"source_record_index", "record_index", "record_no",
|
||
"record_number", "source_index", "index"
|
||
]),
|
||
"object_class_col": find_col(cols, ["object_class", "class_name", "object"]),
|
||
"payload_cols": [
|
||
c for c in cols
|
||
if any(x in c.lower() for x in (
|
||
"raw", "payload", "data", "blob", "json",
|
||
"hex", "content", "record"
|
||
))
|
||
],
|
||
}
|
||
|
||
|
||
def load_invalid_areas(cur, dataset_id: int, wkb_fn: str) -> List[Dict[str, Any]]:
|
||
cur.execute(f"""
|
||
SELECT
|
||
g.feature_id,
|
||
{wkb_fn}(g.geom) AS geom_wkb,
|
||
f.*,
|
||
c.cell_name,
|
||
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
|
||
ORDER BY g.feature_id
|
||
""", (dataset_id,))
|
||
out = []
|
||
for r in cur:
|
||
try:
|
||
geom = load_wkb(to_bytes(r["geom_wkb"]))
|
||
if not geom.is_valid:
|
||
reason = explain_validity(geom)
|
||
p = parse_invalid_point(reason)
|
||
x = dict(r)
|
||
x["_geom"] = geom
|
||
x["_reason"] = reason
|
||
x["_invalid_point"] = p
|
||
out.append(x)
|
||
except Exception:
|
||
continue
|
||
return out
|
||
|
||
|
||
def load_invalid_coverages(cur, dataset_id: int, wkb_fn: str) -> List[Dict[str, Any]]:
|
||
cur.execute(f"""
|
||
SELECT
|
||
v.*,
|
||
{wkb_fn}(v.geom) AS geom_wkb,
|
||
c.cell_name,
|
||
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,))
|
||
out = []
|
||
for r in cur:
|
||
try:
|
||
geom = load_wkb(to_bytes(r["geom_wkb"]))
|
||
if not geom.is_valid:
|
||
reason = explain_validity(geom)
|
||
p = parse_invalid_point(reason)
|
||
x = dict(r)
|
||
x["_geom"] = geom
|
||
x["_reason"] = reason
|
||
x["_invalid_point"] = p
|
||
out.append(x)
|
||
except Exception:
|
||
continue
|
||
return out
|
||
|
||
|
||
def load_all_areas_for_cells(
|
||
cur,
|
||
dataset_id: int,
|
||
cell_ids: Sequence[int],
|
||
wkb_fn: str,
|
||
) -> Dict[int, List[Dict[str, Any]]]:
|
||
out = defaultdict(list)
|
||
ids = sorted(set(int(x) for x in cell_ids))
|
||
for pos in range(0, len(ids), 100):
|
||
batch = ids[pos:pos+100]
|
||
ph = ",".join(["%s"] * len(batch))
|
||
sql = f"""
|
||
SELECT
|
||
g.feature_id,
|
||
{wkb_fn}(g.geom) AS geom_wkb,
|
||
f.object_class,
|
||
f.source_record_index,
|
||
f.coverage_id,
|
||
f.geometry_issue,
|
||
f.geometry_repaired,
|
||
f.raw_geotype,
|
||
f.cell_id,
|
||
c.cell_name
|
||
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 f.cell_id IN ({ph})
|
||
ORDER BY f.cell_id, f.id
|
||
"""
|
||
cur.execute(sql, [dataset_id] + batch)
|
||
for r in cur.fetchall():
|
||
try:
|
||
r["_geom"] = load_wkb(to_bytes(r["geom_wkb"]))
|
||
out[int(r["cell_id"])].append(r)
|
||
except Exception:
|
||
pass
|
||
return out
|
||
|
||
|
||
def dump_attributes(cur, feature_ids: Sequence[int], db: str) -> List[Dict[str, Any]]:
|
||
if not feature_ids:
|
||
return []
|
||
|
||
attr_cols = col_names(get_columns(cur, db, "chart_feature_attribute"))
|
||
dict_cols = col_names(get_columns(cur, db, "chart_dictionary_attribute"))
|
||
|
||
feature_col = find_col(attr_cols, ["feature_id"])
|
||
attr_id_col = find_col(attr_cols, ["attribute_id", "attr_id", "dictionary_attribute_id"])
|
||
if not feature_col:
|
||
return []
|
||
|
||
# 为最大兼容性，先原样导出 feature_attribute。
|
||
raw_rows = []
|
||
ids = list(sorted(set(int(x) for x in feature_ids)))
|
||
for pos in range(0, len(ids), 500):
|
||
b = ids[pos:pos+500]
|
||
ph = ",".join(["%s"] * len(b))
|
||
raw_rows.extend(fetch_all(
|
||
cur,
|
||
f"SELECT * FROM chart_feature_attribute WHERE {qi(feature_col)} IN ({ph}) ORDER BY {qi(feature_col)}",
|
||
b
|
||
))
|
||
|
||
# 如果能识别字典主键，则补充 dictionary 字段。
|
||
if attr_id_col:
|
||
dict_pk = find_col(dict_cols, ["id", "attribute_id", "attr_id"])
|
||
if dict_pk:
|
||
dict_ids = sorted({
|
||
int(r[attr_id_col])
|
||
for r in raw_rows
|
||
if r.get(attr_id_col) is not None
|
||
})
|
||
dictionary = {}
|
||
for pos in range(0, len(dict_ids), 500):
|
||
b = dict_ids[pos:pos+500]
|
||
ph = ",".join(["%s"] * len(b))
|
||
for d in fetch_all(
|
||
cur,
|
||
f"SELECT * FROM chart_dictionary_attribute WHERE {qi(dict_pk)} IN ({ph})",
|
||
b
|
||
):
|
||
dictionary[str(d[dict_pk])] = d
|
||
|
||
enriched = []
|
||
for r in raw_rows:
|
||
x = dict(r)
|
||
d = dictionary.get(str(r.get(attr_id_col)))
|
||
if d:
|
||
for k, v in d.items():
|
||
x["dict__" + k] = v
|
||
enriched.append(x)
|
||
return enriched
|
||
|
||
return raw_rows
|
||
|
||
|
||
def dump_raw_records(
|
||
cur,
|
||
db: str,
|
||
invalid_rows: List[Dict[str, Any]],
|
||
raw_link: Dict[str, Any],
|
||
) -> Tuple[List[Dict[str, Any]], Dict[str, Any]]:
|
||
cell_col = raw_link.get("cell_col")
|
||
idx_col = raw_link.get("record_index_col")
|
||
result_meta = {
|
||
"raw_record_total": int(scalar(cur, "SELECT COUNT(*) FROM chart_raw_record") or 0),
|
||
"cell_col": cell_col,
|
||
"record_index_col": idx_col,
|
||
"exact_match_supported": bool(cell_col and idx_col),
|
||
"exact_match_count": 0,
|
||
"nearby_row_count": 0,
|
||
}
|
||
|
||
if not cell_col or not idx_col:
|
||
# schema 能看，但无法安全猜连接逻辑。
|
||
sample = fetch_all(cur, "SELECT * FROM chart_raw_record LIMIT 50")
|
||
result_meta["fallback_sample_count"] = len(sample)
|
||
return sample, result_meta
|
||
|
||
output = []
|
||
seen = set()
|
||
|
||
for r in invalid_rows:
|
||
cell_id = r.get("cell_id")
|
||
src_idx = r.get("source_record_index")
|
||
if cell_id is None or src_idx is None:
|
||
continue
|
||
|
||
# 精确 + 前后两条，仅用于调查；不假定 record_index 必然等价。
|
||
rs = fetch_all(cur, f"""
|
||
SELECT *,
|
||
CASE WHEN {qi(idx_col)}=%s THEN 1 ELSE 0 END AS _exact_source_index
|
||
FROM chart_raw_record
|
||
WHERE {qi(cell_col)}=%s
|
||
AND {qi(idx_col)} BETWEEN %s AND %s
|
||
ORDER BY {qi(idx_col)}
|
||
""", (src_idx, cell_id, int(src_idx)-2, int(src_idx)+2))
|
||
|
||
for x in rs:
|
||
key = json.dumps(x, sort_keys=True, default=str)
|
||
if key in seen:
|
||
continue
|
||
seen.add(key)
|
||
y = dict(x)
|
||
y["_target_feature_id"] = r.get("feature_id")
|
||
y["_target_cell_name"] = r.get("cell_name")
|
||
y["_target_source_record_index"] = src_idx
|
||
output.append(y)
|
||
result_meta["nearby_row_count"] += 1
|
||
if int(x.get("_exact_source_index") or 0):
|
||
result_meta["exact_match_count"] += 1
|
||
|
||
return output, result_meta
|
||
|
||
|
||
def make_geojson(
|
||
invalid_areas: List[Dict[str, Any]],
|
||
invalid_coverages: List[Dict[str, Any]],
|
||
neighbor_rows: List[Dict[str, Any]],
|
||
path: Path,
|
||
):
|
||
features = []
|
||
|
||
def add_geom(geom, props):
|
||
try:
|
||
features.append({
|
||
"type": "Feature",
|
||
"geometry": mapping(geom),
|
||
"properties": props,
|
||
})
|
||
except Exception:
|
||
pass
|
||
|
||
# invalid Area + invalid point
|
||
for r in invalid_areas:
|
||
geom = r["_geom"]
|
||
p = r["_invalid_point"]
|
||
add_geom(geom, {
|
||
"layer": "invalid_area",
|
||
"feature_id": r.get("feature_id"),
|
||
"object_class": r.get("object_class"),
|
||
"cell_name": r.get("cell_name"),
|
||
"source_record_index": r.get("source_record_index"),
|
||
"coverage_id": r.get("coverage_id"),
|
||
"validity": r.get("_reason"),
|
||
})
|
||
if p:
|
||
add_geom(Point(p), {
|
||
"layer": "invalid_point",
|
||
"source": "AREA",
|
||
"feature_id": r.get("feature_id"),
|
||
"cell_name": r.get("cell_name"),
|
||
"validity": r.get("_reason"),
|
||
})
|
||
|
||
# invalid Coverage + invalid point
|
||
for r in invalid_coverages:
|
||
geom = r["_geom"]
|
||
p = r["_invalid_point"]
|
||
add_geom(geom, {
|
||
"layer": "invalid_coverage",
|
||
"coverage_id": r.get("id"),
|
||
"cell_name": r.get("cell_name"),
|
||
"source_record_index": r.get("source_record_index"),
|
||
"validity": r.get("_reason"),
|
||
})
|
||
if p:
|
||
add_geom(Point(p), {
|
||
"layer": "invalid_point",
|
||
"source": "COVERAGE",
|
||
"coverage_id": r.get("id"),
|
||
"cell_name": r.get("cell_name"),
|
||
"validity": r.get("_reason"),
|
||
})
|
||
|
||
# 只加入确实共享 invalid point 的 Area 邻居，避免 GeoJSON 爆大。
|
||
added_neighbors = set()
|
||
for n in neighbor_rows:
|
||
if n.get("relation") not in {"VERTEX", "SEGMENT_INTERIOR", "BOUNDARY", "BOUNDARY_NEAR"}:
|
||
continue
|
||
fid = n.get("neighbor_feature_id")
|
||
geom = n.get("_geom")
|
||
if geom is None:
|
||
continue
|
||
key = (fid, n.get("target_invalid_x"), n.get("target_invalid_y"))
|
||
if key in added_neighbors:
|
||
continue
|
||
added_neighbors.add(key)
|
||
add_geom(geom, {
|
||
"layer": "shared_neighbor_area",
|
||
"feature_id": fid,
|
||
"object_class": n.get("neighbor_object_class"),
|
||
"cell_name": n.get("cell_name"),
|
||
"target_feature_id": n.get("target_feature_id"),
|
||
"relation": n.get("relation"),
|
||
"neighbor_valid": n.get("neighbor_valid"),
|
||
})
|
||
|
||
fc = {"type": "FeatureCollection", "features": features}
|
||
path.write_text(
|
||
json.dumps(fc, 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 geometry source/topology trace v5")
|
||
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("--point-tolerance", type=float, default=1e-8)
|
||
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_trace_v5_dataset_{args.dataset_id}_{stamp}")
|
||
outdir.mkdir(parents=True, exist_ok=True)
|
||
|
||
report = {
|
||
"version": "v5.0",
|
||
"generated_at": datetime.now().isoformat(timespec="seconds"),
|
||
"database": args.database,
|
||
"dataset_id": args.dataset_id,
|
||
"read_only": True,
|
||
"pbf": "OUT_OF_SCOPE",
|
||
"scope_note": (
|
||
"CM93JAPAN 为局部数据集。V5 不要求跨未导入海域连续，"
|
||
"只追 invalid Area/Coverage 的源记录和同 Cell 邻接拓扑。"
|
||
),
|
||
"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/7] 找出 invalid Area / Coverage ...")
|
||
invalid_areas = load_invalid_areas(cur, args.dataset_id, wkb_fn)
|
||
invalid_cov = load_invalid_coverages(cur, args.dataset_id, wkb_fn)
|
||
print(f" Area={len(invalid_areas)} Coverage={len(invalid_cov)}")
|
||
|
||
affected_cells = sorted({
|
||
int(r["cell_id"])
|
||
for r in invalid_areas + invalid_cov
|
||
if r.get("cell_id") is not None
|
||
})
|
||
report["affected_cell_count"] = len(affected_cells)
|
||
|
||
print("[2/7] 加载受影响 Cell 的所有 Area，检查共享 invalid point ...")
|
||
all_areas_by_cell = load_all_areas_for_cells(
|
||
cur, args.dataset_id, affected_cells, wkb_fn
|
||
)
|
||
|
||
neighbor_rows = []
|
||
shared_counter = Counter()
|
||
target_summary = []
|
||
|
||
for target in invalid_areas:
|
||
p = target["_invalid_point"]
|
||
if not p:
|
||
continue
|
||
fid = int(target["feature_id"])
|
||
cell_id = int(target["cell_id"])
|
||
neighbors_at_point = 0
|
||
valid_neighbors_at_point = 0
|
||
invalid_neighbors_at_point = 0
|
||
relation_counts = Counter()
|
||
|
||
for n in all_areas_by_cell.get(cell_id, []):
|
||
nf = int(n["feature_id"])
|
||
if nf == fid:
|
||
continue
|
||
geom = n["_geom"]
|
||
|
||
rel = point_relation_to_area_boundary(
|
||
geom, p, args.point_tolerance
|
||
)
|
||
relation = rel["relation"]
|
||
|
||
if relation in {"VERTEX", "SEGMENT_INTERIOR", "BOUNDARY", "BOUNDARY_NEAR"}:
|
||
neighbors_at_point += 1
|
||
try:
|
||
nvalid = bool(geom.is_valid)
|
||
except Exception:
|
||
nvalid = False
|
||
if nvalid:
|
||
valid_neighbors_at_point += 1
|
||
else:
|
||
invalid_neighbors_at_point += 1
|
||
relation_counts[relation] += 1
|
||
|
||
row = {
|
||
"target_feature_id": fid,
|
||
"target_object_class": target.get("object_class"),
|
||
"target_source_record_index": target.get("source_record_index"),
|
||
"target_invalid_reason": target.get("_reason"),
|
||
"target_invalid_x": p[0],
|
||
"target_invalid_y": p[1],
|
||
"cell_id": cell_id,
|
||
"cell_name": target.get("cell_name"),
|
||
"neighbor_feature_id": nf,
|
||
"neighbor_object_class": n.get("object_class"),
|
||
"neighbor_source_record_index": n.get("source_record_index"),
|
||
"neighbor_coverage_id": n.get("coverage_id"),
|
||
"neighbor_valid": int(nvalid),
|
||
"neighbor_validity_reason": "" if nvalid else explain_validity(geom),
|
||
}
|
||
row.update(rel)
|
||
row["_geom"] = geom
|
||
neighbor_rows.append(row)
|
||
|
||
shared_counter[neighbors_at_point] += 1
|
||
target_summary.append({
|
||
"feature_id": fid,
|
||
"object_class": target.get("object_class"),
|
||
"cell_id": cell_id,
|
||
"cell_name": target.get("cell_name"),
|
||
"source_record_index": target.get("source_record_index"),
|
||
"coverage_id": target.get("coverage_id"),
|
||
"invalid_reason": target.get("_reason"),
|
||
"invalid_x": p[0],
|
||
"invalid_y": p[1],
|
||
"shared_neighbor_count": neighbors_at_point,
|
||
"shared_valid_neighbor_count": valid_neighbors_at_point,
|
||
"shared_invalid_neighbor_count": invalid_neighbors_at_point,
|
||
"neighbor_relations": dict(relation_counts),
|
||
})
|
||
|
||
print(f" shared neighbor rows={len(neighbor_rows)}")
|
||
|
||
print("[3/7] 导出 invalid Feature attributes ...")
|
||
feature_ids = [int(r["feature_id"]) for r in invalid_areas]
|
||
attrs = dump_attributes(cur, feature_ids, args.database)
|
||
write_csv(outdir / "invalid_feature_attributes.csv", attrs)
|
||
|
||
print("[4/7] 追 chart_raw_record / source_record_index ...")
|
||
raw_link = detect_raw_record_link(cur, args.database)
|
||
raw_rows, raw_meta = dump_raw_records(
|
||
cur, args.database, invalid_areas, raw_link
|
||
)
|
||
report["raw_record_link"] = raw_meta
|
||
write_csv(outdir / "chart_raw_record_schema.csv", raw_link["columns"])
|
||
write_csv(outdir / "invalid_raw_record_trace.csv", raw_rows)
|
||
|
||
print("[5/7] 导出 Feature / Coverage 完整元数据 ...")
|
||
feature_export = []
|
||
for r in invalid_areas:
|
||
x = {
|
||
k: v for k, v in r.items()
|
||
if not k.startswith("_") and k != "geom_wkb"
|
||
}
|
||
x["validity_reason"] = r["_reason"]
|
||
p = r["_invalid_point"]
|
||
x["invalid_x"] = p[0] if p else None
|
||
x["invalid_y"] = p[1] if p else None
|
||
feature_export.append(x)
|
||
|
||
coverage_export = []
|
||
for r in invalid_cov:
|
||
x = {
|
||
k: v for k, v in r.items()
|
||
if not k.startswith("_") and k not in {"geom_wkb", "geom"}
|
||
}
|
||
x["validity_reason"] = r["_reason"]
|
||
p = r["_invalid_point"]
|
||
x["invalid_x"] = p[0] if p else None
|
||
x["invalid_y"] = p[1] if p else None
|
||
coverage_export.append(x)
|
||
|
||
write_csv(outdir / "invalid_feature_full_metadata.csv", feature_export)
|
||
write_csv(outdir / "invalid_coverage_full_metadata.csv", coverage_export)
|
||
write_csv(outdir / "invalid_target_neighbor_summary.csv", target_summary)
|
||
|
||
# CSV 不应塞 Shapely object
|
||
neighbor_csv = [
|
||
{k: v for k, v in r.items() if k != "_geom"}
|
||
for r in neighbor_rows
|
||
]
|
||
write_csv(outdir / "invalid_point_shared_area_neighbors.csv", neighbor_csv)
|
||
|
||
print("[6/7] 生成 GeoJSON ...")
|
||
make_geojson(
|
||
invalid_areas, invalid_cov, neighbor_rows,
|
||
outdir / "invalid_geometry_trace.geojson"
|
||
)
|
||
|
||
# 汇总
|
||
target_counts = Counter()
|
||
valid_share = 0
|
||
invalid_share = 0
|
||
no_share = 0
|
||
for t in target_summary:
|
||
n = int(t["shared_neighbor_count"])
|
||
if n == 0:
|
||
no_share += 1
|
||
target_counts["NO_SHARED_AREA_BOUNDARY_AT_INVALID_POINT"] += 1
|
||
else:
|
||
target_counts["HAS_SHARED_AREA_BOUNDARY_AT_INVALID_POINT"] += 1
|
||
if int(t["shared_valid_neighbor_count"]) > 0:
|
||
valid_share += 1
|
||
if int(t["shared_invalid_neighbor_count"]) > 0:
|
||
invalid_share += 1
|
||
|
||
# coverage_id 分布
|
||
feature_cov = Counter(
|
||
"HAS_COVERAGE_ID" if r.get("coverage_id") is not None else "NO_COVERAGE_ID"
|
||
for r in invalid_areas
|
||
)
|
||
|
||
report["summary"] = {
|
||
"invalid_area_count": len(invalid_areas),
|
||
"invalid_coverage_count": len(invalid_cov),
|
||
"affected_cell_count": len(affected_cells),
|
||
"invalid_area_with_shared_area_boundary_at_invalid_point":
|
||
target_counts["HAS_SHARED_AREA_BOUNDARY_AT_INVALID_POINT"],
|
||
"invalid_area_without_shared_area_boundary_at_invalid_point":
|
||
no_share,
|
||
"invalid_area_sharing_point_with_at_least_one_valid_area":
|
||
valid_share,
|
||
"invalid_area_sharing_point_with_at_least_one_invalid_area":
|
||
invalid_share,
|
||
"shared_neighbor_relation_counts": dict(Counter(
|
||
r["relation"] for r in neighbor_rows
|
||
)),
|
||
"invalid_feature_coverage_id_presence": dict(feature_cov),
|
||
"raw_record_trace": raw_meta,
|
||
}
|
||
|
||
print("[7/7] 写报告并压缩 ...")
|
||
(outdir / "invalid_trace_report_v5.json").write_text(
|
||
json.dumps(report, ensure_ascii=False, indent=2, default=str),
|
||
encoding="utf-8"
|
||
)
|
||
|
||
lines = []
|
||
lines.append("CM93 -> MySQL invalid Area/Coverage 源记录 + 邻接拓扑追踪 v5")
|
||
lines.append("=" * 88)
|
||
lines.append(f"Generated: {report['generated_at']}")
|
||
lines.append(f"Dataset ID: {args.dataset_id}")
|
||
lines.append(f"MySQL: {report['mysql_version']}")
|
||
lines.append("")
|
||
lines.append("SCOPE")
|
||
lines.append(" CM93JAPAN 为局部数据集；跨未导入海域连续性不作为 FAIL。")
|
||
lines.append(" 本轮不修 geometry，只调查 invalid point 的共享拓扑与 source_record。")
|
||
lines.append("")
|
||
lines.append("INVALID BASELINE")
|
||
lines.append(f" invalid Area: {len(invalid_areas)}")
|
||
lines.append(f" invalid Coverage: {len(invalid_cov)}")
|
||
lines.append(f" affected Cells: {len(affected_cells)}")
|
||
lines.append("")
|
||
lines.append("SHARED AREA TOPOLOGY AT INVALID POINT")
|
||
lines.append(
|
||
f" invalid Area with >=1 other Area boundary at same point: "
|
||
f"{report['summary']['invalid_area_with_shared_area_boundary_at_invalid_point']}"
|
||
)
|
||
lines.append(
|
||
f" invalid Area with no other Area boundary at same point: "
|
||
f"{report['summary']['invalid_area_without_shared_area_boundary_at_invalid_point']}"
|
||
)
|
||
lines.append(
|
||
f" invalid Area sharing point with >=1 VALID Area: "
|
||
f"{report['summary']['invalid_area_sharing_point_with_at_least_one_valid_area']}"
|
||
)
|
||
lines.append(
|
||
f" invalid Area sharing point with >=1 INVALID Area: "
|
||
f"{report['summary']['invalid_area_sharing_point_with_at_least_one_invalid_area']}"
|
||
)
|
||
lines.append("")
|
||
lines.append("NEIGHBOR RELATIONS")
|
||
for k, v in report["summary"]["shared_neighbor_relation_counts"].items():
|
||
lines.append(f" {k:28s} {v}")
|
||
lines.append("")
|
||
lines.append("FEATURE coverage_id")
|
||
for k, v in feature_cov.items():
|
||
lines.append(f" {k:28s} {v}")
|
||
lines.append("")
|
||
lines.append("RAW RECORD TRACE")
|
||
lines.append(f" chart_raw_record rows: {raw_meta['raw_record_total']}")
|
||
lines.append(f" detected cell column: {raw_meta['cell_col']}")
|
||
lines.append(f" detected record index column: {raw_meta['record_index_col']}")
|
||
lines.append(f" exact source-index match supported:{raw_meta['exact_match_supported']}")
|
||
lines.append(f" exact source-index matches: {raw_meta['exact_match_count']}")
|
||
lines.append(f" exported nearby raw rows: {raw_meta['nearby_row_count']}")
|
||
lines.append("")
|
||
lines.append("OUTPUT OF SPECIAL INTEREST")
|
||
lines.append(" invalid_target_neighbor_summary.csv")
|
||
lines.append(" invalid_point_shared_area_neighbors.csv")
|
||
lines.append(" invalid_feature_full_metadata.csv")
|
||
lines.append(" invalid_feature_attributes.csv")
|
||
lines.append(" chart_raw_record_schema.csv")
|
||
lines.append(" invalid_raw_record_trace.csv")
|
||
lines.append(" invalid_geometry_trace.geojson")
|
||
lines.append("")
|
||
lines.append("INTERPRETATION")
|
||
lines.append(
|
||
" 若 invalid point 被多个正常 Area 共享，优先怀疑共享拓扑/ring assembly 表达，"
|
||
"而不是孤立坐标损坏。"
|
||
)
|
||
lines.append(
|
||
" 若内部 TRUE_SEGMENT_CROSSING 没有任何邻接对象共享该点，"
|
||
"则更值得追 source_record / edge ordering。"
|
||
)
|
||
lines.append(
|
||
" chart_raw_record 若不能按 cell_id + source_record_index 直接关联，"
|
||
"说明 DB 当前可能没有保留足够原始 topology，需要回到 CM93 parser/importer 调查。"
|
||
)
|
||
lines.append("")
|
||
lines.append("NO DATABASE WRITE / NO GEOMETRY REPAIR / NO PBF")
|
||
|
||
(outdir / "invalid_trace_summary_v5.txt").write_text(
|
||
"\n".join(lines), encoding="utf-8"
|
||
)
|
||
|
||
conn.rollback()
|
||
|
||
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()
|