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/_common.py

@@ -0,0 +1,14 @@
+"""Shared paths and constants for the Copenhagen map scripts."""
+from pathlib import Path
+
+HERE = Path(__file__).resolve().parent
+CPH = HERE.parent
+CONFIG = CPH / "config"
+DATA = CPH / "data"
+RAW = DATA / "raw"
+OSM_RAW = RAW / "osm"
+GTFS_RAW = RAW / "gtfs"
+PROCESSED = DATA / "processed"
+OUTPUT = CPH / "output"
+
+UA = "jetlag-maps/0.1 (https://github.com/marvin/jetlag-maps)"

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()

copenhagen/scripts/download_gtfs.py

@@ -0,0 +1,45 @@
+#!/usr/bin/env python3
+"""Download the Rejseplanen GTFS feed for the Copenhagen area.
+
+STATUS: BLOCKED. Rejseplanen (journey planner) does not currently expose a
+public GTFS download. This script documents the intended flow and exits
+non-zero with guidance.
+
+When a feed URL becomes available, set it here and the script will download
+and unzip into data/raw/gtfs/.
+"""
+import sys
+import urllib.request
+import zipfile
+from pathlib import Path
+
+from _common import GTFS_RAW
+
+# Fill in once access is granted. Likely candidates:
+#   - Rejseplanen / DOT open-data portal
+#   - a static GTFS zip provided on request
+GTFS_URL = None  # e.g. "https://.../rejseplanen.zip"
+
+
+def main():
+    if not GTFS_URL:
+        print(
+            "download_gtfs: BLOCKED — no GTFS feed URL configured.\n"
+            "Set GTFS_URL in this script once Rejseplanen grants access, then re-run.\n"
+            "The OSM-only pipeline (build_area -> download_osm -> prepare -> render)\n"
+            "is fully functional without GTFS; GTFS only enriches colours/stops.",
+            file=sys.stderr,
+        )
+        sys.exit(2)
+    GTFS_RAW.mkdir(parents=True, exist_ok=True)
+    zip_path = GTFS_RAW / "feed.zip"
+    print(f"downloading {GTFS_URL} -> {zip_path} ...", flush=True)
+    urllib.request.urlretrieve(GTFS_URL, zip_path)
+    with zipfile.ZipFile(zip_path) as z:
+        z.extractall(GTFS_RAW)
+    zip_path.unlink(missing_ok=True)
+    print(f"extracted GTFS txt files into {GTFS_RAW}")
+
+
+if __name__ == "__main__":
+    main()

copenhagen/scripts/download_osm.py

