Source Copenhagen transit data from Rejseplanen GTFS

- Replace OSM Overpass transit data with the Rejseplanen GTFS feed
  (routes keyed by (agency, short name); styles from modes.yaml)
- Draw one shape per (style, ref, direction), choosing the shape that
  serves the most in-area stops so lines pass the stops we show
- Collapse stops by names (verified unambiguous); prune stops whose
  serving refs have no drawn line within 300 m
- Patch Københavns Havn + Nordhavn into the area polygon so ferry
  routes and sub-harbour metro tunnels survive clipping; 100 m buffer
  closes relation boundary slivers
- Drop OSM download pipeline; keep tiled basemap

copenhagen/config/area.json

@@ -5,5 +5,22 @@
     {"name": "Kobenhavns Kommune", "id": 2192363},
     {"name": "Frederiksberg Kommune", "id": 2186660},
     {"name": "Amager", "id": 5175924}
+  ],
+  "waterways": [
+    {
+      "name": "Kobenhavns Havn + Nordhavn",
+      "comment": "The harbour channel and Nordhavn basin belong to the City Pass zone in any practical sense: the metro tunnels under it (M1/M2/M4), harbour ferries 991/992 sail it, and the fade-out should not start mid-harbour. The kommune boundary relations exclude all water, so this polygon is patched in.",
+      "ring": [
+        [12.532, 55.640],
+        [12.545, 55.666],
+        [12.565, 55.680],
+        [12.585, 55.715],
+        [12.612, 55.717],
+        [12.632, 55.685],
+        [12.632, 55.660],
+        [12.600, 55.636],
+        [12.560, 55.634]
+      ]
+    }
   ]
 }

copenhagen/config/modes.yaml

@@ -1,40 +1,34 @@
-# Mapping of map modes to OSM route filters, GTFS route types, and style keys.
-# `style` references a key in styling.yaml.
+# Mapping of map styles to GTFS agency / route_type filters.
 #
-# OSM filter keys map directly to Overpass tag filters. A key suffixed with
-# `_neq` emits a negated tag filter (["k"!="v"]); any other key emits an
-# equality filter (["k"="v"]). The bbox is appended by download_osm.py.
+# Used by prepare.py to classify GTFS routes (lines layer) and the routes
+# serving each stop pole (stops/stations layers). Style keys are referenced
+# by styling.yaml (palette, zorder, widths).
 #
-# S-tog are tagged route=light_rail in OSM (NOT route=train); they are scoped
-# by network="Takst Sjælland". The separate `light_rail` mode catches other
-# light rail (e.g. Hovedstadens Letbane, opened Aug 2026) via network!=Takst.
-# GTFS route_type: 0=tram,1=subway,2=rail,3=bus,4=ferry.
+# `agencies` matches GTFS agency_name from agency.txt.
+# `route_types` optionally restricts GTFS route_type values, including the
+# extended codes used in this feed (109 = S-tog suburban rail, 700/715 =
+# Movia bus service types). Without `route_types`, all of the agency's
+# routes match (e.g. Movia including harbour ferries 991/992).
+#
+# ORDER MATTERS: the styles listed here are also the priority order when a
+# pole is served by several modes (e.g. S-tog and regional trains share
+# platforms at hub stations — the higher style wins).
 modes:
-  subway:
-    osm:
-      route: subway
-    gtfs_route_type: 1
-    style: metro
+  metro:
+    agencies: [Metroselskabet]
   s_tog:
-    osm:
-      route: light_rail
-      network: "Takst Sjælland"
-    gtfs_route_type: 2
-    style: s_tog
+    agencies: [DSB S-tog]
+    route_types: [109, 2]
   light_rail:
-    osm:
-      route: light_rail
-      network_neq: "Takst Sjælland"
-    gtfs_route_type: 2
-    style: light_rail
+    agencies: [Hovedstadens Letbane]
   regional:
-    osm:
-      route: train
-    gtfs_route_type: 2
-    style: regional
+    agencies: [DSB, Lokaltog A/S, Skånetrafiken, Snälltåget AB, DSB Vores Tog]
+    route_types: [2]
+  ferry:
+    # harbour buses 991/992 (route_type 4); their shapes cross the harbour
+    # water and are intentionally NOT clipped to the land area polygon
+    agencies: [Movia]
+    route_types: [4]
   bus:
-    osm:
-      route: bus
-      network: Movia
-    gtfs_route_type: 3
-    style: bus
+    agencies: [Movia]
+    route_types: [3, 700, 715]

copenhagen/config/styling.yaml

@@ -4,13 +4,15 @@
 
 zorder:
   bus: 1
-  s_tog: 2
-  light_rail: 3
-  regional: 4
-  metro: 5
+  ferry: 2
+  s_tog: 3
+  light_rail: 4
+  regional: 5
+  metro: 6
 
 line_width:           # body width in points (outline = body + outline_width)
   bus: 0.8
+  ferry: 1.4
   s_tog: 2.0
   light_rail: 2.0
   regional: 1.8
@@ -21,14 +23,16 @@ outline_color: "#2b2b2b"
 
 alpha:
   bus: 0.45
+  ferry: 1.0
   s_tog: 0.95
   light_rail: 0.95
   regional: 0.9
   metro: 1.0
 
-# Fallback palette when no GTFS route_color and no OSM colour tag are present.
-# Metro / S-tog values mirror the OSM colour tags (verified). Regional trains
-# have no OSM colour, so the palette supplies one.
+# Line colours. The Rejseplanen GTFS feed leaves route_color empty for the
+# Copenhagen operators, so in practice everything resolves to this palette.
+# Metro / S-tog values mirror the official line colours; regional trains and
+# the light rail have no operator colour, so the palette supplies one.
 palette:
   metro:
     M1: "#008d41"
@@ -52,6 +56,8 @@ palette:
     C: "#16a085"
     S: "#2a6fb5"
     default: "#c8a80e"
+  ferry:
+    default: "#2aa8e0"   # harbour-bus blue
 
 bus_filters:
   # ref-based bus categorisation (regex applied to route ref)
@@ -73,11 +79,8 @@ basemap:
 stations:
   show: true              # draw station markers (symbols without text)
   styles: [metro, s_tog, regional, bus]
-  stop_cluster_m:          # merge same-name stops within this distance (per mode)
-    metro: 200             # covers platform/entrance spread (~170m max)
-    s_tog: 400             # covers platform/entrance spread (~380m max)
-    regional: 500          # covers platform spread (~475m max)
-    bus: 100               # opposite sides of road, wide boulevards
+  # NOTE: stops/stations arrive pre-collapsed from prepare.py (name-keyed,
+  # one point per stop/station) — no clustering configuration lives here.
   marker:
     shape: circle         # circle | square
     fill: white

copenhagen/scripts/_common.py

@@ -6,7 +6,6 @@ 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"

copenhagen/scripts/build_area.py

@@ -3,7 +3,9 @@
 
 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.
+roles), builds polygons, dissolves the union, and patches in any
+config/area.json "waterways" polygons (the harbour is excluded from the
+kommune boundaries but belongs to the City Pass zone in practice).
 
 Outputs (in data/processed):
     area.geojson  (EPSG:4326, human-readable + portable)
@@ -130,6 +132,23 @@ def main():
     if not dissolved.is_valid:
         dissolved = dissolved.buffer(0)
 
