From a6715eaff5228b6eadda3e1ee07aee3d1d3dcc17 Mon Sep 17 00:00:00 2001 From: marvin Date: Thu, 17 Sep 2026 21:07:40 +0200 Subject: [PATCH] Source Copenhagen transit data from Rejseplanen GTFS MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - 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 | 17 + copenhagen/config/modes.yaml | 60 ++-- copenhagen/config/styling.yaml | 27 +- copenhagen/scripts/_common.py | 1 - copenhagen/scripts/build_area.py | 21 +- copenhagen/scripts/download_gtfs.py | 35 +- copenhagen/scripts/download_osm.py | 302 ----------------- copenhagen/scripts/gtfs_to_geopackage.py | 149 +++++++-- copenhagen/scripts/prepare.py | 407 +++++++++++++++-------- copenhagen/scripts/render.py | 218 +++--------- pyproject.toml | 1 - uv.lock | 35 -- 12 files changed, 522 insertions(+), 751 deletions(-) delete mode 100644 copenhagen/scripts/download_osm.py diff --git a/copenhagen/config/area.json b/copenhagen/config/area.json index 72b1876..1f05ace 100644 --- a/copenhagen/config/area.json +++ b/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] + ] + } ] } diff --git a/copenhagen/config/modes.yaml b/copenhagen/config/modes.yaml index 04836ea..84d0f65 100644 --- a/copenhagen/config/modes.yaml +++ b/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] diff --git a/copenhagen/config/styling.yaml b/copenhagen/config/styling.yaml index da1276e..b43c4a6 100644 --- a/copenhagen/config/styling.yaml +++ b/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 diff --git a/copenhagen/scripts/_common.py b/copenhagen/scripts/_common.py index a54c0db..6541fcf 100644 --- a/copenhagen/scripts/_common.py +++ b/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" diff --git a/copenhagen/scripts/build_area.py b/copenhagen/scripts/build_area.py index d01df34..3a54379 100644 --- a/copenhagen/scripts/build_area.py +++ b/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" diff --git a/copenhagen/scripts/download_gtfs.py b/copenhagen/scripts/download_gtfs.py index 58680ed..e8b3c0b 100644 --- a/copenhagen/scripts/download_gtfs.py +++ b/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) diff --git a/copenhagen/scripts/download_osm.py b/copenhagen/scripts/download_osm.py deleted file mode 100644 index e797817..0000000 --- a/copenhagen/scripts/download_osm.py +++ /dev/null @@ -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() diff --git a/copenhagen/scripts/gtfs_to_geopackage.py b/copenhagen/scripts/gtfs_to_geopackage.py index 55b4436..acb6129 100644 --- a/copenhagen/scripts/gtfs_to_geopackage.py +++ b/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__": diff --git a/copenhagen/scripts/prepare.py b/copenhagen/scripts/prepare.py index 72ea5d2..be6bffb 100644 --- a/copenhagen/scripts/prepare.py +++ b/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") diff --git a/copenhagen/scripts/render.py b/copenhagen/scripts/render.py index 9ef3404..2722ea9 100644 --- a/copenhagen/scripts/render.py +++ b/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") diff --git a/pyproject.toml b/pyproject.toml index 6013be9..1cb5847 100644 --- a/pyproject.toml +++ b/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", diff --git a/uv.lock b/uv.lock index 71cfbb2..be5e719 100644 --- a/uv.lock +++ b/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" -- 2.54.0