pbf/coastline/export_fukuoka_saga_map_assets.py

#!/usr/bin/env python3
"""
NavSea coastline subset export v1.0

把福冈/佐贺海岸样例栅格导出为浏览器可直接加载的 GeoJSON 资产：
- blocked cells (WGS84 GeoJSON)
- coastline lines (WGS84 GeoJSON)
- manifest JSON
"""

from __future__ import annotations

import argparse
import json
import math
import sqlite3
from pathlib import Path

from shapely import from_wkt
from shapely.geometry import mapping, Polygon, LineString


RADIUS = 6378137.0


def mercator_to_lonlat(x: float, y: float) -> tuple[float, float]:
    lon = math.degrees(x / RADIUS)
    lat = math.degrees(2.0 * math.atan(math.exp(y / RADIUS)) - math.pi / 2.0)
    return lon, lat


def cell_polygon_lonlat(minx: float, miny: float, maxx: float, maxy: float) -> Polygon:
    lon1, lat1 = mercator_to_lonlat(minx, miny)
    lon2, lat2 = mercator_to_lonlat(maxx, maxy)
    return Polygon(
        [
            (lon1, lat1),
            (lon2, lat1),
            (lon2, lat2),
            (lon1, lat2),
            (lon1, lat1),
        ]
    )


def line_lonlat_from_wkt(wkt: str) -> LineString:
    geom = from_wkt(wkt)
    coords = [mercator_to_lonlat(x, y) for x, y in geom.coords]
    return LineString(coords)


def main() -> None:
    parser = argparse.ArgumentParser(description="Export coastline grid assets to GeoJSON")
    parser.add_argument(
        "--source-db",
        default="out/coastline/fukuoka_saga_test/fukuoka_saga_coast_grid.sqlite",
        help="输入 SQLite",
    )
    parser.add_argument(
        "--output-dir",
        default="src/pbf/coastline",
        help="输出目录",
    )
    parser.add_argument(
        "--slug",
        default="fukuoka_saga",
        help="输出文件前缀，如 fukuoka_saga / north_kyushu",
    )
    args = parser.parse_args()

    project_root = Path(__file__).resolve().parent.parent
    source_db = Path(args.source_db)
    if not source_db.is_absolute():
        source_db = project_root / source_db
    out_dir = Path(args.output_dir)
    if not out_dir.is_absolute():
        out_dir = project_root / out_dir
    out_dir.mkdir(parents=True, exist_ok=True)

    if not source_db.exists():
        raise SystemExit(f"missing source db: {source_db}")

    conn = sqlite3.connect(source_db)
    conn.row_factory = sqlite3.Row

    manifest = {
        "version": "v1.0",
        "source_db": str(source_db),
        "grid_file": f"{args.slug}_blocked_cells.geojson",
        "coastline_file": f"{args.slug}_coastline_lines.geojson",
        "source_scope": conn.execute("SELECT value FROM meta WHERE key='source_scope'").fetchone()[0],
        "source_code": conn.execute("SELECT value FROM meta WHERE key='source_code'").fetchone()[0],
        "cell_size_m": float(conn.execute("SELECT value FROM meta WHERE key='cell_size_m'").fetchone()[0]),
        "buffer_m": float(conn.execute("SELECT value FROM meta WHERE key='buffer_m'").fetchone()[0]),
        "grid_bounds_mercator": json.loads(
            conn.execute("SELECT value FROM meta WHERE key='grid_bounds_mercator'").fetchone()[0]
        ),
    }

    blocked_features = []
    blocked_rows = conn.execute(
        """
        SELECT cell_id, ix, iy, minx, miny, maxx, maxy, state, state_name, source_mask, blocked_permille
        FROM grid_cell
        WHERE state = 2
        ORDER BY ix, iy
        """
    ).fetchall()

    for row in blocked_rows:
        geom = cell_polygon_lonlat(row["minx"], row["miny"], row["maxx"], row["maxy"])
        blocked_features.append(
            {
                "type": "Feature",
                "properties": {
                    "cell_id": row["cell_id"],
                    "ix": row["ix"],
                    "iy": row["iy"],
                    "state": row["state"],
                    "state_name": row["state_name"],
                    "source_mask": row["source_mask"],
                    "blocked_permille": row["blocked_permille"],
                },
                "geometry": mapping(geom),
            }
        )

    coastline_features = []
    coastline_rows = conn.execute(
        """
        SELECT source_code, prefecture, curve_id, geom_wkt
        FROM coastline_feature
        ORDER BY feature_id
        """
    ).fetchall()

    for row in coastline_rows:
        geom = line_lonlat_from_wkt(row["geom_wkt"])
        coastline_features.append(
            {
                "type": "Feature",
                "properties": {
                    "source_code": row["source_code"],
                    "prefecture": row["prefecture"],
                    "curve_id": row["curve_id"],
                },
                "geometry": mapping(geom),
            }
        )

    blocked_geojson = {
        "type": "FeatureCollection",
        "name": "fukuoka_saga_blocked_cells",
        "features": blocked_features,
    }
    coastline_geojson = {
        "type": "FeatureCollection",
        "name": "fukuoka_saga_coastline_lines",
        "features": coastline_features,
    }

    (out_dir / f"{args.slug}_blocked_cells.geojson").write_text(
        json.dumps(blocked_geojson, ensure_ascii=False),
        encoding="utf-8",
    )
    (out_dir / f"{args.slug}_coastline_lines.geojson").write_text(
        json.dumps(coastline_geojson, ensure_ascii=False),
        encoding="utf-8",
    )

    lonlat_min = mercator_to_lonlat(manifest["grid_bounds_mercator"][0], manifest["grid_bounds_mercator"][1])
    lonlat_max = mercator_to_lonlat(manifest["grid_bounds_mercator"][2], manifest["grid_bounds_mercator"][3])
    manifest["grid_bounds_lonlat"] = [lonlat_min[0], lonlat_min[1], lonlat_max[0], lonlat_max[1]]
    manifest["blocked_count"] = len(blocked_features)
    manifest["coastline_count"] = len(coastline_features)

    (out_dir / f"{args.slug}_manifest.json").write_text(
        json.dumps(manifest, ensure_ascii=False, indent=2),
        encoding="utf-8",
    )

    conn.close()
    print("已导出:")
    print(" -", out_dir / f"{args.slug}_blocked_cells.geojson")
    print(" -", out_dir / f"{args.slug}_coastline_lines.geojson")
    print(" -", out_dir / f"{args.slug}_manifest.json")


if __name__ == "__main__":
    main()