pbf/coastline/render_fish_port_full_20m_aoi_png.py

#!/usr/bin/env python3
from __future__ import annotations

import argparse
import json
import math
from pathlib import Path

from PIL import Image, ImageDraw
from shapely.geometry import box, shape
from shapely.ops import transform, unary_union
from shapely.prepared import prep


RADIUS = 6378137.0
MAX_MERCATOR_LAT = 85.0511287798066
DEFAULT_BBOX = (129.73110007140482, 32.9464631795233, 131.19176072338317, 34.25430274729531)


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 bbox_intersects(a, b) -> bool:
    return not (a[2] < b[0] or a[0] > b[2] or a[3] < b[1] or a[1] > b[3])


def iter_polygons(geom):
    gt = geom.geom_type
    if gt == "Polygon":
        yield geom
    elif gt == "MultiPolygon":
        for g in geom.geoms:
            yield g


def draw_cell(draw: ImageDraw.ImageDraw, minx: float, maxy: float, px_size: float, x: float, y: float, color):
    x0 = (x - minx) / px_size
    y0 = (maxy - (y + px_size)) / px_size
    x1 = x0 + 1
    y1 = y0 + 1
    draw.rectangle((x0, y0, x1, y1), fill=color)


def main() -> None:
    parser = argparse.ArgumentParser(description="Render full Fukuoka/Saga fish-port 20m red/green cell raster")
    parser.add_argument("--aoi-source", default="src/pbf/port-overlay/kyushu/port_overlay_kyushu_areas.geojson")
    parser.add_argument("--output-png", default="src/pbf/coastline-only/fukuoka_saga_fish_port_navsplit_full_20m_aoi/fish_port_navsplit_fukuoka_saga_full_20m_aoi.png")
    parser.add_argument("--bbox", nargs=4, type=float, default=list(DEFAULT_BBOX))
    parser.add_argument("--cell-size-m", type=float, default=20.0)
    parser.add_argument("--coast-buffer-m", type=float, default=12.0)
    args = parser.parse_args()

    project_root = Path(__file__).resolve().parent.parent
    aoi_path = Path(args.aoi_source)
    if not aoi_path.is_absolute():
        aoi_path = project_root / aoi_path
    if not aoi_path.exists():
        raise SystemExit(f"missing input: {aoi_path}")

    data = json.loads(aoi_path.read_text(encoding="utf-8"))
    bbox = tuple(args.bbox)
    px_size = args.cell_size_m

    selected = []
    for feature in data.get("features", []):
        props = feature.get("properties") or {}
        if props.get("canonical_object_type") != "漁港":
            continue
        geom = shape(feature["geometry"])
        if not bbox_intersects(geom.bounds, bbox):
            continue
        selected.append(transform(lambda x, y, z=None: (mercator_x(x), mercator_y(y)), geom))

    if not selected:
        raise SystemExit("no selected fish-port AOI polygons")

    merged = unary_union(selected)
    prepared = prep(merged)

    minlon, minlat, maxlon, maxlat = bbox
    minx, miny = mercator_x(minlon), mercator_y(minlat)
    maxx, maxy = mercator_x(maxlon), mercator_y(maxlat)
    width = max(1, math.ceil((maxx - minx) / px_size))
    height = max(1, math.ceil((maxy - miny) / px_size))

    image = Image.new("RGBA", (width, height), (0, 0, 0, 0))
    draw = ImageDraw.Draw(image, "RGBA")

    cell_count = 0
    land_count = 0
    sea_count = 0

    # 只对渔港 AOI 内的格子着色，外部保持透明。
    for poly in iter_polygons(merged):
        pb = poly.bounds
        start_x = math.floor(pb[0] / px_size) * px_size
        start_y = math.floor(pb[1] / px_size) * px_size
        end_x = math.ceil(pb[2] / px_size) * px_size
        end_y = math.ceil(pb[3] / px_size) * px_size
        y = start_y
        while y < end_y:
            x = start_x
            while x < end_x:
                cell = box(x, y, x + px_size, y + px_size)
                if prepared.intersects(cell):
                    # 海岸线 buffer 之外视为海面可进入，之内视为陆基不可进入。
                    # 这里暂用 AOI 内的海岸缓冲近似，渲染为红/绿小格。
                    draw_cell(draw, minx, maxy, px_size, x, y, (90, 210, 110, 110))
                    sea_count += 1
                    cell_count += 1
                x += px_size
            y += px_size

    # 让港内不可进入区域更明显：再叠一层海岸缓冲红色边缘。
    # 这里直接沿 AOI 外廓作一个保守缓冲带，后续可再细化到真正 coastline 交集。
    outline = merged.buffer(args.coast_buffer_m, cap_style=2, join_style=2)
    outline_p = prep(outline)
    for poly in iter_polygons(outline):
        pb = poly.bounds
        start_x = math.floor(pb[0] / px_size) * px_size
        start_y = math.floor(pb[1] / px_size) * px_size
        end_x = math.ceil(pb[2] / px_size) * px_size
        end_y = math.ceil(pb[3] / px_size) * px_size
        y = start_y
        while y < end_y:
            x = start_x
            while x < end_x:
                cell = box(x, y, x + px_size, y + px_size)
                if outline_p.intersects(cell):
                    draw_cell(draw, minx, maxy, px_size, x, y, (240, 100, 100, 140))
                    land_count += 1
                x += px_size
            y += px_size

    out_path = Path(args.output_png)
    if not out_path.is_absolute():
        out_path = project_root / out_path
    out_path.parent.mkdir(parents=True, exist_ok=True)
    image.save(out_path)

    manifest = {
        "aoi_source": str(aoi_path),
        "bbox_lonlat": list(bbox),
        "pixel_size_m": px_size,
        "selected_polygon_count": len(selected),
        "image_width": width,
        "image_height": height,
        "cell_count": cell_count,
        "sea_count": sea_count,
        "land_count": land_count,
        "png_file": out_path.name,
    }
    (out_path.parent / "fish_port_navsplit_fukuoka_saga_full_20m_aoi_raster_manifest.json").write_text(
        json.dumps(manifest, ensure_ascii=False, indent=2),
        encoding="utf-8",
    )
    print("已导出:")
    print(" -", out_path)
    print(" -", out_path.parent / "fish_port_navsplit_fukuoka_saga_full_20m_aoi_raster_manifest.json")


if __name__ == "__main__":
    main()