+    # patch waterways into the zone: kommune boundaries exclude water, but
+    # the harbour is practically part of the City Pass area (metro tunnels
+    # beneath it, harbour ferries sail it). See area.json -> waterways.
+    for w in area_cfg.get("waterways", []):
+        wp = Polygon(w["ring"])
+        print(f"adding waterway: {w['name']}", flush=True)
+        dissolved = unary_union([dissolved, wp])
+
+    # Small outward buffer (100 m) to close boundary slivers: the three
+    # relation outlines don't abut perfectly, leaving metre-wide cracks that
+    # otherwise fragment lines clipped to the area (bridge nicks included).
+    dissolved = (
+        gpd.GeoSeries([dissolved], crs="EPSG:4326")
+        .to_crs("EPSG:25832").buffer(100)
+        .to_crs("EPSG:4326").iloc[0]
+    )
+
     PROCESSED.mkdir(parents=True, exist_ok=True)
     gdf = gpd.GeoDataFrame(
         {"name": ["City Pass area"]}, geometry=[dissolved], crs="EPSG:4326"

copenhagen/scripts/download_gtfs.py

@@ -1,40 +1,27 @@
 #!/usr/bin/env python3
-"""Download the Rejseplanen GTFS feed for the Copenhagen area.
+"""Download the Rejseplanen GTFS feed (all of Denmark) and unzip into data/raw/gtfs/.
 
-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/.
+Source: https://www.rejseplanen.info/labs — a static GTFS zip published by
+Rejseplanen covering DSB, DSB S-tog, Metroselskabet, Movia, Lokaltog,
+Hovedstadens Letbane, Skånetrafiken, etc. Nationwide feed; gtfs_to_geopackage.py
+filters it down to the Copenhagen area.
 """
-import sys
 import urllib.request
 import zipfile
-from pathlib import Path
 
-from _common import GTFS_RAW
+from _common import GTFS_RAW, UA
 
-# 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"
+GTFS_URL = "https://www.rejseplanen.info/labs/GTFS.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)
+    req = urllib.request.Request(GTFS_URL, headers={"User-Agent": UA})
+    with urllib.request.urlopen(req) as r, zip_path.open("wb") as f:
+        while chunk := r.read(1 << 20):
+            f.write(chunk)
     with zipfile.ZipFile(zip_path) as z:
         z.extractall(GTFS_RAW)
     zip_path.unlink(missing_ok=True)

copenhagen/scripts/download_osm.py

@@ -1,302 +0,0 @@
-#!/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 platform node member, carrying ref/name/mode/route_refs
-(route refs of routes using this stop, for disambiguation).
-"""
-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", "stop_entry_only", "stop_exit_only"}
-# bus stops: use platform nodes (where passengers wait)
-# rail stops: use stop nodes (platforms are often areas, stop nodes are reliable points)
-BUS_STOP_ROLES = PLATFORM_ROLES
-RAIL_STOP_ROLES = STOP_ROLES | PLATFORM_ROLES
-
-
-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):
-    is_bus = mode == "bus"
-    stop_roles = BUS_STOP_ROLES if is_bus else RAIL_STOP_ROLES
-
-    nodes = {}
-    for e in data["elements"]:
-        if e["type"] == "node":
-            nodes[e["id"]] = e
-
-    lines = []
-    # collect route_refs per stop node: {node_id: set(route_refs)}
-    node_routes = {}
-    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
-        route_ref = tags.get("ref")
-        # collect line geometries
-        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": route_ref,
-            "name": tags.get("name"),
-            "network": tags.get("network"),
-            "colour": tags.get("colour"),
-            "operator": tags.get("operator"),
-            "route": tags.get("route"),
-            "geometry": geom,
-        })
-        # record route_ref on each stop member node
-        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:
-                continue
-            node_routes.setdefault(nid, set())
-            if route_ref:
-                node_routes[nid].add(route_ref)
-
-    # build stops from collected nodes
-    stops = []
-    stop_seen = set()
-    for nid, routes in node_routes.items():
-        if 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"),
-            "route_refs": ";".join(sorted(routes)),
-            "geometry": Point(nd["lon"], nd["lat"]),
-        })
-    return lines, stops
-
-
-def stations_query(bbox):
-    s, n, w, e = bbox
-    return f"""[out:json][timeout:180];
-(
-  node["railway"="station"]({s},{w},{n},{e});
-  node["public_transport"="station"]({s},{w},{n},{e});
-  way["railway"="station"]({s},{w},{n},{e});
-  way["public_transport"="station"]({s},{w},{n},{e});
-);
-out center;
-"""
-
-
-def parse_stations(data):
-    """Stations (nodes or area centroids) with proper human-readable names."""
-    rows = []
-    seen = set()
-    for e in data["elements"]:
-        if e["type"] not in ("node", "way"):
-            continue
-        tags = e.get("tags", {})
-        name = tags.get("name")
-        if not name:
-            continue
-        if e["type"] == "node":
-            geom = Point(e["lon"], e["lat"])
-        else:
-            c = e.get("center")
-            if not c:
-                continue
-            geom = Point(c["lon"], c["lat"])
-        key = (name, round(geom.x, 5), round(geom.y, 5))
-        if key in seen:
-            continue
-        seen.add(key)
-        rows.append({
-            "name": name,
-            "railway": tags.get("railway"),
-            "public_transport": tags.get("public_transport"),
-            "station": tags.get("station"),
-            "subway": tags.get("subway"),
-            "light_rail": tags.get("light_rail"),
-            "train": tags.get("train"),
-            "uic_ref": tags.get("uic_ref"),
-            "geometry": geom,
-        })
-    return rows
-
-
-def orphan_bus_stops_query(bbox):
-    s, n, w, e = bbox
-    return f"""[out:json][timeout:180];
-(
-  node["highway"="bus_stop"]["public_transport"="platform"]({s},{w},{n},{e});
-);
-out body;
-"""
-
-
-def parse_orphan_bus_stops(data, known_node_ids):
-    """Parse standalone bus platform nodes not already captured by routes."""
-    stops = []
-    for e in data["elements"]:
-        if e["type"] != "node" or e["id"] in known_node_ids:
-            continue
-        tags = e.get("tags", {})
-        if tags.get("public_transport") != "platform":
-            continue
-        name = tags.get("name")
-        if not name:
-            continue
-        stops.append({
-            "mode": "bus",
-            "name": name,
-            "ref": tags.get("ref"),
-            "public_transport": tags.get("public_transport"),
-            "railway": tags.get("railway"),
-            "route_refs": "",  # orphan — no route membership
-            "geometry": Point(e["lon"], e["lat"]),
-        })
-    return 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)
-    bbox = (s, n, w, e)
-
-    OSM_RAW.mkdir(parents=True, exist_ok=True)
-    bus_stop_node_ids = set()
-    for mode, cfg in modes.items():
-        ql = build_query(cfg, bbox)
-        (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)
-        # track bus stop node ids to avoid dupes with orphans
-        if mode == "bus":
-            for e in data["elements"]:
-                if e["type"] == "node" and e.get("tags", {}).get("public_transport") == "platform":
-                    bus_stop_node_ids.add(e["id"])
-        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")
-
-    # orphan bus stops (platform nodes not in any route relation)
-    oql = orphan_bus_stops_query(bbox)
-    (OSM_RAW / "orphan_bus_stops.overpassql").write_text(oql)
-    print("[orphan_bus_stops] querying Overpass...", flush=True)
-    odata = overpass(oql)
-    orphans = parse_orphan_bus_stops(odata, bus_stop_node_ids)
-    print(f"  -> {len(orphans)} orphan bus stops (not in route relations)", flush=True)
-    if orphans:
-        # append to bus.gpkg stops layer
-        existing = gpd.read_file(OSM_RAW / "bus.gpkg", layer="stops")
-        combined = gpd.GeoDataFrame(
-            __import__("pandas").concat([existing, gpd.GeoDataFrame(orphans, crs="EPSG:4326")],
-                                        ignore_index=True),
-            crs="EPSG:4326"
-        )
-        # rewrite stops layer
-        import tempfile, shutil
-        tmp = OSM_RAW / "bus_tmp.gpkg"
-        lines_gdf = gpd.read_file(OSM_RAW / "bus.gpkg", layer="lines")
-        lines_gdf.to_file(tmp, driver="GPKG", layer="lines")
-        combined.to_file(tmp, driver="GPKG", layer="stops")
-        shutil.move(str(tmp), str(OSM_RAW / "bus.gpkg"))
-        print(f"  -> bus stops layer now: {len(combined)} total", flush=True)
-
-    # stations (named station nodes/areas) for labelling
-    sql = stations_query(bbox)
-    (OSM_RAW / "stations.overpassql").write_text(sql)
-    print("[stations] querying Overpass...", flush=True)
-    sdata = overpass(sql)
-    st = parse_stations(sdata)
-    print(f"  -> {len(st)} named stations", flush=True)
-    if st:
-        gpd.GeoDataFrame(st, crs="EPSG:4326").to_file(
-            OSM_RAW / "stations.gpkg", driver="GPKG", layer="stations"
-        )
-    print("done.")
-
-
-if __name__ == "__main__":
-    main()

