305 lines
11 KiB
Python
305 lines
11 KiB
Python
#!/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()
|