Copenhagen Phase 1: scaffolding + area + OSM download

- pyproject.toml with all deps (geopandas, contextily, partridge, etc.)
- config/: area.json (3 OSM boundary relations), modes.yaml (subway/s_tog/
  light_rail/regional/bus, S-tog scoped to route=light_rail network=Takst),
  styling.yaml (zorder/widths/colours/palette/mask/labels)
- build_area.py: fetch 3 relations via OSM API, stitch rings, dissolve to
  area.gpkg/area.geojson (City Pass approximation)
- download_osm.py: per-mode Overpass (mirror+backoff, out geom + recurse for
  stops), writes {mode}.gpkg (lines+stops) + query text
- download_gtfs.py / gtfs_to_geopackage.py: stubs (blocked on Rejseplanen)
- data/ gitignored (regenerable); output/ tracked

copenhagen/scripts/build_area.py

@@ -0,0 +1,149 @@
+#!/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()