copenhagen/scripts/gtfs_to_geopackage.py

@@ -1,41 +1,150 @@
 #!/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)
+Inputs (data/raw/gtfs): routes.txt, trips.txt, stops.txt, shapes.txt
+Outputs (data/processed):
+    gtfs_shapes.gpkg   (shapes.txt -> LineString per shape_id, EPSG:4326,
+                        each shape tagged with its route_id)
     gtfs_stops.gpkg    (stops.txt -> Point per stop, EPSG:4326)
-    route_colors.csv   (route_id, route_short_name, route_type, route_color)
+    route_colors.csv   (route_id, agency_id, route_short_name, route_long_name,
+                        route_type, route_color, route_text_color)
 
-Uses partridge for fast GTFS parsing. prepare.py reads these when present and
-falls back to OSM-only data when absent.
+The feed is nationwide, so everything is pre-filtered to the Copenhagen area:
+stops by location, shapes to those intersecting the area polygon (or a
+fallback bbox when data/processed/area.gpkg does not exist yet), and routes
+to those having at least one surviving shape. prepare.py clips precisely to
+the area later.
+
+Note: this feed leaves route_color empty for the Copenhagen operators, so
+prepare.py will mostly fall through to OSM colour tags / the styling palette.
+stop_times.txt is not needed here and is intentionally not parsed (220 MB).
 """
 import sys
-from pathlib import Path
+
+import geopandas as gpd
+import pandas as pd
+from shapely.geometry import LineString, box
+from shapely.geometry.base import BaseGeometry
 
 from _common import GTFS_RAW, PROCESSED
 
+# Generous Copenhagen bbox (covers København, Frederiksberg, Amager, and
+# immediate surroundings: airport, Hellerup, Lyngby, Brøndby, ...).
+CPH_BBOX = box(12.30, 55.55, 12.75, 55.85)
+
+
+def load_area() -> BaseGeometry:
+    """City Pass area polygon if build_area has run, else the fallback bbox."""
+    area_path = PROCESSED / "area.gpkg"
+    if area_path.exists():
+        return gpd.read_file(area_path).geometry.union_all()
+    print(f"no {area_path}; using fallback Copenhagen bbox", flush=True)
+    return CPH_BBOX
+
+
+def build_stops(area: BaseGeometry) -> gpd.GeoDataFrame:
+    stops = pd.read_csv(
+        GTFS_RAW / "stops.txt",
+        usecols=["stop_id", "stop_name", "stop_lat", "stop_lon",
+                 "location_type", "parent_station"],
+        dtype={"stop_id": str, "parent_station": str},
+    )
+    minx, miny, maxx, maxy = area.bounds
+    in_bbox = (
+        stops["stop_lat"].between(miny, maxy)
+        & stops["stop_lon"].between(minx, maxx)
+    )
+    stops = stops[in_bbox].copy()
+    g = gpd.GeoDataFrame(
+        stops,
+        geometry=gpd.points_from_xy(stops["stop_lon"], stops["stop_lat"]),
+        crs="EPSG:4326",
+    )
+    return g[g.intersects(area)].drop(columns=["stop_lat", "stop_lon"])
+
+
+def build_shapes(area: BaseGeometry) -> gpd.GeoDataFrame:
+    pts = pd.read_csv(GTFS_RAW / "shapes.txt", dtype={"shape_id": str})
+    minx, miny, maxx, maxy = area.bounds
+    hits = pts[
+        pts["shape_pt_lat"].between(miny, maxy)
+        & pts["shape_pt_lon"].between(minx, maxx)
+    ]
+    pts = pts[pts["shape_id"].isin(set(hits["shape_id"].unique()))]
+    pts = pts.sort_values(["shape_id", "shape_pt_sequence"])
+    lines = (
+        pts.groupby("shape_id")
+        .apply(
+            lambda g: LineString(zip(g["shape_pt_lon"], g["shape_pt_lat"])),
+            include_groups=False,
+        )
+        .rename("geometry")
+        .reset_index()
+    )
+    g = gpd.GeoDataFrame(lines, geometry="geometry", crs="EPSG:4326")
+    return g[g.intersects(area)]
+
+
+def build_routes(
+    shapes: gpd.GeoDataFrame,
+) -> tuple[gpd.GeoDataFrame, pd.DataFrame]:
+    """Tag shapes with route_id and keep routes having >= 1 surviving shape."""
+    trips = pd.read_csv(
+        GTFS_RAW / "trips.txt",
+        usecols=["route_id", "shape_id"],
+        dtype=str,
+    ).dropna(subset=["shape_id"])
+    shape_to_route = (
+        trips[trips["shape_id"].isin(set(shapes["shape_id"]))]
+        .drop_duplicates("shape_id")
+        .set_index("shape_id")["route_id"]
+    )
+    shapes = shapes.copy()
+    shapes["route_id"] = shapes["shape_id"].map(shape_to_route)
+
+    routes = pd.read_csv(GTFS_RAW / "routes.txt", dtype=str).fillna("")
+    routes = routes[
+        routes["route_id"].isin(set(shape_to_route.unique()))
+    ].copy()
+    return shapes, routes[
+        [
+            "route_id", "agency_id", "route_short_name", "route_long_name",
+            "route_type", "route_color", "route_text_color",
+        ]
+    ]
+
 
 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.",
+            "gtfs_to_geopackage: no GTFS data found in "
+            f"{GTFS_RAW}. Run download_gtfs.py first.",
             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)
+    PROCESSED.mkdir(parents=True, exist_ok=True)
+    area = load_area()
+
+    stops = build_stops(area)
+    print(f"stops in area: {len(stops)}", flush=True)
+
+    shapes = build_shapes(area)
+    print(f"shapes in area: {len(shapes)}", flush=True)
+
+    shapes, routes = build_routes(shapes)
+    print(f"routes in area: {len(routes)}", flush=True)
+
+    stops.to_file(PROCESSED / "gtfs_stops.gpkg", driver="GPKG", layer="stops")
+    shapes.to_file(
+        PROCESSED / "gtfs_shapes.gpkg", driver="GPKG", layer="shapes"
+    )
+    routes.to_csv(PROCESSED / "route_colors.csv", index=False)
+    print(
+        "wrote gtfs_stops.gpkg, gtfs_shapes.gpkg, route_colors.csv "
+        f"into {PROCESSED}",
+        flush=True,
+    )
 
 
 if __name__ == "__main__":

copenhagen/scripts/prepare.py

@@ -1,25 +1,44 @@
 #!/usr/bin/env python3
-"""Merge OSM layers into a single master GeoPackage, enriched for rendering.
+"""Merge GTFS-derived layers into a single master GeoPackage, enriched for rendering.
 
-Reads data/raw/osm/{mode}.gpkg (lines + stops) and data/processed/area.gpkg,
-and writes data/processed/master.gpkg with layers "lines" and "stops".
+Reads (from data/processed):
+    gtfs_shapes.gpkg   route-tagged line geometries (gtfs_to_geopackage.py)
+    gtfs_stops.gpkg    raw stop poles (one row per physical pole)
+    route_colors.csv   route metadata (agency_id, short/long name, type, colour)
+    area.gpkg          City Pass area polygon (build_area.py, incl. waterways)
+and from data/raw/gtfs: agency.txt, trips.txt, stop_times.txt.
 
-Enrichment:
-  - style      : from config/modes.yaml (mode -> style key)
-  - colour     : GTFS route_color (if present) -> OSM colour tag -> palette
-  - bus_category : A / C / S / regular (regex from styling.yaml)
-Lines and stops are clipped to the City Pass area polygon.
+Writes data/processed/master.gpkg with layers:
+    lines     one feature per (style, ref, direction) — the shape that serves
+              the most in-area stops (so drawn lines pass the drawn stops)
+    stops     bus/ferry stops, one point per stop name, pruned to stops
+              actually served by a drawn line (see STOP_LINE_MARGIN_M)
+    stations  rail-family stations, one point per (name, style)
+
+Stop names in the Rejseplanen feed are unambiguous per location (verified:
+max spread within a name is ~350 m), so every merge is key-based on names —
+no distance-based clustering anywhere.
+
+Style classification comes from config/modes.yaml (agency + route_type);
+file order is the priority when a pole is served by several modes.
+Colours: GTFS route_color -> styling.yaml palette (the feed leaves
+route_color empty for the Copenhagen operators, so the palette wins).
 """
 import re
