#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
CM95 唐津 AOI Presence Diagnostic v2

只读诊断：
1) chart_cell bbox 是否覆盖唐津 z10/881/410；
2) chart_feature bbox 是否命中；
3) 实际 Point/Line/Area geometry 是否通过 MBRIntersects 命中；
4) Sounding point 是否命中；
5) 输出离唐津中心最近的 Cell / Feature，帮助判断 Dataset 1 是否根本没有唐津细节数据。

不修改数据库。
"""

from __future__ import annotations
import argparse
import getpass
import json
import math
from pathlib import Path
from datetime import datetime

import pymysql
from pymysql.cursors import DictCursor

LON_MIN = 129.7265625
LAT_MIN = 33.4314413355753
LON_MAX = 130.078125
LAT_MAX = 33.7243396617476
CENTER_LON = (LON_MIN + LON_MAX) / 2
CENTER_LAT = (LAT_MIN + LAT_MAX) / 2


def args():
    p = argparse.ArgumentParser()
    p.add_argument("--host", default="127.0.0.1")
    p.add_argument("--port", type=int, default=3306)
    p.add_argument("--user", default="root")
    p.add_argument("--password")
    p.add_argument("--database", default="cm93_chart")
    p.add_argument("--dataset-id", type=int, default=1)
    p.add_argument("--output", default="cm95_karatsu_presence_diagnostic_v2.json")
    return p.parse_args()


def fetch(cur, sql, params=()):
    cur.execute(sql, params)
    return list(cur.fetchall())


def scalar(cur, sql, params=()):
    cur.execute(sql, params)
    r = cur.fetchone()
    return None if not r else next(iter(r.values()))


def main():
    a = args()
    if a.password is None:
        a.password = getpass.getpass("MySQL password: ")

    conn = pymysql.connect(
        host=a.host, port=a.port, user=a.user, password=a.password,
        database=a.database, charset="utf8mb4",
        cursorclass=DictCursor, autocommit=False,
        connect_timeout=20, read_timeout=900,
    )

    report = {
        "version": "v2",
        "generated_at": datetime.now().isoformat(timespec="seconds"),
        "dataset_id": a.dataset_id,
        "target": {
            "tile": "z10/881/410",
            "bbox": [LON_MIN, LAT_MIN, LON_MAX, LAT_MAX],
            "center": [CENTER_LON, CENTER_LAT],
        },
        "read_only": True,
    }

    poly_wkt = (
        f"POLYGON(({LON_MIN} {LAT_MIN},"
        f"{LON_MAX} {LAT_MIN},"
        f"{LON_MAX} {LAT_MAX},"
        f"{LON_MIN} {LAT_MAX},"
        f"{LON_MIN} {LAT_MIN}))"
    )

    try:
        with conn.cursor() as cur:
            try:
                cur.execute("SET SESSION TRANSACTION READ ONLY")
                cur.execute("START TRANSACTION READ ONLY")
            except Exception:
                cur.execute("START TRANSACTION")

            report["dataset"] = fetch(cur, """
                SELECT id,name,source_type,source_root,parser_version,normalized_srid
                FROM chart_dataset WHERE id=%s
            """, (a.dataset_id,))

            report["global_feature_bbox"] = fetch(cur, """
                SELECT
                    MIN(bbox_min_lon) min_lon,
                    MIN(bbox_min_lat) min_lat,
                    MAX(bbox_max_lon) max_lon,
                    MAX(bbox_max_lat) max_lat,
                    SUM(bbox_min_lon IS NOT NULL) bbox_rows,
                    COUNT(*) total_rows
                FROM chart_feature
                WHERE dataset_id=%s
            """, (a.dataset_id,))

            report["global_cell_bbox"] = fetch(cur, """
                SELECT
                    MIN(bbox_min_lon) min_lon,
                    MIN(bbox_min_lat) min_lat,
                    MAX(bbox_max_lon) max_lon,
                    MAX(bbox_max_lat) max_lat,
                    SUM(bbox_min_lon IS NOT NULL) bbox_rows,
                    COUNT(*) total_rows
                FROM chart_cell
                WHERE dataset_id=%s
            """, (a.dataset_id,))

            bbox_params = (
                a.dataset_id,
                LON_MIN, LON_MAX, LAT_MIN, LAT_MAX,
            )

            report["cell_bbox_overlap"] = fetch(cur, """
                SELECT id,source_relative_path,scale_code,native_scale,
                       bbox_min_lon,bbox_min_lat,bbox_max_lon,bbox_max_lat,
                       feature_count
                FROM chart_cell
                WHERE dataset_id=%s
                  AND bbox_max_lon >= %s
                  AND bbox_min_lon <= %s
                  AND bbox_max_lat >= %s
                  AND bbox_min_lat <= %s
                ORDER BY native_scale,source_relative_path
            """, bbox_params)

            report["feature_bbox_count"] = scalar(cur, """
                SELECT COUNT(*)
                FROM chart_feature
                WHERE dataset_id=%s
                  AND bbox_max_lon >= %s
                  AND bbox_min_lon <= %s
                  AND bbox_max_lat >= %s
                  AND bbox_min_lat <= %s
            """, bbox_params)

            # Direct actual geometry checks. These bypass chart_feature bbox metadata.
            direct = {}
            for kind, table in [
                ("point", "chart_point_geometry"),
                ("line", "chart_line_geometry"),
                ("area", "chart_area_geometry"),
            ]:
                try:
                    direct[kind] = scalar(cur, f"""
                        SELECT COUNT(*)
                        FROM {table} g
                        JOIN chart_feature f ON f.id=g.feature_id
                        WHERE f.dataset_id=%s
                          AND MBRIntersects(
                              g.geom,
                              ST_GeomFromText(%s,4326)
                          )
                    """, (a.dataset_id, poly_wkt))
                except Exception as e:
                    direct[kind] = {"error": str(e)}

            report["direct_geometry_mbr_counts"] = direct

            try:
                report["sounding_bbox_count"] = scalar(cur, """
                    SELECT COUNT(*)
                    FROM chart_sounding_point s
                    JOIN chart_feature f ON f.id=s.feature_id
                    WHERE f.dataset_id=%s
                      AND ST_X(s.geom) BETWEEN %s AND %s
                      AND ST_Y(s.geom) BETWEEN %s AND %s
                """, (
                    a.dataset_id,
                    LON_MIN, LON_MAX, LAT_MIN, LAT_MAX,
                ))
            except Exception as e:
                report["sounding_bbox_count"] = {"error": str(e)}

            # Nearest cell by bbox center, excluding pathological world-scale cells when possible.
            report["nearest_cells"] = fetch(cur, """
                SELECT id,source_relative_path,scale_code,native_scale,
                       bbox_min_lon,bbox_min_lat,bbox_max_lon,bbox_max_lat,
                       feature_count,
                       SQRT(
                         POW(((bbox_min_lon+bbox_max_lon)/2)-%s,2) +
                         POW(((bbox_min_lat+bbox_max_lat)/2)-%s,2)
                       ) AS center_distance_deg
                FROM chart_cell
                WHERE dataset_id=%s
                  AND bbox_min_lon IS NOT NULL
                  AND bbox_min_lat IS NOT NULL
                  AND bbox_max_lon IS NOT NULL
                  AND bbox_max_lat IS NOT NULL
                  AND (bbox_max_lon-bbox_min_lon) < 20
                  AND (bbox_max_lat-bbox_min_lat) < 20
                ORDER BY center_distance_deg ASC
                LIMIT 30
            """, (CENTER_LON, CENTER_LAT, a.dataset_id))

            report["nearest_features"] = fetch(cur, """
                SELECT f.id,f.cell_id,c.source_relative_path,
                       f.source_record_index,f.object_class_id,f.object_class,
                       f.geometry_type,f.scale_code,f.native_scale,
                       f.bbox_min_lon,f.bbox_min_lat,f.bbox_max_lon,f.bbox_max_lat,
                       SQRT(
                         POW(((f.bbox_min_lon+f.bbox_max_lon)/2)-%s,2) +
                         POW(((f.bbox_min_lat+f.bbox_max_lat)/2)-%s,2)
                       ) AS center_distance_deg
                FROM chart_feature f
                JOIN chart_cell c ON c.id=f.cell_id
                WHERE f.dataset_id=%s
                  AND f.bbox_min_lon IS NOT NULL
                  AND f.bbox_min_lat IS NOT NULL
                  AND f.bbox_max_lon IS NOT NULL
                  AND f.bbox_max_lat IS NOT NULL
                ORDER BY center_distance_deg ASC
                LIMIT 50
            """, (CENTER_LON, CENTER_LAT, a.dataset_id))

            # Coarse 2-degree spatial occupancy by feature bbox center.
            report["coarse_feature_grid"] = fetch(cur, """
                SELECT
                    FLOOR(((bbox_min_lon+bbox_max_lon)/2)/2)*2 AS lon_bin,
                    FLOOR(((bbox_min_lat+bbox_max_lat)/2)/2)*2 AS lat_bin,
                    COUNT(*) AS feature_count
                FROM chart_feature
                WHERE dataset_id=%s
                  AND bbox_min_lon IS NOT NULL
                  AND bbox_min_lat IS NOT NULL
                  AND bbox_max_lon IS NOT NULL
                  AND bbox_max_lat IS NOT NULL
                  AND ((bbox_min_lon+bbox_max_lon)/2) BETWEEN 120 AND 160
                  AND ((bbox_min_lat+bbox_max_lat)/2) BETWEEN 20 AND 50
                GROUP BY lon_bin,lat_bin
                ORDER BY lat_bin,lon_bin
            """, (a.dataset_id,))

            conn.rollback()
    finally:
        conn.close()

    # Derive diagnosis.
    fcnt = report.get("feature_bbox_count")
    d = report.get("direct_geometry_mbr_counts", {})
    direct_sum = 0
    direct_errors = []
    for k, v in d.items():
        if isinstance(v, int):
            direct_sum += v
        else:
            direct_errors.append((k, v))

    if isinstance(fcnt, int) and fcnt > 0:
        diagnosis = "AOI_PRESENT_BY_FEATURE_BBOX"
    elif direct_sum > 0:
        diagnosis = "AOI_PRESENT_IN_GEOMETRY_BUT_FEATURE_BBOX_QUERY_FAILED"
    else:
        diagnosis = "NO_CM95_GEOMETRY_FOUND_IN_KARATSU_AOI"

    report["diagnosis"] = diagnosis
    report["direct_geometry_total"] = direct_sum
    report["direct_geometry_errors"] = direct_errors

    out = Path(a.output).resolve()
    out.write_text(json.dumps(report, ensure_ascii=False, indent=2, default=str), encoding="utf-8")

    print("CM95 唐津 AOI Presence Diagnostic v2")
    print("=" * 72)
    print("Target bbox       :", [LON_MIN, LAT_MIN, LON_MAX, LAT_MAX])
    print("Cell bbox overlap :", len(report["cell_bbox_overlap"]))
    print("Feature bbox count:", report["feature_bbox_count"])
    print("Direct geometry   :", report["direct_geometry_mbr_counts"])
    print("Soundings         :", report["sounding_bbox_count"])
    print("Diagnosis         :", diagnosis)
    print("")
    print("Nearest cells:")
    for r in report["nearest_cells"][:10]:
        print(
            f"  {r['source_relative_path']} scale={r['native_scale']} "
            f"bbox=({r['bbox_min_lon']},{r['bbox_min_lat']},"
            f"{r['bbox_max_lon']},{r['bbox_max_lat']}) "
            f"dist={float(r['center_distance_deg']):.3f}°"
        )
    print("")
    print("Output:", out)
    print("把 JSON 上传给 ChatGPT。")


if __name__ == "__main__":
    main()