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