-from pathlib import Path
+import sys
 
 import geopandas as gpd
 import pandas as pd
 import yaml
 
-from _common import CONFIG, OSM_RAW, PROCESSED
+from _common import CONFIG, GTFS_RAW, PROCESSED
 
-MIXED_CRS_WARN = "mixed CRS"
+LENGTH_CRS = "EPSG:25832"  # UTM 32N — metres, for length/centroid computations
+STATION_STYLES = {"metro", "s_tog", "light_rail", "regional"}
+# A stop is drawn only if a drawn line of one of its serving refs passes
+# within this distance. Routes have variants; without this check a stop can
+# end up further from the map than the variant we chose not to draw.
+STOP_LINE_MARGIN_M = 300
 
 
 def load_modes():
@@ -30,173 +49,269 @@ def load_styling():
     return yaml.safe_load((CONFIG / "styling.yaml").read_text())
 
 
-def collect_lines(modes):
-    frames = []
-    for mode, cfg in modes.items():
-        p = OSM_RAW / f"{mode}.gpkg"
-        if not p.exists():
-            continue
-        g = gpd.read_file(p, layer="lines")
-        if g.empty:
-            continue
-        g["mode"] = mode
-        g["style"] = cfg["style"]
-        g["gtfs_route_type"] = cfg.get("gtfs_route_type")
-        frames.append(g)
-    if not frames:
-        return gpd.GeoDataFrame(columns=["mode", "ref", "style", "geometry"], crs="EPSG:4326")
-    return pd.concat(frames, ignore_index=True)
+def load_area():
+    area = gpd.read_file(PROCESSED / "area.gpkg")
+    return area.geometry.union_all()
 
 
-def collect_stops(modes):
-    frames = []
-    for mode, cfg in modes.items():
-        p = OSM_RAW / f"{mode}.gpkg"
-        if not p.exists():
+def classify(agency, route_type, modes):
+    """Map a GTFS (agency_name, route_type) pair to a style key, or None."""
+    for style, cfg in modes.items():
+        if agency not in cfg.get("agencies", []):
             continue
-        try:
-            g = gpd.read_file(p, layer="stops")
-        except Exception:
+        types = cfg.get("route_types")
+        if types and route_type not in types:
             continue
-        if g.empty:
-            continue
-        g["mode"] = mode
-        g["style"] = cfg["style"]
-        frames.append(g)
-    if not frames:
-        return None
-    return pd.concat(frames, ignore_index=True)
+        return style
+    return None
 
 
-def gtfs_colour_map():
-    """Return {(route_short_name, route_type): route_color} if GTFS data exists."""
-    csv = PROCESSED / "route_colors.csv"
-    if not csv.exists():
-        return {}
-    df = pd.read_csv(csv)
-    out = {}
-    for _, r in df.iterrows():
-        name = str(r.get("route_short_name") or "").strip()
-        if name:
-            out[(name, r.get("route_type"))] = r.get("route_color")
-    return out
+def load_routes(modes):
+    """route_colors.csv joined with agency names and classified by style."""
+    agencies = pd.read_csv(GTFS_RAW / "agency.txt", dtype=str)
+    agency_names = dict(zip(agencies["agency_id"], agencies["agency_name"]))
+    routes = pd.read_csv(PROCESSED / "route_colors.csv", dtype=str)
+    routes["route_type"] = pd.to_numeric(routes["route_type"], errors="coerce")
+    routes["agency_name"] = routes["agency_id"].map(agency_names)
+    routes["style"] = [
+        classify(a, t, modes)
+        for a, t in zip(routes["agency_name"], routes["route_type"])
+    ]
+    return routes
 
 
 def categorise_bus(ref, patterns):
-    if ref is None:
-        return "regular"
-    s = str(ref)
+    s = "" if ref is None else str(ref)
     for cat, pat in patterns.items():
         if re.search(pat, s):
             return cat
     return "regular"
 
 
-def resolve_colour(row, palette, gtfs):
-    style = row["style"]
-    ref = row.get("ref")
-    colour = row.get("colour")
-    # 1. GTFS (if available)
-    if gtfs and isinstance(ref, str):
-        c = gtfs.get((ref, row.get("gtfs_route_type")))
-        if c:
-            return str(c)
-    # 2. OSM colour tag (only accept real hex/named strings, not NA/None)
-    if isinstance(colour, str) and colour.strip():
-        return colour
-    # 3. palette
+def resolve_colour(style, ref, bus_category, gtfs_colour, palette):
+    # 1. GTFS route_color (mostly empty in this feed)
+    if isinstance(gtfs_colour, str) and gtfs_colour.strip():
+        return gtfs_colour
+    # 2. palette
     p = palette.get(style, {})
     if style == "bus":
