Mask out Skyttehøj and everything south of it (lat 55.628)

Add south_bound_lat to modes.yaml and a load_south_bound() helper in prepare.py. When set, the area geometry is intersected with the north half-plane above the cutoff, so all downstream clipping (lines, stops, stations) cascades automatically. Removed: 94 bus stops (including Skyttehøj), 1 station (Vestamager), and line segments south of the boundary (Metro M1 tail, buses 32-36).
Add bus whitelist, route exclusions, crow-fly shape filter, and KML export
2026-09-18 00:32:02 +02:00 · 2026-09-18 00:24:37 +02:00 · 2026-09-17 21:07:40 +02:00
15 changed files with 1149 additions and 750 deletions
@@ -4,10 +4,11 @@ copenhagen/data/processed/**
 !copenhagen/data/raw/.gitkeep
 !copenhagen/data/processed/.gitkeep
-# Generated map outputs (regenerable via render.py)
+# Generated map outputs (regenerable via render.py / export_google_mymaps.py)
 copenhagen/output/*.png
 copenhagen/output/*.svg
 copenhagen/output/*.pdf
 copenhagen/output/*.kml
 # Python
 .venv/
@@ -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]
      ]
    }
  ]
 }
@@ -0,0 +1,251 @@
 <?xml version='1.0' encoding='utf-8'?>
 <kml xmlns="http://www.opengis.net/kml/2.2">
  <Document>
    <name>Coastline (from Hide+Seek Copenhagen City Pass map)</name>
    <StyleMap id="line-01579B-3131-nodesc">
      <Pair>
        <key>normal</key>
        <styleUrl>#line-01579B-3131-nodesc-normal</styleUrl>
      </Pair>
      <Pair>
        <key>highlight</key>
        <styleUrl>#line-01579B-3131-nodesc-highlight</styleUrl>
      </Pair>
    </StyleMap>
    <Style id="line-01579B-3131-nodesc-highlight">
      <LineStyle>
        <color>ff9b5701</color>
        <width>4.6965</width>
      </LineStyle>
      <BalloonStyle>
        <text>&lt;h3&gt;$[name]&lt;/h3&gt;</text>
      </BalloonStyle>
    </Style>
    <Style id="line-01579B-3131-nodesc-normal">
      <LineStyle>
        <color>ff9b5701</color>
        <width>3.131</width>
      </LineStyle>
      <BalloonStyle>
        <text>&lt;h3&gt;$[name]&lt;/h3&gt;</text>
      </BalloonStyle>
    </Style>
    <Folder>
      <name>Coastline</name>
      <Placemark>
        <name>Coastline</name>
        <styleUrl>#line-01579B-3131-nodesc</styleUrl>
        <LineString>
          <tessellate>1</tessellate>
          <coordinates>
            12.5843171,55.7397957,0
            12.5827293,55.7360025,0
            12.5824718,55.7315324,0
            12.5817422,55.731484,0
            12.5813989,55.7304449,0
            12.5821284,55.7290433,0
            12.5866346,55.7285359,0
            12.5878509,55.7274296,0
            12.5895354,55.7273934,0
            12.5903811,55.7271523,0
            12.6023887,55.7256926,0
            12.6158183,55.7206703,0
            12.620303,55.7236674,0
            12.6216977,55.7241629,0
            12.6243907,55.7247369,0
            12.6262897,55.7249967,0
            12.6284247,55.7251719,0
            12.6303022,55.7251779,0
            12.6324909,55.7250873,0
            12.6324695,55.7249665,0
            12.637909,55.7230934,0
            12.6386434,55.7218282,0
            12.6177346,55.7122447,0
            12.6173913,55.7123233,0
            12.6116943,55.709694,0
            12.6113267,55.7091432,0
            12.6115091,55.7090676,0
            12.6128441,55.7077879,0
            12.612709,55.7073944,0
            12.612709,55.7065058,0
            12.61295,55.7056923,0
            12.6132772,55.7049427,0
            12.6128646,55.7039017,0
            12.6146831,55.7036991,0
            12.6153002,55.7030916,0
            12.6146212,55.7020159,0
            12.6161983,55.7016803,0
            12.6164274,55.6978865,0
            12.619534,55.6980419,0
            12.6289389,55.6973909,0
            12.6295397,55.6974605,0
            12.6305757,55.6946383,0
            12.6328654,55.6905452,0
            12.6359244,55.6877876,0
            12.6341066,55.6871202,0
            12.6367352,55.6835607,0
            12.6369351,55.6835927,0
            12.6370405,55.6835331,0
            12.6387568,55.6845909,0
            12.6386816,55.6846722,0
            12.6388157,55.6847159,0
            12.6392421,55.6843356,0
            12.6390959,55.684283,0
            12.6388532,55.6844572,0
            12.637179,55.6834458,0
            12.6372853,55.6833646,0
            12.6370415,55.6832048,0
            12.6391502,55.6804415,0
            12.6397198,55.6804556,0
            12.6407945,55.6788023,0
            12.6446818,55.6791733,0
            12.6481294,55.6760777,0
            12.6477002,55.6758901,0
            12.6482367,55.6754545,0
            12.6482152,55.6750794,0
            12.6465406,55.6743208,0
            12.6448562,55.6719189,0
            12.6458325,55.6714712,0
            12.6459398,55.6711989,0
            12.6457537,55.6708966,0
            12.6449598,55.6705669,0
            12.6441852,55.6693045,0
            12.6439062,55.6691804,0
            12.6432464,55.6691078,0
            12.6420252,55.6671961,0
            12.6419078,55.6671067,0
            12.6416769,55.6670523,0
            12.639605,55.6655071,0
            12.6396426,55.6652772,0
            12.639571,55.6651256,0
            12.6388107,55.6643574,0
            12.6389868,55.6640928,0
            12.6392577,55.6638084,0
            12.6395473,55.6635557,0
            12.6396974,55.6632324,0
            12.6397753,55.6629964,0
            12.6398423,55.6628557,0
            12.6403183,55.6622263,0
            12.6404072,55.6620556,0
            12.6406298,55.6617788,0
            12.6408639,55.661419,0
            12.6409444,55.6613494,0
            12.641049,55.6611951,0
            12.641113,55.6611449,0
            12.6415267,55.6606773,0
            12.641626,55.6606054,0
            12.6417306,55.6604957,0
            12.641803,55.6604564,0
            12.6422784,55.6599927,0
            12.6423696,55.6599458,0
            12.6427853,55.6595902,0
            12.642856,55.6594981,0
            12.6431243,55.659318,0
            12.6436914,55.6588662,0
            12.6438497,55.6587679,0
            12.6444571,55.658322,0
            12.6446935,55.658102,0
            12.6448129,55.6580369,0
            12.6450637,55.6578107,0
            12.6452367,55.657732,0
            12.6454375,55.6575592,0
            12.6456441,55.6574465,0
            12.6457118,55.657381,0
            12.6459499,55.6572384,0
            12.6464378,55.6568159,0
            12.6469394,55.6565511,0
            12.6473741,55.656235,0
            12.6476763,55.6560435,0
            12.6478747,55.6559391,0
            12.6480178,55.6558298,0
            12.6483299,55.6556358,0
            12.6484803,55.6555661,0
            12.648821,55.6554397,0
            12.6492219,55.655188,0
            12.649856,55.6547974,0
            12.6500907,55.6547195,0
            12.6504703,55.654637,0
            12.6520371,55.6546291,0
            12.6520451,55.6545564,0
            12.6505229,55.6545288,0
            12.6500079,55.6543775,0
            12.6497504,55.6541565,0
            12.6493642,55.6540082,0
            12.6488855,55.6534695,0
            12.6485422,55.6528339,0
            12.6482847,55.6525433,0
            12.6480272,55.652053,0
            12.6481384,55.6512885,0
            12.647838,55.6509979,0
            12.6479453,55.6507679,0
            12.6480847,55.6506529,0
            12.6482242,55.6503804,0
            12.6481906,55.649885,0
            12.6482577,55.649549,0
            12.6485414,55.6491014,0
            12.6486309,55.6489064,0
            12.6488052,55.6486401,0
            12.6489045,55.6485477,0
            12.6490386,55.648516,0
            12.6491539,55.6484191,0
            12.6502279,55.6483093,0
            12.6504773,55.6482442,0
            12.6506651,55.648161,0
            12.6507885,55.6480263,0
            12.6508232,55.6478915,0
            12.6507293,55.6476554,0
            12.6506341,55.6475652,0
            12.650433,55.6474623,0
            12.6502318,55.6474138,0
            12.6499394,55.6474002,0
            12.6485447,55.6475591,0
            12.6483221,55.6473064,0
            12.6480517,55.6467941,0
            12.6479123,55.6463279,0
            12.6478908,55.6459889,0
            12.6481054,55.6455166,0
            12.6484007,55.6451585,0
            12.6485601,55.6447979,0
            12.6488874,55.6444679,0
            12.6500944,55.6444286,0
            12.6505128,55.6443408,0
            12.6507434,55.6441894,0
            12.651591,55.6432419,0
            12.6525459,55.6434296,0
            12.6531499,55.6433718,0
            12.6534396,55.643172,0
            12.6549523,55.6415493,0
            12.6548772,55.6413071,0
            12.6533607,55.640507,0
            12.6547246,55.6390978,0
            12.6564626,55.6388799,0
            12.6574282,55.6382682,0
            12.6566614,55.6368416,0
            12.6569832,55.6364147,0
            12.657974,55.6358075,0
            12.6578238,55.6355471,0
            12.6568636,55.6353685,0
            12.6570406,55.6350324,0
            12.6568199,55.6344235,0
            12.6570774,55.6340874,0
            12.657882,55.6340419,0
            12.6596845,55.6338481,0
            12.661562,55.6338179,0
            12.6637722,55.6339874,0
            12.6649201,55.6340419,0
            12.6671088,55.6340359,0
            12.6694657,55.6338899,0
            12.6714935,55.6336659,0
            12.6735105,55.6333267,0
            12.6766638,55.6325319,0
            12.677919,55.6320049,0
            12.6793996,55.6310904,0
            12.6805918,55.6297217,0
            12.6809995,55.6286073,0
            12.6811927,55.6248094,0
          </coordinates>
        </LineString>
      </Placemark>
    </Folder>
  </Document>
 </kml>
@@ -1,40 +1,56 @@
-# Mapping of map modes to OSM route filters, GTFS route types, and style keys.
+# Mapping of map styles to GTFS agency / route_type filters.
 # `style` references a key in styling.yaml.
 #
-# OSM filter keys map directly to Overpass tag filters. A key suffixed with
+# Used by prepare.py to classify GTFS routes (lines layer) and the routes
-# `_neq` emits a negated tag filter (["k"!="v"]); any other key emits an
+# serving each stop pole (stops/stations layers). Style keys are referenced
-# equality filter (["k"="v"]). The bbox is appended by download_osm.py.
+# by styling.yaml (palette, zorder, widths).
 #
-# S-tog are tagged route=light_rail in OSM (NOT route=train); they are scoped
+# `agencies` matches GTFS agency_name from agency.txt.
-# by network="Takst Sjælland". The separate `light_rail` mode catches other
+# `route_types` optionally restricts GTFS route_type values, including the
-# light rail (e.g. Hovedstadens Letbane, opened Aug 2026) via network!=Takst.
+# extended codes used in this feed (109 = S-tog suburban rail, 700/715 =
-# GTFS route_type: 0=tram,1=subway,2=rail,3=bus,4=ferry.
+# 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).
 # Hard exclusions: specific routes to drop entirely, keyed by (agency,
 # short name) — the pairing rule used everywhere else in this feed.
 # 083 (Snälltåget night train, Malmö-Stockholm): its in-area shapes are
 # sparse crow-fly stubs that draw as straight cuts across the map; the
 # long-distance train is not relevant for the city game, so skip it.
 # Bus whitelist: only these buses are included — there are far more city
 # bus routes than the game scale needs. Match on bus category (per
 # styling.yaml bus_filters) or on specific refs. Remove this section to
 # include all buses; with empty lists no buses are included.
 bus_whitelist:
  categories: [A]   # all "A" city buses (1A, 2A, 4A, 6A, 9A, ...)
  refs: ["21", "30", "31", "32", "33", "34", "35", "36", "37", "38", "39", "68", "77"]
 exclude:
  - agency: Snälltåget AB
    ref: "083"
 # South bound: drop every stop/station south of this latitude (inclusive),
 # and clip transit lines at this boundary. Used to cut the map short of
 # Skyttehøj (Amager Landevej) and everything south of it on Amager.
 south_bound_lat: 55.628
 modes:
-  subway:
+  metro:
-    osm:
+    agencies: [Metroselskabet]
      route: subway
    gtfs_route_type: 1
    style: metro
  s_tog:
-    osm:
+    agencies: [DSB S-tog]
-      route: light_rail
+    route_types: [109, 2]
      network: "Takst Sjælland"
    gtfs_route_type: 2
    style: s_tog
  light_rail:
-    osm:
+    agencies: [Hovedstadens Letbane]
      route: light_rail
      network_neq: "Takst Sjælland"
    gtfs_route_type: 2
    style: light_rail
  regional:
-    osm:
+    agencies: [DSB, Lokaltog A/S, Skånetrafiken, Snälltåget AB, DSB Vores Tog]
-      route: train
+    route_types: [2]
-    gtfs_route_type: 2
+  ferry:
-    style: regional
+    # 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:
+    agencies: [Movia]
-      route: bus
+    route_types: [3, 700, 715]
      network: Movia
    gtfs_route_type: 3
    style: bus
@@ -4,13 +4,15 @@
 zorder:
  bus: 1
-  s_tog: 2
+  ferry: 2
-  light_rail: 3
+  s_tog: 3
-  regional: 4
+  light_rail: 4
-  metro: 5
+  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.
+# Line colours. The Rejseplanen GTFS feed leaves route_color empty for the
-# Metro / S-tog values mirror the OSM colour tags (verified). Regional trains
+# Copenhagen operators, so in practice everything resolves to this palette.
-# have no OSM colour, so the palette supplies one.
+# 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)
+  # NOTE: stops/stations arrive pre-collapsed from prepare.py (name-keyed,
-    metro: 200             # covers platform/entrance spread (~170m max)
+  # one point per stop/station) — no clustering configuration lives here.
    s_tog: 400             # covers platform/entrance spread (~380m max)
    regional: 500          # covers platform spread (~475m max)
    bus: 100               # opposite sides of road, wide boulevards
  marker:
    shape: circle         # circle | square
    fill: white
@@ -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"
@@ -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"
@@ -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
+Source: https://www.rejseplanen.info/labs — a static GTFS zip published by
-public GTFS download. This script documents the intended flow and exits
+Rejseplanen covering DSB, DSB S-tog, Metroselskabet, Movia, Lokaltog,
-non-zero with guidance.
+Hovedstadens Letbane, Skånetrafiken, etc. Nationwide feed; gtfs_to_geopackage.py
-
+filters it down to the Copenhagen area.
 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
+from _common import GTFS_RAW, UA
-# Fill in once access is granted. Likely candidates:
+GTFS_URL = "https://www.rejseplanen.info/labs/GTFS.zip"
 #   - 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)
+    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)
@@ -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()
@@ -0,0 +1,257 @@
 #!/usr/bin/env python3
 """Export the master GeoPackage to a Google My Maps-importable KML.
 Reads (from data/processed): master.gpkg (lines, stops, stations), area.gpkg
 Read (from config): coastline.kml — the coastline layer extracted verbatim
 from the reference map
 (https://www.google.com/maps/d/u/0/kml?mid=17T6ZDYbGz72h_eteL-CxdkiaXcgDLEA&forcekml=1)
 Writes: output/copenhagen-mymaps.kml
 My Maps turns each top-level KML Folder into a layer on import. We emit:
    Bus stops / Bus routes
    Train stations / Train routes   (metro + S-tog + light rail + regional)
    Ferry stops / Ferry routes
    City Pass boundary (outline + faint fill)
    Coastline (verbatim from the reference map's layer)
 Colours match the PNG map: routes use colour_final from prepare.py; stop
 icons use styling.yaml palette defaults (train stations match the reference
 map's red pin). Import manually in Google My Maps:
 Create map → Import → upload the .kml.
 """
 import re
 import sys
 import xml.etree.ElementTree as ET
 import geopandas as gpd
 import yaml
 from _common import CONFIG, OUTPUT, PROCESSED
 KML = "http://www.opengis.net/kml/2.2"
 ET.register_namespace("", KML)
 # My Maps stock pin (same as the reference map uses); tinted via IconStyle.
 ICON_HREF = "https://www.gstatic.com/mapspro/images/stock/503-wht-blank_maps.png"
 ICON_SCALE = 1.0
 TRAIN_STATION_COLOR = "#C2185B"  # reference map's station-pin red
 TRAIN_STYLES = ["metro", "s_tog", "light_rail", "regional"]
 # route line widths in pixels, per style
 LINE_WIDTH = {
    "metro": 4.5,
    "s_tog": 4.25,
    "light_rail": 4.25,
    "regional": 3.5,
    "ferry": 3.75,
    "bus": 2.75,
 }
 BOUNDARY_COLOR = "#006064"
 BOUNDARY_LINE_WIDTH = 3.0
 BOUNDARY_FILL_ALPHA = 0x54  # ~33 %
 class Styles:
    """Registry of unique (kind, color, width) -> KML <Style> ids."""
    def __init__(self):
        self._map = {}
    def get(self, kind, color_hex, width=0.0):
        key = (kind, color_hex, width)
        if key not in self._map:
            self._map[key] = f"s{len(self._map):03d}"
        return self._map[key]
    def elements(self):
        """Yield <Style> elements for everything registered so far."""
        by_id = sorted(self._map.items(), key=lambda kv: kv[1])
        for (kind, color_hex, width), sid in by_id:
            st = ET.Element(f"{{{KML}}}Style", id=sid)
            if kind == "icon":
                el = ET.SubElement(st, f"{{{KML}}}IconStyle")
                el.append(_tex("color", kml_color(color_hex)))
                el.append(_tex("scale", str(ICON_SCALE)))
                icon = ET.SubElement(el, f"{{{KML}}}Icon")
                icon.append(_tex("href", ICON_HREF))
            else:
                el = ET.SubElement(st, f"{{{KML}}}LineStyle")
                el.append(_tex("color", kml_color(color_hex)))
                el.append(_tex("width", str(width)))
                if kind == "poly":
                    poly = ET.SubElement(st, f"{{{KML}}}PolyStyle")
                    poly.append(_tex("color", kml_color(color_hex, BOUNDARY_FILL_ALPHA)))
            balloon = ET.SubElement(st, f"{{{KML}}}BalloonStyle")
            balloon.append(_tex("text", "<h3>$[name]</h3>"))
            yield st
 def _tex(tag, value):
    el = ET.Element(f"{{{KML}}}{tag}")
    el.text = value
    return el
 def kml_color(hex_color, alpha=0xFF):
    """#RRGGBB -> KML AABBGGRR."""
    h = hex_color.lstrip("#")
    return f"{alpha:02x}{h[4:6]}{h[2:4]}{h[0:2]}"
 def placemark(folder_el, name, sid, geom_el):
    pm = ET.SubElement(folder_el, f"{{{KML}}}Placemark")
    pm.append(_tex("name", str(name)))
    pm.append(_tex("styleUrl", f"#{sid}"))
    pm.append(geom_el)
 def point_el(geom):
    pt = ET.Element(f"{{{KML}}}Point")
    pt.append(_tex("coordinates", f"{geom.x:.7f},{geom.y:.7f}"))
    return pt
 def linestring_el(geom):
    ls = ET.Element(f"{{{KML}}}LineString")
    ls.append(_tex("tessellate", "1"))
    ls.append(_tex("coordinates", " ".join(f"{x:.7f},{y:.7f}" for x, y in geom.coords)))
    return ls
 def polygon_el(geom):
    pg = ET.Element(f"{{{KML}}}Polygon")
    pg.append(_tex("tessellate", "1"))
    outer = ET.SubElement(pg, f"{{{KML}}}outerBoundaryIs")
    ring = ET.SubElement(outer, f"{{{KML}}}LinearRing")
    ring.append(_tex("coordinates", " ".join(f"{x:.7f},{y:.7f}" for x, y in geom.exterior.coords)))
    for inner in geom.interiors:
        b = ET.SubElement(pg, f"{{{KML}}}innerBoundaryIs")
        ring = ET.SubElement(b, f"{{{KML}}}LinearRing")
        ring.append(_tex("coordinates", " ".join(f"{x:.7f},{y:.7f}" for x, y in inner.coords)))
    return pg
 def multi_el(geom, single):
    """Wrap (Multi)Geometry into a KML element (MultiGeometry if needed)."""
    geoms = list(getattr(geom, "geoms", [geom]))
    if len(geoms) == 1:
        return single(geoms[0])
    mg = ET.Element(f"{{{KML}}}MultiGeometry")
    for g in geoms:
        mg.append(single(g))
    return mg
 def add_folder(document, name):
    f = ET.SubElement(document, f"{{{KML}}}Folder")
    f.append(_tex("name", name))
    return f
 def add_stops_folder(document, name, gdf, tint, styles):
    f = add_folder(document, name)
    sid = styles.get("icon", tint)
    for _, row in gdf.sort_values("name").iterrows():
        placemark(f, row["name"], sid, point_el(row.geometry))
    return len(gdf)
 def add_routes_folder(document, name, gdf, styles):
    f = add_folder(document, name)
    for _, row in gdf.iterrows():
        sid = styles.get("line", row["colour_final"], LINE_WIDTH[row["style"]])
        placemark(f, f'{row["ref"]} {row["name"]}', sid, multi_el(row.geometry, linestring_el))
    return len(gdf)
 def add_boundary_folder(document, area_geom, styles):
    f = add_folder(document, "City Pass boundary")
    sid = styles.get("poly", BOUNDARY_COLOR, BOUNDARY_LINE_WIDTH)
    placemark(f, "City Pass boundary", sid, multi_el(area_geom, polygon_el))
    return 1
 def append_coastline(document):
    """Append the verbatim coastline layer stored in config/coastline.kml."""
    src = ET.parse(CONFIG / "coastline.kml").getroot().find(f"{{{KML}}}Document")
    n = 0
    for el in src:
        if el.tag == f"{{{KML}}}Folder":
            document.append(el)
            n = len(el.findall(f"{{{KML}}}Placemark"))
        elif el.tag in (f"{{{KML}}}Style", f"{{{KML}}}StyleMap"):
            document.append(el)
    return n
 def natural_key(ref):
    return tuple(int(t) if t.isdigit() else t for t in re.split(r"(\d+)", str(ref)))
 def main():
    lines = gpd.read_file(PROCESSED / "master.gpkg", layer="lines")
    stops = gpd.read_file(PROCESSED / "master.gpkg", layer="stops")
    stations = gpd.read_file(PROCESSED / "master.gpkg", layer="stations")
    area = gpd.read_file(PROCESSED / "area.gpkg").geometry.union_all()
    styling = yaml.safe_load((CONFIG / "styling.yaml").read_text())
    palette = styling["palette"]
    kml = ET.Element(f"{{{KML}}}kml")
    document = ET.SubElement(kml, f"{{{KML}}}Document")
    document.append(_tex("name", "Kurragömma Köpenhamn"))
    styles = Styles()
    counts = {}
    counts["Bus stops"] = add_stops_folder(
        document, "Bus stops", stops[stops["style"] == "bus"],
        palette["bus"]["default"], styles,
    )
    bus_lines = lines[lines["style"] == "bus"].copy()
    bus_lines["_k"] = bus_lines["ref"].map(natural_key)
    counts["Bus routes"] = add_routes_folder(
        document, "Bus routes", bus_lines.sort_values("_k"), styles
    )
    counts["Train stations"] = add_stops_folder(
        document, "Train stations", stations, TRAIN_STATION_COLOR, styles
    )
    train = lines[lines["style"].isin(TRAIN_STYLES)].copy()
    train["_order"] = train["style"].map(TRAIN_STYLES.index)
    train["_k"] = train["ref"].map(natural_key)
    counts["Train routes"] = add_routes_folder(
        document, "Train routes", train.sort_values(["_order", "_k"]), styles
    )
    counts["Ferry stops"] = add_stops_folder(
        document, "Ferry stops", stops[stops["style"] == "ferry"],
        palette["ferry"]["default"], styles,
    )
    ferry_lines = lines[lines["style"] == "ferry"].sort_values("ref")
    counts["Ferry routes"] = add_routes_folder(document, "Ferry routes", ferry_lines, styles)
    counts["City Pass boundary"] = add_boundary_folder(document, area, styles)
    counts["Coastline"] = append_coastline(document)
    # <Style> elements belong before the first <Folder> (KML resolves style
    # ids regardless of order; this is conventional + nicer on the eyes).
    insert_at = list(document).index(document.find(f"{{{KML}}}name")) + 1
    for offset, st in enumerate(styles.elements()):
        document.insert(insert_at + offset, st)
    out = OUTPUT / "copenhagen-mymaps.kml"
    tree = ET.ElementTree(kml)
    ET.indent(tree, space="  ")
    tree.write(out, encoding="utf-8", xml_declaration=True)
    print(f"wrote {out}")
    for n, c in counts.items():
        print(f"  {n}: {c} placemarks")
 if __name__ == "__main__":
    main()
@@ -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).
+Inputs (data/raw/gtfs): routes.txt, trips.txt, stops.txt, shapes.txt
-
+Outputs (data/processed):
-Intended outputs (data/processed):
+    gtfs_shapes.gpkg   (shapes.txt -> LineString per shape_id, EPSG:4326,
-    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
+The feed is nationwide, so everything is pre-filtered to the Copenhagen area:
-falls back to OSM-only data when absent.
+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 "
+            "gtfs_to_geopackage: no GTFS data found in "
-            f"{GTFS_RAW}. Run download_gtfs.py first once a feed URL is set.",
+            f"{GTFS_RAW}. Run download_gtfs.py first.",
            file=sys.stderr,
        )
        sys.exit(2)
-    # TODO (when GTFS available):
+    PROCESSED.mkdir(parents=True, exist_ok=True)
-    #   import partridge as pt
+    area = load_area()
-    #   import geopandas as gpd
+
-    #   from shapely.geometry import LineString, Point
+    stops = build_stops(area)
-    #   feed = pt.load_geo_feed(str(GTFS_RAW))
+    print(f"stops in area: {len(stops)}", flush=True)
-    #   shapes -> gtfs_shapes.gpkg (LineString per shape_id)
+
-    #   stops  -> gtfs_stops.gpkg
+    shapes = build_shapes(area)
-    #   routes -> route_colors.csv (route_id, route_short_name, route_type, route_color)
+    print(f"shapes in area: {len(shapes)}", flush=True)
-    print("gtfs_to_geopackage: not yet implemented (GTFS feed unavailable).")
+
-    sys.exit(2)
+    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__":
@@ -1,202 +1,414 @@
 #!/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,
+Reads (from data/processed):
-and writes data/processed/master.gpkg with layers "lines" and "stops".
+    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:
+Writes data/processed/master.gpkg with layers:
-  - style      : from config/modes.yaml (mode -> style key)
+    lines     one feature per (style, ref, direction) — the shape that serves
-  - colour     : GTFS route_color (if present) -> OSM colour tag -> palette
+              the most in-area stops (so drawn lines pass the drawn stops)
-  - bus_category : A / C / S / regular (regex from styling.yaml)
+    stops     bus/ferry stops, one point per stop name, pruned to stops
-Lines and stops are clipped to the City Pass area polygon.
+              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 shapely.geometry import box
-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
 # Shapes sparser than this (points per km) are "crow-fly" placeholders from
 # the feed (a handful of points for a several-hundred-km line) — they render
 # as straight cuts across the map. Real rail/bus geometry is >= 0.5 pts/km.
 MIN_SHAPE_POINTS_PER_KM = 0.05
 def load_modes():
-    return yaml.safe_load((CONFIG / "modes.yaml").read_text())["modes"]
+    return _load_modes_yaml().get("modes", {})
 def load_exclude():
    """Hard route blacklist from modes.yaml: list of (agency, ref)."""
    return [(e["agency"], str(e["ref"]))
            for e in _load_modes_yaml().get("exclude", [])]
 def load_bus_whitelist():
    """Bus whitelist from modes.yaml: dict with categories/refs sets.
    Returns None when the section is absent (= include all buses).
    """
    wl = _load_modes_yaml().get("bus_whitelist")
    if wl is None:
        return None
    return {
        "categories": set(wl.get("categories") or []),
        "refs": {str(r) for r in (wl.get("refs") or [])},
    }
 def load_south_bound():
    """South latitude cutoff from modes.yaml, or None.
    When set, stops/stations south of (and transit lines below) this
    latitude are dropped. Used to trim the map at a boundary.
    """
    return _load_modes_yaml().get("south_bound_lat")
 def _load_modes_yaml():
    return yaml.safe_load((CONFIG / "modes.yaml").read_text())
 def load_styling():
    return yaml.safe_load((CONFIG / "styling.yaml").read_text())
-def collect_lines(modes):
+def load_area():
-    frames = []
+    area = gpd.read_file(PROCESSED / "area.gpkg")
-    for mode, cfg in modes.items():
+    return area.geometry.union_all()
        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 collect_stops(modes):
+def classify(agency, route_type, modes):
-    frames = []
+    """Map a GTFS (agency_name, route_type) pair to a style key, or None."""
-    for mode, cfg in modes.items():
+    for style, cfg in modes.items():
-        p = OSM_RAW / f"{mode}.gpkg"
+        if agency not in cfg.get("agencies", []):
        if not p.exists():
            continue
-        try:
+        types = cfg.get("route_types")
-            g = gpd.read_file(p, layer="stops")
+        if types and route_type not in types:
        except Exception:
            continue
-        if g.empty:
+        return style
            continue
        g["mode"] = mode
        g["style"] = cfg["style"]
        frames.append(g)
    if not frames:
    return None
    return pd.concat(frames, ignore_index=True)
-def gtfs_colour_map():
+def load_routes(modes, exclude=()):
-    """Return {(route_short_name, route_type): route_color} if GTFS data exists."""
+    """route_colors.csv joined with agency names and classified by style."""
-    csv = PROCESSED / "route_colors.csv"
+    agencies = pd.read_csv(GTFS_RAW / "agency.txt", dtype=str)
-    if not csv.exists():
+    agency_names = dict(zip(agencies["agency_id"], agencies["agency_name"]))
-        return {}
+    routes = pd.read_csv(PROCESSED / "route_colors.csv", dtype=str)
-    df = pd.read_csv(csv)
+    routes["route_type"] = pd.to_numeric(routes["route_type"], errors="coerce")
-    out = {}
+    routes["agency_name"] = routes["agency_id"].map(agency_names)
-    for _, r in df.iterrows():
+    routes["style"] = [
-        name = str(r.get("route_short_name") or "").strip()
+        classify(a, t, modes)
-        if name:
+        for a, t in zip(routes["agency_name"], routes["route_type"])
-            out[(name, r.get("route_type"))] = r.get("route_color")
+    ]
-    return out
+    if exclude:
        mask = [
            (a, str(r)) in set(exclude)
            for a, r in zip(routes["agency_name"], routes["route_short_name"])
        ]
        n = sum(mask)
        if n:
            routes.loc[mask, "style"] = None
            print(f"excluded by modes.yaml: {n} route(s) "
                  f"({sorted(set(zip(routes.loc[mask, 'agency_name'], routes.loc[mask, 'route_short_name'])) )})",
                  flush=True)
    return routes
 def categorise_bus(ref, patterns):
-    if ref is None:
+    s = "" if ref is None else str(ref)
        return "regular"
    s = str(ref)
    for cat, pat in patterns.items():
        if re.search(pat, s):
            return cat
    return "regular"
-def resolve_colour(row, palette, gtfs):
+def resolve_colour(style, ref, bus_category, gtfs_colour, palette):
-    style = row["style"]
+    # 1. GTFS route_color (mostly empty in this feed)
-    ref = row.get("ref")
+    if isinstance(gtfs_colour, str) and gtfs_colour.strip():
-    colour = row.get("colour")
+        return gtfs_colour
-    # 1. GTFS (if available)
+    # 2. palette
    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
    p = palette.get(style, {})
    if style == "bus":
-        cat = row.get("bus_category")
+        return p.get(bus_category) or p.get("default")
-        if not isinstance(cat, str):
+    return p.get(ref) or p.get("default")
-            cat = "regular"
+
-        return p.get(cat) or p.get("default")
+
-    key = ref if isinstance(ref, str) else None
+def shape_points_per_km(shapes_25832):
-    return p.get(key) or p.get("default")
+    """Point density per shape. Geometries are projected (metres)."""
    def n_pts(geom):
        geoms = getattr(geom, "geoms", [geom])
        return sum(len(g.coords) for g in geoms)
    pts = shapes_25832.geometry.map(n_pts)
    return pts / (shapes_25832.geometry.length / 1000.0).clip(lower=1e-6)
 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")
    # The feed contains a few "crow-fly" placeholder shapes for long-distance
    # trains (e.g. Snälltåget, 6-8 points for ~700 km). Those draw as
    # straight lines across the map and can shadow the proper rail-geometry
    # shape in coverage comparison. Drop below a point-density floor;
    # real shapes are >= 0.5 pts/km, placeholders are ~0.01 pts/km, and even
    # the 0.44 km ferry 993 (a handful of points over 440 m) stays well above.
    shapes = shapes.to_crs(LENGTH_CRS)
    pts_km = shape_points_per_km(shapes)
    degenerate = pts_km < MIN_SHAPE_POINTS_PER_KM
    if degenerate.any():
        print(f"dropping {int(degenerate.sum())} degenerate (crow-fly) "
              f"shapes (<{MIN_SHAPE_POINTS_PER_KM:g} pts/km)", flush=True)
    shapes = shapes[~degenerate].to_crs("EPSG:4326")
    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"]
+    area_geom = load_area()
    bus_patterns = styling["bus_filters"]["categories"]
-    area = gpd.read_file(PROCESSED / "area.gpkg")
+    south_bound = load_south_bound()
-    area_geom = area.geometry.union_all()
+    if south_bound is not None:
        north = box(-180, south_bound, 180, 90)
        area_geom = area_geom.intersection(north)
        print(f"south bound: clipping area at lat {south_bound}", flush=True)
-    lines = collect_lines(modes)
+    routes = load_routes(modes, load_exclude())
    print(f"collected {len(lines)} raw lines across modes", flush=True)
-    # bus categorisation (only meaningful for bus style)
+    bus_wl = load_bus_whitelist()
-    lines["bus_category"] = [
+    if bus_wl is not None:
-        categorise_bus(r, bus_patterns) if s == "bus" else None
+        patterns = styling["bus_filters"]["categories"]
-        for r, s in zip(lines.get("ref"), lines["style"])
+        is_bus = routes["style"] == "bus"
-    ]
+        keep = routes["route_short_name"].map(
            lambda r: str(r) in bus_wl["refs"]
                      or categorise_bus(r, patterns) in bus_wl["categories"]
        )
        routes.loc[is_bus & ~keep, "style"] = None
        kept = sorted(
            routes.loc[is_bus & keep, "route_short_name"].unique(),
            key=lambda s: [int(t) if t.isdigit() else t
                           for t in re.split(r"(\d+)", str(s))],
        )
        print(f"bus whitelist: kept {len(kept)} of "
              f"{int(is_bus.sum())} bus refs: {', '.join(kept)}",
              flush=True)
-    gtfs = gtfs_colour_map()
+    print(f"routes classified: {routes['style'].notna().sum()} of "
-    if gtfs:
+          f"{len(routes)} map to a style", flush=True)
        print(f"GTFS colours loaded: {len(gtfs)} routes", flush=True)
    else:
        print("no GTFS colour data; using OSM colour -> palette", flush=True)
-    lines["colour_final"] = [
+    area_stop_ids = set(
-        resolve_colour(r, palette, gtfs) for _, r in lines.iterrows()
+        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,
    )
-    # clip to area
+    lines = build_lines(routes, styling, area_geom, trips, st, area_stop_ids)
-    lines = lines[~lines.geometry.isna() & lines.geometry.is_valid]
+    print(f"lines: {len(lines)} features (one per style/ref/direction)",
-    clipped = gpd.clip(lines, area_geom)
+          flush=True)
    print(f"clipped to area: {len(clipped)} line features remain", flush=True)
-    out_cols = [
+    pole_styles, pole_refs = build_pole_classes(
-        "mode", "style", "ref", "name", "network", "operator", "route",
+        routes, modes, trips, st, area_stop_ids
-        "colour", "bus_category", "colour_final", "geometry",
+    )
-    ]
+    stops, stations = build_stops_and_stations(
-    for c in out_cols:
+        pole_styles, pole_refs, area_geom, lines
-        if c not in clipped.columns:
+    )
-            clipped[c] = None
+    print(f"stops: {len(stops)} (bus + ferry); stations: {len(stations)}",
-    clipped = clipped.set_geometry("geometry")
+          flush=True)
    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")
+    lines.to_file(out, driver="GPKG", layer="lines")
-    print(f"wrote lines layer: {len(clipped)} features", flush=True)
+    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")
@@ -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 +
+Reads data/processed/master.gpkg (layers: lines, stops, stations — all
-config/styling.yaml, applies CLI filters, and composes a printable map:
+GTFS-derived by prepare.py) + area.gpkg + config/styling.yaml, applies CLI
-  - Carto Positron (no labels) basemap via contextily (EPSG:3857)
+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."""
+    """Style key for a station row, or None. Stations are pre-classified."""
-    if r.get("subway") == "yes" or r.get("station") == "subway":
+    s = r.get("style")
-        return "metro"
+    return s if s in RAIL_STYLES else None
    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
 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)
+        pts_by_style.setdefault(s, []).append((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)
    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
+        pts = pts_by_style[s]
        # 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)
        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
+    A bus stop whose name exactly matches a rail station is skipped — the
-    station marker represents that stop. Stops are grouped by shared
+    station marker represents it.
    route membership: same-name stops on disjoint routes are distinct.
    """
-    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 = rail_names or set()
        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
    mk = cfg.get("marker", {})
-    sizes = mk.get("size", {})
+    sz = mk.get("size", {}).get("bus", 1.5)
    sz = sizes.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 mk.get("shape", "circle") == "circle" else "s"
    marker = "o" if shape == "circle" else "s"
    top_z = max(styling["zorder"].values()) + 1
-    # group by name, then by route membership; skip rail station names
+    pts = [(r.geometry.x, r.geometry.y)
-    centroids = []
+           for _, r in bus.iterrows()
-    skipped = 0
+           if not (isinstance(r.get("name"), str) and r.get("name") in rail_names)]
-    for name, grp in bus.groupby("name"):
+    ax.scatter([p[0] for p in pts], [p[1] for p in pts],
        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],
               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):
+    # stations are unique per (name, style); dedupe by name for labelling
-        return classify_station(r) == "metro"
+    seen = set()
-
+    uniq = []
    def is_stog(r):
        return classify_station(r) == "s_tog"
    pts = []
    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")
@@ -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",
@@ -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"