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