-        cat = row.get("bus_category")
-        if not isinstance(cat, str):
-            cat = "regular"
-        return p.get(cat) or p.get("default")
-    key = ref if isinstance(ref, str) else None
-    return p.get(key) or p.get("default")
+        return p.get(bus_category) or p.get("default")
+    return p.get(ref) or p.get("default")
+
+
+def build_lines(routes, styling, area_geom, trips, st, area_stop_ids):
+    """One feature per (style, ref, direction).
+
+    Geometry per group: the shape that serves the most in-area stops
+    (tie-break: longest). Picking by length alone can draw a variant that
+    skips stops shown on the map (terminal stubs, short-turn branches).
+    """
+    shapes = gpd.read_file(PROCESSED / "gtfs_shapes.gpkg")
+
+    shape_rows = trips.dropna(subset=["shape_id"]).drop_duplicates("shape_id")
+    shape_direction = dict(zip(shape_rows["shape_id"], shape_rows["direction_id"]))
+
+    # coverage: distinct in-area stops visited per shape
+    trip_shape = dict(zip(trips["trip_id"], trips["shape_id"]))
+    sv = st.assign(shape_id=st["trip_id"].map(trip_shape))
+    sv = sv[sv["stop_id"].isin(area_stop_ids)]
+    coverage = sv.groupby("shape_id")["stop_id"].nunique()
+
+    use = routes.dropna(subset=["style"])[
+        ["route_id", "route_short_name", "route_long_name",
+         "route_color", "style"]
+    ]
+    g = shapes.merge(use, on="route_id", how="inner")
+    g["direction_id"] = g["shape_id"].map(shape_direction)
+
+    g = g.to_crs(LENGTH_CRS)
+    g["_len"] = g.geometry.length
+    g["_cov"] = g["shape_id"].map(coverage).fillna(0)
+    g = (
+        g.sort_values(["_cov", "_len"], ascending=False)
+        .drop_duplicates(["style", "route_short_name", "direction_id"])
+        .to_crs("EPSG:4326")
+    )
+
+    g["ref"] = g["route_short_name"]
+    g["name"] = [
+        ln if isinstance(ln, str) and ln.strip() else ref
+        for ln, ref in zip(g["route_long_name"], g["ref"])
+    ]
+    patterns = styling["bus_filters"]["categories"]
+    g["bus_category"] = [
+        categorise_bus(ref, patterns) if style == "bus" else None
+        for ref, style in zip(g["ref"], g["style"])
+    ]
+    palette = styling["palette"]
+    g["colour_final"] = [
+        resolve_colour(style, ref, cat, gtfs_c, palette)
+        for style, ref, cat, gtfs_c
+        in zip(g["style"], g["ref"], g["bus_category"], g["route_color"])
+    ]
+
+    g = g[["ref", "name", "style", "bus_category", "colour_final", "geometry"]]
+    g = g[~g.geometry.isna() & g.geometry.is_valid]
+    # the area polygon includes patched-in waterways (area.json), so ferry
+    # shapes and sub-harbour metro tunnels survive the clip unfragmented
+    return gpd.clip(g, area_geom)
+
+
+def build_pole_classes(routes, modes, trips, st, area_stop_ids):
+    """Per-pole classification and serving refs, from stop_times.
+
+    Returns (pole_styles, pole_refs):
+      pole_styles[stop_id] = highest-priority style among classified serving
+                             routes (modes.yaml order)
+      pole_refs[stop_id]   = set of serving route short names (classified)
+    """
+    classified = routes.dropna(subset=["style"])
+    route_style = dict(zip(classified["route_id"], classified["style"]))
+    route_ref = dict(
+        zip(classified["route_id"], classified["route_short_name"])
+    )
+    trip_route = dict(zip(trips["trip_id"], trips["route_id"]))
+
+    sv = st[st["stop_id"].isin(area_stop_ids)].copy()
+    sv["route_id"] = sv["trip_id"].map(trip_route)
+    sv = sv.dropna(subset=["route_id"])
+    sv["style"] = sv["route_id"].map(route_style)
+    sv = sv.dropna(subset=["style"])
+
+    order = {s: i for i, s in enumerate(modes)}
+    sv["_p"] = sv["style"].map(order)
+    best = sv.loc[sv.groupby("stop_id")["_p"].idxmin(), ["stop_id", "style"]]
+    pole_styles = dict(zip(best["stop_id"], best["style"]))
+
+    pole_refs = sv.groupby("stop_id")["route_id"].apply(
+        lambda ids: {route_ref[i] for i in ids}
+    ).to_dict()
+    return pole_styles, pole_refs
+
+
+def build_stops_and_stations(pole_styles, pole_refs, area_geom, lines):
+    poles = gpd.read_file(PROCESSED / "gtfs_stops.gpkg")
+    poles["style"] = poles["stop_id"].map(pole_styles)
+    poles["refs"] = poles["stop_id"].map(lambda s: pole_refs.get(s, set()))
+    poles = poles.dropna(subset=["style"])
+
+    # project for accurate centroids and distances
+    poles = poles.to_crs(LENGTH_CRS)
+
+    # stops (everything not rail-family: buses + harbour ferries), one row
+    # per name, union of serving refs across its poles
+    not_rail = poles[~poles["style"].isin(STATION_STYLES)]
+    stop_rows = []
+    for (name, style), grp in not_rail.groupby(["stop_name", "style"]):
+        refs = set().union(*grp["refs"]) if len(grp) else set()
+        stop_rows.append({
+            "name": name,
+            "style": style,
+            "n_poles": len(grp),
+            "refs": refs,
+            "geometry": grp.geometry.union_all().centroid,
+        })
+
+    # stations: one row per (name, style)
+    rail = poles[poles["style"].isin(STATION_STYLES)]
+    station_rows = []
+    for (name, style), grp in rail.groupby(["stop_name", "style"]):
+        station_rows.append({
+            "name": name,
+            "style": style,
+            "n_poles": len(grp),
+            "geometry": grp.geometry.union_all().centroid,
+        })
+
+    area_25832 = gpd.GeoSeries([area_geom], crs="EPSG:4326").to_crs(LENGTH_CRS).union_all()
+
+    stops = gpd.GeoDataFrame(stop_rows, crs=LENGTH_CRS)
+
+    # prune stops not served by any drawn line within STOP_LINE_MARGIN_M
+    if len(stops) and len(lines):
+        line_geom = (
+            lines.to_crs(LENGTH_CRS)
+            .groupby("ref")["geometry"]
+            .agg(lambda g: g.union_all())
+            .to_dict()
+        )
+        keep = []
+        for _, r in stops.iterrows():
+            dists = [
+                geom.distance(r.geometry)
+                for rf in r["refs"]
+                if (geom := line_geom.get(rf)) is not None
+            ]
+            keep.append(bool(dists) and min(dists) <= STOP_LINE_MARGIN_M)
+        n_dropped = (~pd.Series(keep, index=stops.index)).sum()
+        if n_dropped:
+            print(f"pruned {n_dropped} stops not served by a drawn line "
+                  f"(>{STOP_LINE_MARGIN_M} m from nearest)", flush=True)
+        stops = stops[keep].drop(columns=["refs"])
+    else:
+        stops = stops.drop(columns=["refs"])
+
+    stops = gpd.clip(stops, area_25832)
+    stations = gpd.clip(
+        gpd.GeoDataFrame(station_rows, crs=LENGTH_CRS), area_25832
+    )
+    return stops.to_crs("EPSG:4326"), stations.to_crs("EPSG:4326")
 
 
 def main():
+    for dep in ("gtfs_shapes.gpkg", "gtfs_stops.gpkg", "route_colors.csv"):
+        if not (PROCESSED / dep).exists():
+            print(f"prepare: missing {PROCESSED / dep} — "
+                  "run gtfs_to_geopackage.py first.", file=sys.stderr)
+            sys.exit(2)
+
     modes = load_modes()
     styling = load_styling()
-    palette = styling["palette"]
-    bus_patterns = styling["bus_filters"]["categories"]
+    area_geom = load_area()
 
-    area = gpd.read_file(PROCESSED / "area.gpkg")
-    area_geom = area.geometry.union_all()
+    routes = load_routes(modes)
+    print(f"routes classified: {routes['style'].notna().sum()} of "
+          f"{len(routes)} map to a style", flush=True)
 
-    lines = collect_lines(modes)
-    print(f"collected {len(lines)} raw lines across modes", flush=True)
+    area_stop_ids = set(
+        gpd.read_file(PROCESSED / "gtfs_stops.gpkg")["stop_id"]
+    )
+    trips = pd.read_csv(
+        GTFS_RAW / "trips.txt",
+        usecols=["trip_id", "route_id", "shape_id", "direction_id"],
+        dtype=str,
+    )
+    st = pd.read_csv(
+        GTFS_RAW / "stop_times.txt",
+        usecols=["trip_id", "stop_id"],
+        dtype=str,
+    )
 
-    # bus categorisation (only meaningful for bus style)
-    lines["bus_category"] = [
-        categorise_bus(r, bus_patterns) if s == "bus" else None
-        for r, s in zip(lines.get("ref"), lines["style"])
-    ]
+    lines = build_lines(routes, styling, area_geom, trips, st, area_stop_ids)
+    print(f"lines: {len(lines)} features (one per style/ref/direction)",
+          flush=True)
 
