#!/usr/bin/env python3 from __future__ import annotations import argparse import datetime as dt import json import math from pathlib import Path from typing import Iterable import mapbox_vector_tile import pymysql from shapely.geometry import box, shape from shapely.prepared import prep DB_NAME = "navsea_japan_coast_grid" DB_USER = "root" DB_PASSWORD = "2chi9ks2" DB_HOST = "localhost" DB_SOCKET = "/tmp/mysql.sock" DEFAULT_TILE_ROOT = Path("/home/wwwroot/pbf-delivery-full-20260418-rebuild") DEFAULT_TILE_Z = 12 DEFAULT_OUT_DIR = "src/pbf/coastline-mysql/japan_national" DEFAULT_SOURCE_DESC = "pbf-delivery-full-20260418-rebuild z12" HAZARD_CELL_M = 50.0 HAZARD_LAYER_NAME = "hazard_50m" HAZARD_DESCRIPTION = "全国 PBF 海上障碍 50m 黄格" HAZARD_LAYERS = ( "navigation_hazard_area", "fixed_fishing_gear_area", "anchor_caution_hazard_area", "navigation_hazard_point", "anchor_caution_hazard_point", "navigation_marks", ) BREAKWATER_LAYER = "baseline_area" POINT_LAYERS = { "navigation_hazard_point", "anchor_caution_hazard_point", "navigation_marks", } HAZARD_CANONICAL_OBJECT_TYPES = { "breakwater", } RADIUS = 6378137.0 MAX_MERCATOR_LAT = 85.0511287798066 def mercator_x(lon: float) -> float: return RADIUS * math.radians(lon) def mercator_y(lat: float) -> float: lat = max(min(lat, MAX_MERCATOR_LAT), -MAX_MERCATOR_LAT) return RADIUS * math.log(math.tan(math.pi / 4.0 + math.radians(lat) / 2.0)) def lon_from_mercator(x: float) -> float: return math.degrees(x / RADIUS) def lat_from_mercator(y: float) -> float: return math.degrees(2.0 * math.atan(math.exp(y / RADIUS)) - math.pi / 2.0) def rect_geojson(min_lon: float, min_lat: float, max_lon: float, max_lat: float) -> dict: return { "type": "Feature", "geometry": { "type": "Polygon", "coordinates": [ [ [min_lon, min_lat], [max_lon, min_lat], [max_lon, max_lat], [min_lon, max_lat], [min_lon, min_lat], ] ], }, } def local_name(tag: str) -> str: return tag.split("}", 1)[-1] def load_json(path: Path) -> dict: return json.loads(path.read_text(encoding="utf-8")) def stream_geojson(path: Path, features: Iterable[dict]) -> int: path.parent.mkdir(parents=True, exist_ok=True) count = 0 with path.open("w", encoding="utf-8") as fh: fh.write('{"type":"FeatureCollection","features":[\n') first = True for feature in features: if not first: fh.write(",\n") fh.write(json.dumps(feature, ensure_ascii=False)) first = False count += 1 fh.write("\n]}\n") return count def mysql_connect(database: str | None = None): kwargs = { "host": DB_HOST, "user": DB_USER, "password": DB_PASSWORD, "charset": "utf8mb4", "autocommit": False, "cursorclass": pymysql.cursors.Cursor, } if DB_SOCKET: kwargs["unix_socket"] = DB_SOCKET if database: kwargs["database"] = database return pymysql.connect(**kwargs) def ensure_database(name: str) -> None: conn = mysql_connect() try: with conn.cursor() as cur: cur.execute( f"CREATE DATABASE IF NOT EXISTS `{name}` " "DEFAULT CHARACTER SET utf8mb4 COLLATE utf8mb4_unicode_ci" ) conn.commit() finally: conn.close() def ensure_schema(conn) -> None: stmts = [ """ CREATE TABLE IF NOT EXISTS navsea_grid_layer_meta ( layer_name VARCHAR(32) NOT NULL PRIMARY KEY, description VARCHAR(255) NOT NULL, source_desc TEXT NOT NULL, cell_size_m DOUBLE NOT NULL, feature_count BIGINT NOT NULL, bbox_min_lon DOUBLE NOT NULL, bbox_min_lat DOUBLE NOT NULL, bbox_max_lon DOUBLE NOT NULL, bbox_max_lat DOUBLE NOT NULL, export_file VARCHAR(255) DEFAULT NULL, created_at TIMESTAMP NOT NULL DEFAULT CURRENT_TIMESTAMP ) ENGINE=InnoDB DEFAULT CHARSET=utf8mb4 """, """ CREATE TABLE IF NOT EXISTS navsea_grid_cell ( id BIGINT UNSIGNED NOT NULL AUTO_INCREMENT PRIMARY KEY, layer_name VARCHAR(32) NOT NULL, cell_id VARCHAR(64) NOT NULL, row_idx INT NOT NULL, col_idx INT NOT NULL, cell_size_m DOUBLE NOT NULL, state_name VARCHAR(32) NOT NULL, class_name VARCHAR(32) NOT NULL, source_name VARCHAR(191) NOT NULL, min_lon DOUBLE NOT NULL, min_lat DOUBLE NOT NULL, max_lon DOUBLE NOT NULL, max_lat DOUBLE NOT NULL, created_at TIMESTAMP NOT NULL DEFAULT CURRENT_TIMESTAMP, UNIQUE KEY uniq_layer_cell (layer_name, cell_id), KEY idx_layer_state (layer_name, state_name), KEY idx_layer_rowcol (layer_name, row_idx, col_idx) ) ENGINE=InnoDB DEFAULT CHARSET=utf8mb4 """, ] with conn.cursor() as cur: for stmt in stmts: cur.execute(stmt) conn.commit() def replace_layer_rows(conn, layer_name: str) -> None: with conn.cursor() as cur: cur.execute("DELETE FROM navsea_grid_cell WHERE layer_name=%s", (layer_name,)) cur.execute("DELETE FROM navsea_grid_layer_meta WHERE layer_name=%s", (layer_name,)) conn.commit() def insert_cells(cur, rows: list[tuple]) -> None: cur.executemany( """ INSERT INTO navsea_grid_cell (layer_name, cell_id, row_idx, col_idx, cell_size_m, state_name, class_name, source_name, min_lon, min_lat, max_lon, max_lat) VALUES (%s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s) """, rows, ) def upsert_meta( cur, *, layer_name: str, description: str, source_desc: str, cell_size_m: float, feature_count: int, bbox: tuple[float, float, float, float], export_file: str, ) -> None: cur.execute( """ INSERT INTO navsea_grid_layer_meta (layer_name, description, source_desc, cell_size_m, feature_count, bbox_min_lon, bbox_min_lat, bbox_max_lon, bbox_max_lat, export_file) VALUES (%s, %s, %s, %s, %s, %s, %s, %s, %s, %s) ON DUPLICATE KEY UPDATE description=VALUES(description), source_desc=VALUES(source_desc), cell_size_m=VALUES(cell_size_m), feature_count=VALUES(feature_count), bbox_min_lon=VALUES(bbox_min_lon), bbox_min_lat=VALUES(bbox_min_lat), bbox_max_lon=VALUES(bbox_max_lon), bbox_max_lat=VALUES(bbox_max_lat), export_file=VALUES(export_file) """, ( layer_name, description, source_desc, cell_size_m, feature_count, bbox[0], bbox[1], bbox[2], bbox[3], export_file, ), ) def tile_mercator_bounds(z: int, x: int, y: int) -> tuple[float, float, float, float]: n = 2**z lon_left = x / n * 360.0 - 180.0 lon_right = (x + 1) / n * 360.0 - 180.0 lat_top = math.degrees(math.atan(math.sinh(math.pi * (1 - 2 * y / n)))) lat_bottom = math.degrees(math.atan(math.sinh(math.pi * (1 - 2 * (y + 1) / n)))) return mercator_x(lon_left), mercator_y(lat_bottom), mercator_x(lon_right), mercator_y(lat_top) def transform_geom_from_tile(geom, bounds_m: tuple[float, float, float, float], extent: int): minx, miny, maxx, maxy = bounds_m dx = maxx - minx dy = maxy - miny def coord(x: float, y: float): mx = minx + (x / extent) * dx my = miny + (y / extent) * dy return mx, my def walk(obj): if isinstance(obj[0], (int, float)): return coord(obj[0], obj[1]) return [walk(item) for item in obj] return walk(geom) def iter_source_tiles(tile_root: Path, zoom: int) -> list[Path]: zoom_root = tile_root / str(zoom) if not zoom_root.exists(): raise SystemExit(f"missing tile zoom root: {zoom_root}") return sorted(zoom_root.glob("*/*.pbf")) def iter_hazard_cells(tile_root: Path, zoom: int) -> Iterable[tuple[int, int]]: hazard_cells: set[tuple[int, int]] = set() tiles = iter_source_tiles(tile_root, zoom) for index, p in enumerate(tiles, start=1): try: tx = int(p.parent.name) ty = int(p.stem) except ValueError: continue try: tile = mapbox_vector_tile.decode(p.read_bytes()) except Exception as exc: print(f"[hazard_50m] 跳过无法解析的 tile: {p} ({exc.__class__.__name__})") continue if index % 5000 == 0: print(f"[hazard_50m] 已扫描 tile {index}/{len(tiles)}: {p.parent.parent.name}/{tx}/{ty}.pbf") for layer_name in HAZARD_LAYERS + (BREAKWATER_LAYER,): layer = tile.get(layer_name) if not layer: continue extent = int(layer.get("extent") or 1048576) minx, miny, maxx, maxy = tile_mercator_bounds(zoom, tx, ty) for feature in layer.get("features", []): props = feature.get("properties") or {} if layer_name == BREAKWATER_LAYER and props.get("canonical_object_type") not in HAZARD_CANONICAL_OBJECT_TYPES: continue if ( props.get("canonical_object_type") == "fish_reef" or props.get("chart_symbol_code") == "fish_reef" or props.get("class_name") == "魚礁" ): continue geom = feature.get("geometry") or {} coords = geom.get("coordinates") if not coords: continue shp_coords = transform_geom_from_tile(coords, (minx, miny, maxx, maxy), extent) shp = shape({"type": geom.get("type"), "coordinates": shp_coords}) if shp.is_empty: continue if layer_name in POINT_LAYERS: shp = shp.buffer(25.0) minx2, miny2, maxx2, maxy2 = shp.bounds start_x = math.floor(minx2 / HAZARD_CELL_M) * HAZARD_CELL_M start_y = math.floor(miny2 / HAZARD_CELL_M) * HAZARD_CELL_M end_x = math.ceil(maxx2 / HAZARD_CELL_M) * HAZARD_CELL_M end_y = math.ceil(maxy2 / HAZARD_CELL_M) * HAZARD_CELL_M prep_geom = prep(shp) cy = start_y while cy < end_y: cx = start_x while cx < end_x: cell = box(cx, cy, cx + HAZARD_CELL_M, cy + HAZARD_CELL_M) if prep_geom.intersects(cell): hazard_cells.add((int(round(cx / HAZARD_CELL_M)), int(round(cy / HAZARD_CELL_M)))) cx += HAZARD_CELL_M cy += HAZARD_CELL_M return sorted(hazard_cells, key=lambda item: (item[1], item[0])) def cell_bbox_from_mercator(ix: int, iy: int, cell_m: float) -> tuple[float, float, float, float]: minx = ix * cell_m miny = iy * cell_m maxx = minx + cell_m maxy = miny + cell_m return ( lon_from_mercator(minx), lat_from_mercator(miny), lon_from_mercator(maxx), lat_from_mercator(maxy), ) def import_hazard_layer( conn, *, layer_name: str, description: str, source_desc: str, cell_size_m: float, tile_root: Path, zoom: int, export_path: Path, ) -> tuple[int, tuple[float, float, float, float]]: cells = iter_hazard_cells(tile_root, zoom) batch: list[tuple] = [] count = 0 export_bbox = [float("inf"), float("inf"), float("-inf"), float("-inf")] with conn.cursor() as cur: print(f"[{layer_name}] 开始扫描全国 PBF 危险层：{tile_root} / z{zoom}") with export_path.open("w", encoding="utf-8") as fh: fh.write('{"type":"FeatureCollection","features":[\n') first = True for ix, iy in cells: min_lon, min_lat, max_lon, max_lat = cell_bbox_from_mercator(ix, iy, cell_size_m) export_bbox[0] = min(export_bbox[0], min_lon) export_bbox[1] = min(export_bbox[1], min_lat) export_bbox[2] = max(export_bbox[2], max_lon) export_bbox[3] = max(export_bbox[3], max_lat) cell_id = f"{ix}:{iy}" batch.append( ( layer_name, cell_id, iy, ix, cell_size_m, "HAZARD_50M", "HAZARD", source_desc, min_lon, min_lat, max_lon, max_lat, ) ) feature = { "type": "Feature", "properties": { "layer_name": layer_name, "cell_id": cell_id, "row": iy, "col": ix, "cell_size_m": cell_size_m, "state_name": "HAZARD_50M", "class_name": "HAZARD", "source_name": source_desc, }, "geometry": rect_geojson(min_lon, min_lat, max_lon, max_lat)["geometry"], } if not first: fh.write(",\n") fh.write(json.dumps(feature, ensure_ascii=False)) first = False count += 1 if len(batch) >= 5000: insert_cells(cur, batch) conn.commit() batch.clear() if count % 20000 == 0: print(f"[{layer_name}] 已处理 {count} 个格子 ...") if batch: insert_cells(cur, batch) conn.commit() fh.write("\n]}\n") if count == 0: raise RuntimeError(f"{layer_name} 扫描结果为空，请检查 tile_root={tile_root} zoom={zoom}") print(f"[{layer_name}] 导入完成：{count} 条，导出 {export_path}") return count, (export_bbox[0], export_bbox[1], export_bbox[2], export_bbox[3]) def parse_args() -> argparse.Namespace: parser = argparse.ArgumentParser(description="重算全国 50x50 海上障碍格并写入统一 MySQL 库。") parser.add_argument("--db-name", default=DB_NAME, help="MySQL 数据库名") parser.add_argument("--tile-root", type=Path, default=DEFAULT_TILE_ROOT, help="全国 PBF 根目录") parser.add_argument("--zoom", type=int, default=DEFAULT_TILE_Z, help="PBF 瓦片 zoom，默认 12") parser.add_argument("--out-dir", default=DEFAULT_OUT_DIR, help="GeoJSON 输出目录") parser.add_argument("--source-desc", default=DEFAULT_SOURCE_DESC, help="写入元数据的来源说明") return parser.parse_args() def main() -> None: args = parse_args() project_root = Path(__file__).resolve().parent.parent out_dir = Path(args.out_dir) if not out_dir.is_absolute(): out_dir = project_root / out_dir out_dir.mkdir(parents=True, exist_ok=True) ensure_database(args.db_name) conn = mysql_connect(args.db_name) try: print(f"数据库 {args.db_name} 已连接，准备初始化表结构 ...") ensure_schema(conn) replace_layer_rows(conn, HAZARD_LAYER_NAME) hazard_export = out_dir / "hazard_50m_grid.geojson" hazard_count, hazard_bbox = import_hazard_layer( conn, layer_name=HAZARD_LAYER_NAME, description=HAZARD_DESCRIPTION, source_desc=args.source_desc, cell_size_m=HAZARD_CELL_M, tile_root=args.tile_root, zoom=args.zoom, export_path=hazard_export, ) with conn.cursor() as cur: upsert_meta( cur, layer_name=HAZARD_LAYER_NAME, description=HAZARD_DESCRIPTION, source_desc=args.source_desc, cell_size_m=HAZARD_CELL_M, feature_count=hazard_count, bbox=hazard_bbox, export_file=str(hazard_export.relative_to(project_root)), ) conn.commit() manifest = { "database": args.db_name, "generated_at": dt.datetime.now().isoformat(timespec="seconds"), "layer": { "name": HAZARD_LAYER_NAME, "count": hazard_count, "bbox_lonlat": list(hazard_bbox), "cell_size_m": HAZARD_CELL_M, "export_file": str(hazard_export.relative_to(project_root)), "source_desc": args.source_desc, "tile_root": str(args.tile_root), "zoom": args.zoom, }, } manifest_path = out_dir / "manifest.json" manifest_path.write_text(json.dumps(manifest, ensure_ascii=False, indent=2), encoding="utf-8") print("已导入 MySQL 数据库:", args.db_name) print(" - hazard_50m:", hazard_count, hazard_export) print(" - manifest:", manifest_path) finally: conn.close() if __name__ == "__main__": main()