#!/usr/bin/env python3 # -*- coding: utf-8 -*- r""" CM95 / CM93 MySQL → PBF 测试 AOI 自动候选选择器 v1 目标: 从已导入 MySQL 的 CM93 Dataset 中,自动寻找适合做 CM95→NavSea PBF 第一轮对照测试的 z10 XYZ tile。 优先选择同时包含: - 陆地/岸线 - 水深区 - 等深线 - 灯标/航标 - 危险物 - 较细比例尺数据 脚本只读数据库,不修改任何内容。 执行: python .\cm95_select_pbf_test_aoi_v1.py ` --host 127.0.0.1 ` --port 3306 ` --user root ` --database cm93_chart ` --dataset-id 1 输出: - cm95_pbf_aoi_candidates_v1.csv - cm95_pbf_aoi_candidates_v1.json - cm95_pbf_aoi_candidates_v1.txt """ from __future__ import annotations import argparse import csv import getpass import json import math from collections import Counter, defaultdict from datetime import datetime from pathlib import Path import pymysql from pymysql.cursors import SSDictCursor Z = 10 LAND = {"LNDARE", "COALNE", "SLCONS"} DEPTH = {"DEPARE"} CONTOUR = {"DEPCNT"} LIGHT_NAV_EXACT = { "LIGHTS", "LITFLT", "LITVES", "LNDMRK", "FOGSIG", "RTPBCN", "TOPMAR", "DAYMAR", } LIGHT_NAV_PREFIX = ("BOY", "BCN") HAZARD = { "OBSTRN", "UWTROC", "WRECKS", "SBDARE", "CTNARE", "MORFAC", } # native_scale 数值越小,通常越详细。 DETAIL_SCALE_LIMIT = 50000 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", default=None) p.add_argument("--database", default="cm93_chart") p.add_argument("--dataset-id", type=int, default=1) p.add_argument("--min-lon", type=float, default=120.0) p.add_argument("--max-lon", type=float, default=160.0) p.add_argument("--min-lat", type=float, default=40.0) p.add_argument("--max-lat", type=float, default=50.0) p.add_argument("--top", type=int, default=40) p.add_argument("--output-dir", default=".") return p.parse_args() def lonlat_to_tile(lon, lat, z): lat = max(min(lat, 85.05112878), -85.05112878) n = 2 ** z x = int((lon + 180.0) / 360.0 * n) lat_rad = math.radians(lat) y = int( (1.0 - math.asinh(math.tan(lat_rad)) / math.pi) / 2.0 * n ) x = max(0, min(n - 1, x)) y = max(0, min(n - 1, y)) return x, y def tile_bounds(z, x, y): n = 2 ** z min_lon = x / n * 360.0 - 180.0 max_lon = (x + 1) / n * 360.0 - 180.0 def lat(ty): return math.degrees( math.atan(math.sinh(math.pi * (1 - 2 * ty / n))) ) max_lat = lat(y) min_lat = lat(y + 1) return min_lon, min_lat, max_lon, max_lat def is_navmark(cls): return cls in LIGHT_NAV_EXACT or cls.startswith(LIGHT_NAV_PREFIX) def classify(cls): flags = set() if cls in LAND: flags.add("land") if cls in DEPTH: flags.add("depth") if cls in CONTOUR: flags.add("contour") if is_navmark(cls): flags.add("nav") if cls in HAZARD: flags.add("hazard") return flags def main(): a = args() if a.password is None: a.password = getpass.getpass("MySQL password: ") outdir = Path(a.output_dir).resolve() outdir.mkdir(parents=True, exist_ok=True) conn = pymysql.connect( host=a.host, port=a.port, user=a.user, password=a.password, database=a.database, charset="utf8mb4", cursorclass=SSDictCursor, autocommit=False, connect_timeout=20, read_timeout=900, ) tiles = {} try: with conn.cursor() as cur: try: cur.execute("SET SESSION TRANSACTION READ ONLY") except Exception: pass try: cur.execute("START TRANSACTION READ ONLY") except Exception: cur.execute("START TRANSACTION") cur.execute( """ SELECT object_class, geometry_type, native_scale, bbox_min_lon, bbox_min_lat, bbox_max_lon, bbox_max_lat 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_max_lon >= %s AND bbox_min_lon <= %s AND bbox_max_lat >= %s AND bbox_min_lat <= %s """, ( a.dataset_id, a.min_lon, a.max_lon, a.min_lat, a.max_lat, ), ) scanned = 0 for row in cur: scanned += 1 lon = ( float(row["bbox_min_lon"]) + float(row["bbox_max_lon"]) ) / 2.0 lat = ( float(row["bbox_min_lat"]) + float(row["bbox_max_lat"]) ) / 2.0 if not ( a.min_lon <= lon <= a.max_lon and a.min_lat <= lat <= a.max_lat ): continue x, y = lonlat_to_tile(lon, lat, Z) key = (x, y) if key not in tiles: tiles[key] = { "feature_count": 0, "detail_feature_count": 0, "classes": Counter(), "geometry": Counter(), "flags": Counter(), "native_scales": Counter(), } t = tiles[key] cls = str(row["object_class"] or "") geom = str(row["geometry_type"] or "") scale = int(row["native_scale"] or 999999999) t["feature_count"] += 1 t["classes"][cls] += 1 t["geometry"][geom] += 1 t["native_scales"][scale] += 1 if scale <= DETAIL_SCALE_LIMIT: t["detail_feature_count"] += 1 for flag in classify(cls): t["flags"][flag] += 1 conn.rollback() finally: conn.close() rows = [] for (x, y), t in tiles.items(): flags = t["flags"] required_present = sum( 1 for name in ("land", "depth", "contour", "nav", "hazard") if flags[name] > 0 ) hard_gate_present = ( int(flags["nav"] > 0) + int(flags["hazard"] > 0) + int(flags["depth"] > 0) ) # 先奖励类别闭环,再奖励细尺度和总要素。 score = ( required_present * 1_000_000 + hard_gate_present * 250_000 + min(t["detail_feature_count"], 10000) * 100 + min(t["feature_count"], 100000) ) min_lon, min_lat, max_lon, max_lat = tile_bounds(Z, x, y) rows.append({ "score": score, "z": Z, "x": x, "y": y, "tile_path": f"{Z}/{x}/{y}.pbf", "center_lon": (min_lon + max_lon) / 2, "center_lat": (min_lat + max_lat) / 2, "min_lon": min_lon, "min_lat": min_lat, "max_lon": max_lon, "max_lat": max_lat, "feature_count": t["feature_count"], "detail_feature_count_le_50000": t["detail_feature_count"], "required_groups_present": required_present, "land_count": flags["land"], "depth_count": flags["depth"], "contour_count": flags["contour"], "navmark_light_count": flags["nav"], "hazard_count": flags["hazard"], "point_count": t["geometry"]["POINT"], "line_count": t["geometry"]["LINE"], "area_count": t["geometry"]["AREA"], "top_classes": [ {"class": k, "count": v} for k, v in t["classes"].most_common(15) ], "native_scales": [ {"scale": k, "count": v} for k, v in sorted(t["native_scales"].items()) ], }) rows.sort( key=lambda r: ( r["score"], r["detail_feature_count_le_50000"], r["feature_count"], ), reverse=True, ) top = rows[: a.top] csv_path = outdir / "cm95_pbf_aoi_candidates_v1.csv" json_path = outdir / "cm95_pbf_aoi_candidates_v1.json" txt_path = outdir / "cm95_pbf_aoi_candidates_v1.txt" csv_fields = [ "rank", "score", "z", "x", "y", "tile_path", "center_lon", "center_lat", "min_lon", "min_lat", "max_lon", "max_lat", "feature_count", "detail_feature_count_le_50000", "required_groups_present", "land_count", "depth_count", "contour_count", "navmark_light_count", "hazard_count", "point_count", "line_count", "area_count", "top_classes", "native_scales", ] with csv_path.open("w", encoding="utf-8-sig", newline="") as f: w = csv.DictWriter(f, fieldnames=csv_fields) w.writeheader() for i, row in enumerate(top, 1): r = dict(row) r["rank"] = i r["top_classes"] = json.dumps( r["top_classes"], ensure_ascii=False ) r["native_scales"] = json.dumps( r["native_scales"], ensure_ascii=False ) w.writerow(r) json_path.write_text( json.dumps( { "version": "v1", "generated_at": datetime.now().isoformat(timespec="seconds"), "dataset_id": a.dataset_id, "search_bbox": [ a.min_lon, a.min_lat, a.max_lon, a.max_lat ], "z": Z, "candidate_count": len(rows), "top": top, }, ensure_ascii=False, indent=2, ), encoding="utf-8", ) lines = [ "CM95 → NavSea PBF 第一测试 AOI 候选 v1", "=" * 86, f"Dataset ID: {a.dataset_id}", f"扫描范围: lon {a.min_lon}..{a.max_lon}, lat {a.min_lat}..{a.max_lat}", "", ] for i, r in enumerate(top[:15], 1): lines += [ ( f"#{i:02d} z{r['z']}/{r['x']}/{r['y']} " f"center=({r['center_lon']:.5f},{r['center_lat']:.5f}) " f"features={r['feature_count']} " f"detail={r['detail_feature_count_le_50000']} " f"groups={r['required_groups_present']}/5" ), ( f" land={r['land_count']} depth={r['depth_count']} " f"contour={r['contour_count']} " f"nav={r['navmark_light_count']} " f"hazard={r['hazard_count']}" ), ] lines += [ "", "下一步:把 JSON 或 CSV 上传给 ChatGPT。", "ChatGPT 会选一个最适合的 tile,并告诉你需要提供的现有 NavSea PBF z10/z11/z12 路径。", ] txt_path.write_text("\n".join(lines), encoding="utf-8") print("\n".join(lines)) print("") print("输出:") print(csv_path) print(json_path) print(txt_path) if __name__ == "__main__": main()