-    gtfs = gtfs_colour_map()
-    if gtfs:
-        print(f"GTFS colours loaded: {len(gtfs)} routes", flush=True)
-    else:
-        print("no GTFS colour data; using OSM colour -> palette", flush=True)
+    pole_styles, pole_refs = build_pole_classes(
+        routes, modes, trips, st, area_stop_ids
+    )
+    stops, stations = build_stops_and_stations(
+        pole_styles, pole_refs, area_geom, lines
+    )
+    print(f"stops: {len(stops)} (bus + ferry); stations: {len(stations)}",
+          flush=True)
 
-    lines["colour_final"] = [
-        resolve_colour(r, palette, gtfs) for _, r in lines.iterrows()
-    ]
-
-    # clip to area
-    lines = lines[~lines.geometry.isna() & lines.geometry.is_valid]
-    clipped = gpd.clip(lines, area_geom)
-    print(f"clipped to area: {len(clipped)} line features remain", flush=True)
-
-    out_cols = [
-        "mode", "style", "ref", "name", "network", "operator", "route",
-        "colour", "bus_category", "colour_final", "geometry",
-    ]
-    for c in out_cols:
-        if c not in clipped.columns:
-            clipped[c] = None
-    clipped = clipped.set_geometry("geometry")
-    clipped = gpd.GeoDataFrame(clipped[out_cols], crs="EPSG:4326")
-
-    PROCESSED.mkdir(parents=True, exist_ok=True)
-    # drop existing master.gpkg so layer overwrite is clean
     out = PROCESSED / "master.gpkg"
     if out.exists():
         out.unlink()
-    clipped.to_file(out, driver="GPKG", layer="lines")
-    print(f"wrote lines layer: {len(clipped)} features", flush=True)
+    lines.to_file(out, driver="GPKG", layer="lines")
+    print(f"wrote lines layer: {len(lines)} features", flush=True)
+    stops.to_file(out, driver="GPKG", layer="stops")
+    print(f"wrote stops layer: {len(stops)} features", flush=True)
+    stations.to_file(out, driver="GPKG", layer="stations")
+    print(f"wrote stations layer: {len(stations)} features", flush=True)
 
-    # stops (stop_positions along routes)
-    stops = collect_stops(modes)
-    if stops is not None and not stops.empty:
-        stops = stops[~stops.geometry.isna() & stops.geometry.is_valid]
-        stops_clipped = gpd.clip(stops, area_geom)
-        keep = ["mode", "style", "name", "ref", "public_transport", "railway", "route_refs", "geometry"]
-        for c in keep:
-            if c not in stops_clipped.columns:
-                stops_clipped[c] = None
-        stops_clipped = gpd.GeoDataFrame(stops_clipped[keep], crs="EPSG:4326")
-        stops_clipped.to_file(out, driver="GPKG", layer="stops")
-        print(f"wrote stops layer: {len(stops_clipped)} features", flush=True)
-    else:
-        print("no stops to write", flush=True)
-
-    # stations (named station nodes/areas) for labelling
-    st_path = OSM_RAW / "stations.gpkg"
-    if st_path.exists():
-        st = gpd.read_file(st_path, layer="stations")
-        if not st.empty:
-            st = st[~st.geometry.isna() & st.geometry.is_valid]
-            st_clipped = gpd.clip(st, area_geom)
-            st_clipped = gpd.GeoDataFrame(st_clipped, crs="EPSG:4326")
-            st_clipped.to_file(out, driver="GPKG", layer="stations")
-            print(f"wrote stations layer: {len(st_clipped)} named stations", flush=True)
-    else:
-        print("no stations file; labelling will be limited", flush=True)
-
-    # summary by style
     print("\nsummary by style:")
-    for style, grp in clipped.groupby("style"):
+    for style, grp in lines.groupby("style"):
         refs = grp["ref"].dropna().unique()
         print(f"  {style:10s}: {len(grp):4d} feats, {len(refs):3d} refs")
 

copenhagen/scripts/render.py