@@ -0,0 +1,153 @@
+#!/usr/bin/env python3
+"""Download transit route lines + stops from Overpass, per mode.
+
+For each mode in config/modes.yaml, runs an Overpass query (with mirror
+fallback + exponential backoff) and writes:
+    data/raw/osm/{mode}.gpkg          (layers: "lines", "stops", EPSG:4326)
+    data/raw/osm/{mode}.overpassql    (the exact query text)
+
+Lines: one MultiLineString per route relation (concatenation of member way
+geometries), carrying ref/name/network/colour/operator/route/mode tags.
+Stops: one Point per stop/platform node member, carrying ref/name/mode.
+"""
+import sys
+import time
+from pathlib import Path
+
+import geopandas as gpd
+import requests
+import yaml
+from shapely.geometry import LineString, MultiLineString, Point
+
+from _common import CONFIG, OSM_RAW, PROCESSED, UA
+
+MIRRORS = [
+    "https://overpass-api.de/api/interpreter",
+    "https://overpass.kumi.systems/api/interpreter",
+    "https://overpass.private.coffee/api/interpreter",
+]
+
+PLATFORM_ROLES = {"platform", "platform_entry", "platform_exit"}
+STOP_ROLES = {"stop", "stop_entry", "stop_exit", "platform", "platform_entry", "platform_exit"}
+
+
+def build_query(mode_cfg, bbox):
+    s, n, w, e = bbox
+    filters = []
+    for k, v in mode_cfg["osm"].items():
+        if k.endswith("_neq"):
+            filters.append(f'["{k[:-4]}"!="{v}"]')
+        else:
+            filters.append(f'["{k}"="{v}"]')
+    filt = "".join(filters)
+    return f"""[out:json][timeout:180];
+relation{filt}({s},{w},{n},{e});
+out geom;
+>;
+out body qt;
+"""
+
+
+def overpass(ql):
+    last_err = None
+    for mi, mirror in enumerate(MIRRORS):
+        for attempt in range(4):
+            try:
+                r = requests.post(
+                    mirror,
+                    data={"data": ql},
+                    headers={"User-Agent": UA},
+                    timeout=300,
+                )
+                if r.status_code == 429 or r.status_code >= 500:
+                    raise RuntimeError(f"HTTP {r.status_code}")
+                r.raise_for_status()
+                return r.json()
+            except Exception as e:
+                last_err = e
+                wait = 2 ** (attempt + mi)
+                print(f"  [{mirror}] attempt {attempt+1} failed: {e}; retry in {wait}s", flush=True)
+                time.sleep(wait)
+    raise RuntimeError(f"all mirrors failed: {last_err}")
+
+
+def parse(data, mode):
+    nodes = {}
+    for e in data["elements"]:
+        if e["type"] == "node":
+            nodes[e["id"]] = e
+
+    lines, stops = [], []
+    stop_seen = set()
+    for e in data["elements"]:
+        if e["type"] != "relation" or "tags" not in e:
+            continue
+        tags = e["tags"]
+        if tags.get("route") is None:
+            continue
+        segs = []
+        for m in e["members"]:
+            if m["type"] == "way" and "geometry" in m and m.get("role") not in PLATFORM_ROLES:
+                coords = [(g["lon"], g["lat"]) for g in m["geometry"]]
+                if len(coords) >= 2:
+                    segs.append(LineString(coords))
+        geom = MultiLineString(segs) if len(segs) > 1 else (segs[0] if segs else None)
+        if geom is None:
+            continue
+        lines.append({
+            "mode": mode,
+            "ref": tags.get("ref"),
+            "name": tags.get("name"),
+            "network": tags.get("network"),
+            "colour": tags.get("colour"),
+            "operator": tags.get("operator"),
+            "route": tags.get("route"),
+            "geometry": geom,
+        })
+        for m in e["members"]:
+            if m["type"] != "node" or m.get("role") not in STOP_ROLES:
+                continue
+            nid = m["ref"]
+            if nid not in nodes or nid in stop_seen:
+                continue
+            nd = nodes[nid]
+            if "lat" not in nd or "lon" not in nd:
+                continue
+            stop_seen.add(nid)
+            ntags = nd.get("tags", {})
+            stops.append({
+                "mode": mode,
+                "name": ntags.get("name") or ntags.get("public_transport") or "stop",
+                "ref": ntags.get("ref"),
+                "public_transport": ntags.get("public_transport"),
+                "railway": ntags.get("railway"),
+                "geometry": Point(nd["lon"], nd["lat"]),
+            })
+    return lines, stops
+
+
+def main():
+    modes = yaml.safe_load((CONFIG / "modes.yaml").read_text())["modes"]
+    area = gpd.read_file(PROCESSED / "area.gpkg")
+    w, s, e, n = area.total_bounds  # (minx, miny, maxx, maxy) = (west, south, east, north)
+
+    OSM_RAW.mkdir(parents=True, exist_ok=True)
+    for mode, cfg in modes.items():
+        ql = build_query(cfg, (s, n, w, e))
+        (OSM_RAW / f"{mode}.overpassql").write_text(ql)
+        print(f"[{mode}] querying Overpass bbox=({s:.4f},{w:.4f},{n:.4f},{e:.4f})...", flush=True)
+        data = overpass(ql)
+        lines, stops = parse(data, mode)
+        print(f"  -> {len(lines)} lines, {len(stops)} stops", flush=True)
+        if not lines:
+            continue
+        lgdf = gpd.GeoDataFrame(lines, crs="EPSG:4326")
+        lgdf.to_file(OSM_RAW / f"{mode}.gpkg", driver="GPKG", layer="lines")
+        if stops:
+            sgdf = gpd.GeoDataFrame(stops, crs="EPSG:4326")
+            sgdf.to_file(OSM_RAW / f"{mode}.gpkg", driver="GPKG", layer="stops")
+    print("done.")
+
+
+if __name__ == "__main__":
+    main()

copenhagen/scripts/gtfs_to_geopackage.py

@@ -0,0 +1,42 @@
+#!/usr/bin/env python3
+"""Convert raw GTFS txt files into portable GeoPackage layers + a colour table.
+
+STATUS: BLOCKED on download_gtfs.py (no feed URL yet).
+
+Intended outputs (data/processed):
+    gtfs_shapes.gpkg   (shapes.txt -> LineString per shape_id, EPSG:4326)
+    gtfs_stops.gpkg    (stops.txt -> Point per stop, EPSG:4326)
+    route_colors.csv   (route_id, route_short_name, route_type, route_color)
+
+Uses partridge for fast GTFS parsing. prepare.py reads these when present and
+falls back to OSM-only data when absent.
+"""
+import sys
+from pathlib import Path
+
+from _common import GTFS_RAW, PROCESSED
+
+
+def main():
+    if not (GTFS_RAW / "routes.txt").exists():
+        print(
+            "gtfs_to_geopackage: BLOCKED — no GTFS data found in "
+            f"{GTFS_RAW}. Run download_gtfs.py first once a feed URL is set.",
+            file=sys.stderr,
+        )
+        sys.exit(2)
+
+    # TODO (when GTFS available):
+    #   import partridge as pt
+    #   import geopandas as gpd
+    #   from shapely.geometry import LineString, Point
+    #   feed = pt.load_geo_feed(str(GTFS_RAW))
+    #   shapes -> gtfs_shapes.gpkg (LineString per shape_id)
+    #   stops  -> gtfs_stops.gpkg
+    #   routes -> route_colors.csv (route_id, route_short_name, route_type, route_color)
+    print("gtfs_to_geopackage: not yet implemented (GTFS feed unavailable).")
+    sys.exit(2)
+
+
+if __name__ == "__main__":
+    main()