jetlag-maps/copenhagen/scripts/build_area.py

#!/usr/bin/env python3
"""Build the City Pass area polygon from three OSM boundary relations.

Fetches relation/{id}/full.json from the OSM API for each relation in
config/area.json, assembles member ways into rings (respecting outer/inner
roles), builds polygons, and dissolves the union into a single multipolygon.

Outputs (in data/processed):
    area.geojson  (EPSG:4326, human-readable + portable)
    area.gpkg     (EPSG:4326, for geopandas/scripts)
"""
import json
import sys
from pathlib import Path

import geopandas as gpd
import requests
from shapely.geometry import MultiPolygon, Polygon
from shapely.ops import unary_union

from _common import CONFIG, PROCESSED, UA

OSM_API = "https://api.openstreetmap.org/api/0.6"


def fetch_relation(rel_id):
    url = f"{OSM_API}/relation/{rel_id}/full.json"
    r = requests.get(url, headers={"User-Agent": UA}, timeout=120)
    r.raise_for_status()
    return r.json()["elements"]


def stitch_rings(ways):
    """Stitch a list of node-coordinate polylines into closed rings.

    `ways` is a list of lists of (lon, lat). Returns a list of closed rings
    (each starting point == ending point).
    """
    eps = 1e-9

    def close(a, b):
        return abs(a[0] - b[0]) < eps and abs(a[1] - b[1]) < eps

    pool = [w[:] for w in ways if w]
    rings = []
    while pool:
        cur = pool.pop(0)
        changed = True
        while changed and not close(cur[0], cur[-1]):
            changed = False
            for i, w in enumerate(pool):
                if close(cur[-1], w[0]):
                    cur = cur + w[1:]
                    pool.pop(i)
                    changed = True
                    break
                if close(cur[-1], w[-1]):
                    cur = cur + w[::-1][1:]
                    pool.pop(i)
                    changed = True
                    break
        rings.append(cur)
    return [r for r in rings if close(r[0], r[-1]) and len(r) >= 4]


def build_polygon(elements):
    nodes = {}
    ways = {}
    relation = None
    for e in elements:
        t = e["type"]
        if t == "node":
            nodes[e["id"]] = (e["lon"], e["lat"])
        elif t == "way":
            ways[e["id"]] = e["nodes"]
        elif t == "relation":
            relation = e
    if relation is None:
        raise ValueError("no relation in element set")

    outer, inner = [], []
    for m in relation["members"]:
        if m["type"] != "way" or m["ref"] not in ways:
            continue
        seq = [nodes[n] for n in ways[m["ref"]] if n in nodes]
        if len(seq) < 2:
            continue
        (inner if m.get("role") == "inner" else outer).append(seq)

    outer_rings = stitch_rings(outer)
    inner_rings = stitch_rings(inner)

    polys = []
    for oring in outer_rings:
        op = Polygon(oring)
        if not op.is_valid:
            op = op.buffer(0)
        holes = []
        for iring in inner_rings:
            ip = Polygon(iring)
            if not ip.is_valid:
                ip = ip.buffer(0)
            if op.contains(ip.representative_point()):
                holes.append(list(ip.exterior.coords))
        try:
            polys.append(Polygon(oring, holes=holes))
        except Exception:
            polys.append(Polygon(oring))
    if not polys:
        raise ValueError("could not assemble any outer rings")
    if len(polys) == 1:
        return polys[0]
    return MultiPolygon(polys)


def main():
    area_cfg = json.loads((CONFIG / "area.json").read_text())
    geoms = []
    for rel in area_cfg["relations"]:
        rid = rel["id"]
        print(f"fetching relation {rid} ({rel['name']})...", flush=True)
        elements = fetch_relation(rid)
        geom = build_polygon(elements)
        if not geom.is_valid:
            geom = geom.buffer(0)
        geoms.append(geom)
        print(f"  -> {geom.geom_type}, area={geom.area:.5f} deg^2", flush=True)

    dissolved = unary_union(geoms)
    if not dissolved.is_valid:
        dissolved = dissolved.buffer(0)

    PROCESSED.mkdir(parents=True, exist_ok=True)
    gdf = gpd.GeoDataFrame(
        {"name": ["City Pass area"]}, geometry=[dissolved], crs="EPSG:4326"
    )
    gdf.to_file(PROCESSED / "area.geojson", driver="GeoJSON")
    gdf.to_file(PROCESSED / "area.gpkg", driver="GPKG", layer="area")
    bounds = dissolved.bounds
    print(
        f"dissolved: {dissolved.geom_type} | bounds "
        f"({bounds[0]:.4f},{bounds[1]:.4f})-({bounds[2]:.4f},{bounds[3]:.4f})",
        flush=True,
    )
    print(f"wrote {PROCESSED/'area.geojson'} and {PROCESSED/'area.gpkg'}")


if __name__ == "__main__":
    main()