@@ -1,14 +1,19 @@
 #!/usr/bin/env python3
 """Render the Copenhagen transit map to PNG/SVG/PDF.
 
-Reads data/processed/master.gpkg (layers: lines, stations) + area.gpkg +
-config/styling.yaml, applies CLI filters, and composes a printable map:
-  - Carto Positron (no labels) basemap via contextily (EPSG:3857)
+Reads data/processed/master.gpkg (layers: lines, stops, stations — all
+GTFS-derived by prepare.py) + area.gpkg + config/styling.yaml, applies CLI
+filters, and composes a printable map:
+  - Carto/Esri grey basemap via contextily (EPSG:3857)
   - two-tone lines: dark casing + colour body, per route
   - z-order: bus (bottom) -> s_tog -> light_rail -> regional -> metro (top)
   - shapeburst fade mask outside the City Pass area
   - station labels (metro + S-tog) via adjustText
 
+All transit data comes from the Rejseplanen GTFS feed. Stop/station layers
+are pre-collapsed (one point per stop name / station) by prepare.py, so
+rendering is pure plotting — no clustering or merging here.
+
 Usage:
   uv run python render.py                       # default: all modes, PNG, markers only
   uv run python render.py --labels              # add station name labels
@@ -25,6 +30,7 @@ import matplotlib
 matplotlib.use("Agg")
 import matplotlib.pyplot as plt
 import numpy as np
+import pandas as pd
 import yaml
 from affine import Affine
 from adjustText import adjust_text
@@ -81,15 +87,11 @@ def filter_lines(lines, args):
             refs = list(dict.fromkeys(bus["ref"].dropna()))
             keep_refs = set(refs[:args.max_bus_routes])
             bus = bus[bus["ref"].isin(keep_refs)]
-        out = gpd.GeoDataFrame(pd_concat([other, bus]), crs=out.crs)
+        out = gpd.GeoDataFrame(pd.concat([other, bus], ignore_index=True),
+                               crs=out.crs)
     return out
 
 
-def pd_concat(frames):
-    import pandas as pd
-    return pd.concat(frames, ignore_index=True)
-
-
 def add_basemap(ax, styling, zoom=None):
     if zoom is None:
         zoom = styling["basemap"]["zoom"]
@@ -158,60 +160,23 @@ def plot_lines(ax, lines, styling):
                  capstyle="round", joinstyle="round")
 
 
+RAIL_STYLES = ("metro", "s_tog", "light_rail", "regional")
+
+
 def classify_station(r):
-    """Return the style key for a station row, or None."""
-    if r.get("subway") == "yes" or r.get("station") == "subway":
-        return "metro"
-    if r.get("light_rail") == "yes" or r.get("station") == "light_rail":
-        return "s_tog"
-    if r.get("railway") == "station" or r.get("train") == "yes":
-        return "regional"
-    return None
-
-
-DISTANCE_CRS = "EPSG:25832"  # UTM 32N — accurate meters for Copenhagen
-
-
-def cluster_stops(coords, threshold_m):
-    """Cluster (x, y) coordinates in EPSG:3857 within threshold_m (real meters).
-
-    Reprojects to UTM 32N for accurate distance computation, clusters, then
-    returns centroids in EPSG:3857 for plotting.
-    """
-    from pyproj import Transformer
-    from scipy.cluster.hierarchy import fcluster, linkage
-    from scipy.spatial.distance import pdist
-
-    coords = np.asarray(coords)
-    if len(coords) <= 1:
-        return coords.tolist() if len(coords) else []
-
-    # reproject to UTM for accurate distances
-    to_utm = Transformer.from_crs(TARGET_CRS, DISTANCE_CRS, always_xy=True)
-    back = Transformer.from_crs(DISTANCE_CRS, TARGET_CRS, always_xy=True)
-    utm = np.array([to_utm.transform(x, y) for x, y in coords])
-
-    dists = pdist(utm)
-    links = linkage(dists, method="single")
-    labels = fcluster(links, t=threshold_m, criterion="distance")
-
-    centroids = []
-    for label in set(labels):
-        members = utm[labels == label]
-        ux, uy = members.mean(axis=0)
-        cx, cy = back.transform(ux, uy)
-        centroids.append((cx, cy))
-    return centroids
+    """Style key for a station row, or None. Stations are pre-classified."""
+    s = r.get("style")
+    return s if s in RAIL_STYLES else None
 
 
 def plot_stations(ax, stations, styling, active_styles):
+    """Draw station markers. Input is pre-collapsed: one row per (name, style)."""
     if stations is None or stations.empty:
         return
     cfg = styling.get("stations", {})
     if not cfg.get("show", True):
         return
     marker_styles = set(cfg.get("styles", ["metro", "s_tog", "regional"]))
-    cluster_cfg = cfg.get("stop_cluster_m", {})
     mk = cfg.get("marker", {})
     sizes = mk.get("size", {})
     shape = mk.get("shape", "circle")
@@ -220,37 +185,25 @@ def plot_stations(ax, stations, styling, active_styles):
     lw = mk.get("linewidth", 0.8)
     marker = "o" if shape == "circle" else "s"
 
-    # group coordinates by style, then by name for clustering
     pts_by_style = {}
-    names_by_style = {}
     for _, r in stations.iterrows():
         s = classify_station(r)
         if s is None or s not in marker_styles or s not in active_styles:
             continue
-        coord = (r.geometry.x, r.geometry.y)
-        name = r.get("name") if isinstance(r.get("name"), str) else None
-        pts_by_style.setdefault(s, []).append(coord)
-        names_by_style.setdefault(s, {}).setdefault(name, []).append(coord)
+        pts_by_style.setdefault(s, []).append((r.geometry.x, r.geometry.y))
 
     zord = styling["zorder"]
     top_z = max(zord.values()) + 1  # all markers above all lines
     for s in sorted(pts_by_style, key=lambda k: zord.get(k, 0)):
-        threshold = cluster_cfg.get(s, 50) if isinstance(cluster_cfg, dict) else cluster_cfg
-        # cluster same-name stations to one marker; keep unnamed as-is
-        centroids = []
-        for name, coords in names_by_style[s].items():
-            if name and len(coords) > 1:
-                centroids.extend(cluster_stops(coords, threshold))
-            else:
-                centroids.extend(coords)
+        pts = pts_by_style[s]
         sz = sizes.get(s, 3.0)
-        ax.scatter([p[0] for p in centroids], [p[1] for p in centroids],
+        ax.scatter([p[0] for p in pts], [p[1] for p in pts],
                    s=sz ** 2, marker=marker, c=fill, edgecolors=edge,
                    linewidths=lw, zorder=top_z, alpha=1.0)
 
 
 def rail_station_names(stations):
-    """Return the set of names for rail stations (metro/s_tog/regional)."""
+    """Set of names of rail stations (metro/s_tog/light_rail/regional)."""
     if stations is None or stations.empty:
         return set()
     names = set()
@@ -262,109 +215,34 @@ def rail_station_names(stations):
     return names
 
 
-def group_by_route_membership(rows, cluster_m):
-    """Group stop rows by (name, shared route_ref) into distinct stops.
-
-    Two same-name stops that share at least one route_ref are the same stop.
-    Stops on disjoint routes are different stops. Within each group,
-    platform duplicates are merged by distance clustering.
-    Returns a list of (x, y) centroids.
-    """
-    if rows is None or rows.empty:
-        return []
-
-    coords = [(r.geometry.x, r.geometry.y) for _, r in rows.iterrows()]
-    route_sets = []
-    for _, r in rows.iterrows():
-        rr = r.get("route_refs")
-        if isinstance(rr, str) and rr:
-            route_sets.append(set(rr.split(";")))
-        else:
-            route_sets.append(set())
-
-    n = len(coords)
-    # union-find: stops sharing a route are connected
-    parent = list(range(n))
-
-    def find(x):
-        while parent[x] != x:
-            parent[x] = parent[parent[x]]
-            x = parent[x]
-        return x
-
-    def union(a, b):
-        ra, rb = find(a), find(b)
-        if ra != rb:
-            parent[ra] = rb
-
-    for i in range(n):
-        for j in range(i + 1, n):
-            if route_sets[i] & route_sets[j]:
-                union(i, j)
-
-    # group indices by connected component
-    components = {}
-    for i in range(n):
-        root = find(i)
-        components.setdefault(root, []).append(i)
-
-    centroids = []
-    for indices in components.values():
-        comp_coords = [coords[i] for i in indices]
-        comp_routes = [route_sets[i] for i in indices]
-        has_routes = any(rs for rs in comp_routes)
-        if len(comp_coords) == 1:
-            centroids.append(comp_coords[0])
-        elif has_routes:
-            # stops sharing a route are the same stop; merge to centroid
-            centroids.append(tuple(np.mean(comp_coords, axis=0)))
-        else:
-            # no route info (orphans); fall back to distance clustering
-            centroids.extend(cluster_stops(comp_coords, cluster_m))
-    return centroids
-
-
 def plot_bus_stops(ax, stops, styling, active_styles, rail_names=None):
-    """Draw small markers for bus stops from the stops layer.
+    """Draw small markers for bus stops (pre-collapsed: one point per name).
 
-    Bus stops whose name matches a rail station are skipped — the rail
-    station marker represents that stop. Stops are grouped by shared
-    route membership: same-name stops on disjoint routes are distinct.
+    A bus stop whose name exactly matches a rail station is skipped — the
+    station marker represents it.
     """
-    if stops is None or stops.empty or "bus" not in active_styles:
+    if stops is None or stops.empty or not ({"bus", "ferry"} & active_styles):
         return
     cfg = styling.get("stations", {})
     if not cfg.get("show", True) or "bus" not in set(cfg.get("styles", [])):
         return
-    bus = stops[stops["style"] == "bus"]
+    bus = stops[stops["style"].isin(["bus", "ferry"])]
     if bus.empty:
         return
 
-    if rail_names is None:
-        rail_names = set()
-
-    cluster_cfg = cfg.get("stop_cluster_m", {})
-    cluster_m = cluster_cfg.get("bus", 50) if isinstance(cluster_cfg, dict) else cluster_cfg
+    rail_names = rail_names or set()
     mk = cfg.get("marker", {})
-    sizes = mk.get("size", {})
-    sz = sizes.get("bus", 1.5)
+    sz = mk.get("size", {}).get("bus", 1.5)
     fill = mk.get("fill", "white")
     edge = mk.get("edge", "#2b2b2b")
     lw = mk.get("linewidth", 0.8)
-    shape = mk.get("shape", "circle")
-    marker = "o" if shape == "circle" else "s"
+    marker = "o" if mk.get("shape", "circle") == "circle" else "s"
     top_z = max(styling["zorder"].values()) + 1
 
