#!/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()