添加全国三层生成与港名查FPC工具
This commit is contained in:
OpenAI Codex
516
coastline/build_japan_hazard_50m_mysql.py
Normal file
516
coastline/build_japan_hazard_50m_mysql.py
Normal file
@@ -0,0 +1,516 @@
|
||||
#!/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()
|
||||
Reference in New Issue
Block a user