-    # group by name, then by route membership; skip rail station names
-    centroids = []
-    skipped = 0
-    for name, grp in bus.groupby("name"):
-        if name and name in rail_names:
-            skipped += len(grp)
-            continue
-        centroids.extend(group_by_route_membership(grp, cluster_m))
-
-    ax.scatter([p[0] for p in centroids], [p[1] for p in centroids],
+    pts = [(r.geometry.x, r.geometry.y)
+           for _, r in bus.iterrows()
+           if not (isinstance(r.get("name"), str) and r.get("name") in rail_names)]
+    ax.scatter([p[0] for p in pts], [p[1] for p in pts],
                s=sz ** 2, marker=marker, c=fill, edgecolors=edge,
                linewidths=lw, zorder=top_z, alpha=0.8)
 
@@ -379,33 +257,21 @@ def label_stations(ax, stations, styling, active_styles):
     fmin = cfg.get("min_fontsize", 5)
     fmax = cfg.get("max_fontsize", 9)
 
-    def is_metro(r):
-        return classify_station(r) == "metro"
-
-    def is_stog(r):
-        return classify_station(r) == "s_tog"
-
-    pts = []
+    # stations are unique per (name, style); dedupe by name for labelling
+    seen = set()
+    uniq = []
     for _, r in stations.iterrows():
         s = classify_station(r)
         if s not in label_styles:
             continue
         name = r.get("name")
-        if not isinstance(name, str) or not name.strip():
+        if not isinstance(name, str) or not name.strip() or name in seen:
             continue
-        pts.append((r.geometry.x, r.geometry.y, name, s == "metro"))
+        seen.add(name)
+        uniq.append((r.geometry.x, r.geometry.y, name, s == "metro"))
 
-    if not pts:
+    if not uniq:
         return
-    # dedupe by name (keep first location)
-    seen = {}
-    uniq = []
-    for x, y, name, is_m in pts:
-        if name in seen:
-            continue
-        seen[name] = True
-        uniq.append((x, y, name, is_m))
-
     texts = []
     for x, y, name, is_m in uniq:
         fs = fmax if is_m else fmin
@@ -485,7 +351,7 @@ def main():
             fontsize=fig_w * 1.1, fontweight="bold", color="#222",
             bbox=dict(boxstyle="round,pad=0.3", fc="white", ec="none", alpha=0.7))
     ax.text(0.01, 0.01,
-            "Data: © OpenStreetMap contributors (ODbL) · Base: Esri, HERE",
+            "Transit: Rejseplanen GTFS · Area: © OpenStreetMap contributors (ODbL) · Base: Esri, HERE",
             transform=ax.transAxes, ha="left", va="bottom",
             fontsize=max(5, fig_w * 0.5), color="#666")
 

pyproject.toml

@@ -9,7 +9,6 @@ dependencies = [
     "matplotlib>=3.9",
     "contextily>=1.6",
     "pyogrio>=0.10",
-    "partridge>=1.1",
     "requests>=2.31",
     "pyyaml>=6.0",
     "adjusttext>=1.2",

uv.lock

@@ -415,15 +415,6 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/57/b0/0e52c878c53f245edd3a11020f20979b3f490f245af532c7cae3027754b5/idna-3.19-py3-none-any.whl", hash = "sha256:815e7be7a7806d54abb586dc943addc79e8b2ee16915059658cbeff4b1b43bf4", size = 68550, upload-time = "2026-08-18T05:14:22.343Z" },
 ]
 
-[[package]]
-name = "isoweek"
-version = "1.3.3"
-source = { registry = "https://pypi.org/simple" }
-sdist = { url = "https://files.pythonhosted.org/packages/9a/79/68c68f92d8d55b3b790224bdba879a8fe77756a42c0b719e2e695756eb34/isoweek-1.3.3.tar.gz", hash = "sha256:73f3f7bac443e05a3ab45c32a72048b0c4f26d53d81462ec4b142c7581d3ffe8", size = 6847, upload-time = "2017-01-04T15:16:31.072Z" }
-wheels = [
-    { url = "https://files.pythonhosted.org/packages/c2/d4/fe7e2637975c476734fcbf53776e650a29680194eb0dd21dbdc020ca92de/isoweek-1.3.3-py2.py3-none-any.whl", hash = "sha256:d3324c497d97f1534669de225ec877964222e4cc773a4a99063086f7a4e342b6", size = 7138, upload-time = "2017-01-04T15:16:33.969Z" },
-]
-
 [[package]]
 name = "jetlag-maps"
 version = "0.1.0"
@@ -433,7 +424,6 @@ dependencies = [
     { name = "contextily" },
     { name = "geopandas" },
     { name = "matplotlib" },
-    { name = "partridge" },
     { name = "pyogrio" },
     { name = "pyyaml" },
     { name = "requests" },
@@ -448,7 +438,6 @@ requires-dist = [
     { name = "contextily", specifier = ">=1.6" },
     { name = "geopandas", specifier = ">=1.1" },
     { name = "matplotlib", specifier = ">=3.9" },
-    { name = "partridge", specifier = ">=1.1" },
     { name = "pyogrio", specifier = ">=0.10" },
     { name = "pyyaml", specifier = ">=6.0" },
     { name = "requests", specifier = ">=2.31" },
@@ -652,15 +641,6 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/b2/d6/de0cc74f8d36976aeca0dd2e9cbf711882ff8e177495115fd82459afdc4d/mercantile-1.2.1-py3-none-any.whl", hash = "sha256:30f457a73ee88261aab787b7069d85961a5703bb09dc57a170190bc042cd023f", size = 14779, upload-time = "2021-04-21T14:42:39.841Z" },
 ]
 
-[[package]]
-name = "networkx"
-version = "3.6.1"
-source = { registry = "https://pypi.org/simple" }
-sdist = { url = "https://files.pythonhosted.org/packages/6a/51/63fe664f3908c97be9d2e4f1158eb633317598cfa6e1fc14af5383f17512/networkx-3.6.1.tar.gz", hash = "sha256:26b7c357accc0c8cde558ad486283728b65b6a95d85ee1cd66bafab4c8168509", size = 2517025, upload-time = "2025-12-08T17:02:39.908Z" }
-wheels = [
-    { url = "https://files.pythonhosted.org/packages/9e/c9/b2622292ea83fbb4ec318f5b9ab867d0a28ab43c5717bb85b0a5f6b3b0a4/networkx-3.6.1-py3-none-any.whl", hash = "sha256:d47fbf302e7d9cbbb9e2555a0d267983d2aa476bac30e90dfbe5669bd57f3762", size = 2068504, upload-time = "2025-12-08T17:02:38.159Z" },
-]
-
 [[package]]
 name = "numpy"
 version = "2.5.3"
@@ -789,21 +769,6 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/d2/cf/6a51b2c38980e04c279fd2fa908a1b0982064e860444acfca4ec2e2c8359/pandas-3.0.5-cp314-cp314t-win_arm64.whl", hash = "sha256:3c5015fd1730fbf883647e88068176c839c102cea883ba1769a6f4593bfc1f8c", size = 9509776, upload-time = "2026-07-22T22:19:26.694Z" },
 ]
 
-[[package]]
-name = "partridge"
-version = "1.1.2"
-source = { registry = "https://pypi.org/simple" }
-dependencies = [
-    { name = "charset-normalizer" },
-    { name = "isoweek" },
-    { name = "networkx" },
-    { name = "pandas" },
-]
-sdist = { url = "https://files.pythonhosted.org/packages/75/76/4bcf12ceb495f9e1d9861ed9211e348dd15fabb86208be0491e7ab1b3a96/partridge-1.1.2.tar.gz", hash = "sha256:9d9ba044f4123a83c6419dde073c15ac71b63301ed398386d6d3101ac8381e34", size = 27579, upload-time = "2023-12-03T23:05:32.26Z" }
-wheels = [
-    { url = "https://files.pythonhosted.org/packages/5b/7e/78009acdcf798f821468894702d2accc23058ba35f56c7e18e42d6bf01ec/partridge-1.1.2-py2.py3-none-any.whl", hash = "sha256:18e1bf451153f6f8ac74f165bf8662dae39ee1c283125f44746717461b43a55f", size = 16216, upload-time = "2023-12-03T23:05:30.021Z" },
-]
-
 [[package]]
 name = "pillow"
 version = "12.3.0"