marvin
dccf04349a
Replaces distance-based clustering with route-membership grouping: - download_osm.py: collect route_refs per stop node; use platform roles for bus stops, stop roles for rail stops - prepare.py: carry route_refs through to master.gpkg - render.py: group_by_route_membership() uses union-find to cluster same-name stops sharing at least one route into one marker; stops on disjoint routes are kept separate. No distance threshold needed. Results: 786 bus platforms -> 465 markers (was 1610 -> 781). Veksøvej: 2 platforms sharing routes 2A+5C -> 1 marker. Peter Bangs Vej: 7 platforms on routes 10/21/22/4A/7A -> 2 markers.
506 lines
18 KiB
Python
506 lines
18 KiB
Python
#!/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)
|
|
- 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
|
|
|
|
Usage:
|
|
uv run python render.py # default: all modes, PNG, markers only
|
|
uv run python render.py --labels # add station name labels
|
|
uv run python render.py --no-stops # no station markers
|
|
uv run python render.py --no-buses --format svg --out map.svg
|
|
uv run python render.py --bus-subset A,C,S # only trunk buses
|
|
uv run python render.py --modes metro,s_tog # rail-only
|
|
"""
|
|
import argparse
|
|
from pathlib import Path
|
|
|
|
import geopandas as gpd
|
|
import matplotlib
|
|
matplotlib.use("Agg")
|
|
import matplotlib.pyplot as plt
|
|
import numpy as np
|
|
import yaml
|
|
from affine import Affine
|
|
from adjustText import adjust_text
|
|
from scipy.ndimage import distance_transform_edt
|
|
|
|
from _common import CONFIG, OUTPUT, PROCESSED, UA
|
|
|
|
TARGET_CRS = "EPSG:3857"
|
|
|
|
|
|
def load_args():
|
|
p = argparse.ArgumentParser(description="Render the Copenhagen transit map.")
|
|
p.add_argument("--modes", default=None,
|
|
help="comma-separated styles to include (metro,s_tog,light_rail,regional,bus)")
|
|
p.add_argument("--no-buses", action="store_true", help="omit bus layer")
|
|
p.add_argument("--bus-subset", default=None,
|
|
help="comma-separated bus categories to keep (A,C,S,regular)")
|
|
p.add_argument("--max-bus-routes", type=int, default=None,
|
|
help="cap on number of distinct bus refs drawn")
|
|
p.add_argument("--dpi", type=int, default=200)
|
|
p.add_argument("--size", default="16",
|
|
help="figure size in inches: 'W' (height auto) or 'WxH'")
|
|
p.add_argument("--format", default="png", help="output format (png/svg/pdf)")
|
|
p.add_argument("--out", default=None, help="output path (default: output/copenhagen.<fmt>)")
|
|
p.add_argument("--no-basemap", action="store_true", help="skip basemap tiles (offline/faster)")
|
|
p.add_argument("--labels", action="store_true", help="draw station name labels (off by default)")
|
|
p.add_argument("--no-stops", action="store_true", help="skip station markers")
|
|
return p.parse_args()
|
|
|
|
|
|
def parse_size(spec, aspect):
|
|
if "x" in spec:
|
|
w, h = spec.split("x")
|
|
return float(w), float(h)
|
|
w = float(spec)
|
|
return w, w / aspect
|
|
|
|
|
|
def filter_lines(lines, args):
|
|
keep = set(lines["style"].unique())
|
|
if args.modes:
|
|
keep &= {m.strip() for m in args.modes.split(",")}
|
|
if args.no_buses:
|
|
keep.discard("bus")
|
|
out = lines[lines["style"].isin(keep)].copy()
|
|
|
|
if "bus" in keep and out["style"].eq("bus").any():
|
|
bus = out[out["style"].eq("bus")]
|
|
other = out[~out["style"].eq("bus")]
|
|
if args.bus_subset:
|
|
cats = {c.strip() for c in args.bus_subset.split(",")}
|
|
bus = bus[bus["bus_category"].isin(cats)]
|
|
if args.max_bus_routes is not None:
|
|
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)
|
|
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"]
|
|
try:
|
|
import contextily as cx
|
|
provider = styling["basemap"]["provider"]
|
|
src = cx.providers
|
|
for part in provider.split("."):
|
|
src = getattr(src, part)
|
|
cx.add_basemap(ax, crs=TARGET_CRS, source=src, zoom=zoom,
|
|
attribution=False, zorder=0,
|
|
headers={"User-Agent": UA})
|
|
return True
|
|
except Exception as e:
|
|
print(f" [basemap] unavailable ({e}); using plain background", flush=True)
|
|
ax.set_facecolor("#f2f2f0")
|
|
return False
|
|
|
|
|
|
def shapeburst_mask(ax, area_3857, extent, styling, res=1600):
|
|
fade_m = styling["mask"]["fade_distance_m"]
|
|
max_alpha = styling["mask"]["max_alpha"]
|
|
color = styling["mask"]["color"]
|
|
x0, x1, y0, y1 = extent
|
|
nx = res
|
|
ny = max(1, int(round(res * (y1 - y0) / (x1 - x0))))
|
|
pxw = (x1 - x0) / nx
|
|
pxh = (y1 - y0) / ny
|
|
transform = Affine.translation(x0, y1) * Affine.scale(pxw, -pxh)
|
|
|
|
from rasterio.features import rasterize
|
|
geom = area_3857.geometry.union_all()
|
|
mask = rasterize([(geom, 1)], out_shape=(ny, nx), transform=transform,
|
|
fill=0, dtype=np.uint8, all_touched=False).astype(bool)
|
|
dist_px = distance_transform_edt(~mask)
|
|
dist_m = dist_px * pxw
|
|
alpha = np.where(mask, 0.0, np.clip(dist_m / fade_m, 0.0, 1.0)) * max_alpha
|
|
|
|
cmap = matplotlib.colors.to_rgba(color)
|
|
rgba = np.zeros((ny, nx, 4))
|
|
rgba[..., :3] = cmap[:3]
|
|
rgba[..., 3] = alpha
|
|
ax.imshow(rgba, extent=(x0, x1, y0, y1), origin="upper",
|
|
aspect="equal", interpolation="bilinear", zorder=0.5)
|
|
|
|
|
|
def plot_lines(ax, lines, styling):
|
|
zord = styling["zorder"]
|
|
widths = styling["line_width"]
|
|
alphas = styling["alpha"]
|
|
outline_w = styling["outline_width"]
|
|
outline_c = styling["outline_color"]
|
|
for style in sorted(zord, key=lambda k: zord[k]):
|
|
grp = lines[lines["style"] == style]
|
|
if grp.empty:
|
|
continue
|
|
w = widths.get(style, 1.0)
|
|
z = zord[style]
|
|
# dark casing
|
|
grp.plot(ax=ax, color=outline_c, linewidth=w + outline_w,
|
|
zorder=z, alpha=1.0, capstyle="round", joinstyle="round")
|
|
# colour body
|
|
body_colors = grp["colour_final"].tolist()
|
|
grp.plot(ax=ax, color=body_colors, linewidth=w,
|
|
zorder=z + 0.1, alpha=alphas.get(style, 1.0),
|
|
capstyle="round", joinstyle="round")
|
|
|
|
|
|
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
|
|
|
|
|
|
def plot_stations(ax, stations, styling, active_styles):
|
|
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")
|
|
fill = mk.get("fill", "white")
|
|
edge = mk.get("edge", "#2b2b2b")
|
|
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)
|
|
|
|
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)
|
|
sz = sizes.get(s, 3.0)
|
|
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=1.0)
|
|
|
|
|
|
def rail_station_names(stations):
|
|
"""Return the set of names for rail stations (metro/s_tog/regional)."""
|
|
if stations is None or stations.empty:
|
|
return set()
|
|
names = set()
|
|
for _, r in stations.iterrows():
|
|
if classify_station(r) is not None:
|
|
name = r.get("name")
|
|
if isinstance(name, str) and name.strip():
|
|
names.add(name)
|
|
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]
|
|
if len(comp_coords) == 1:
|
|
centroids.append(comp_coords[0])
|
|
else:
|
|
# stops sharing a route are the same stop; merge to centroid
|
|
import numpy as np
|
|
centroids.append(tuple(np.mean(comp_coords, axis=0)))
|
|
return centroids
|
|
|
|
|
|
def plot_bus_stops(ax, stops, styling, active_styles, rail_names=None):
|
|
"""Draw small markers for bus stops from the stops layer.
|
|
|
|
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.
|
|
"""
|
|
if stops is None or stops.empty or "bus" not in 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"]
|
|
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
|
|
mk = cfg.get("marker", {})
|
|
sizes = mk.get("size", {})
|
|
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 shape == "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],
|
|
s=sz ** 2, marker=marker, c=fill, edgecolors=edge,
|
|
linewidths=lw, zorder=top_z, alpha=0.8)
|
|
|
|
|
|
def label_stations(ax, stations, styling, active_styles):
|
|
if stations is None or stations.empty:
|
|
return
|
|
cfg = styling.get("labels", {})
|
|
if not cfg.get("show", False):
|
|
return
|
|
label_styles = set(cfg.get("styles", ["metro", "s_tog"])) & 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 = []
|
|
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():
|
|
continue
|
|
pts.append((r.geometry.x, r.geometry.y, name, s == "metro"))
|
|
|
|
if not pts:
|
|
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
|
|
bbox = dict(boxstyle="round,pad=0.15", fc="white", ec="none", alpha=0.7)
|
|
t = ax.text(x, y, name, fontsize=fs, color="#111111",
|
|
zorder=10, ha="center", va="center", bbox=bbox)
|
|
texts.append(t)
|
|
try:
|
|
adjust_text(texts, ax=ax, expand_points=(1.4, 1.6),
|
|
force_text=(0.4, 0.6), lim=300,
|
|
arrowprops=dict(arrowstyle="-", color="#888", lw=0.4))
|
|
except Exception as e:
|
|
print(f" [labels] adjustText failed ({e}); labels may overlap", flush=True)
|
|
|
|
|
|
def main():
|
|
args = load_args()
|
|
styling = yaml.safe_load((CONFIG / "styling.yaml").read_text())
|
|
|
|
area = gpd.read_file(PROCESSED / "area.gpkg").to_crs(TARGET_CRS)
|
|
lines = gpd.read_file(PROCESSED / "master.gpkg", layer="lines").to_crs(TARGET_CRS)
|
|
try:
|
|
stations = gpd.read_file(PROCESSED / "master.gpkg", layer="stations").to_crs(TARGET_CRS)
|
|
except Exception:
|
|
stations = None
|
|
try:
|
|
stops = gpd.read_file(PROCESSED / "master.gpkg", layer="stops").to_crs(TARGET_CRS)
|
|
except Exception:
|
|
stops = None
|
|
|
|
lines = filter_lines(lines, args)
|
|
print(f"rendering {len(lines)} line features", flush=True)
|
|
|
|
bounds = area.total_bounds # (minx, miny, maxx, maxy)
|
|
dx = bounds[2] - bounds[0]
|
|
dy = bounds[3] - bounds[1]
|
|
aspect = dx / dy
|
|
pad = 0.12
|
|
extent = (bounds[0] - pad * dx, bounds[2] + pad * dx,
|
|
bounds[1] - pad * dy, bounds[3] + pad * dy)
|
|
|
|
fig_w, fig_h = parse_size(args.size, aspect)
|
|
fig, ax = plt.subplots(figsize=(fig_w, fig_h), dpi=args.dpi)
|
|
ax.set_xlim(extent[0], extent[1])
|
|
ax.set_ylim(extent[2], extent[3])
|
|
ax.set_aspect("equal")
|
|
ax.axis("off")
|
|
|
|
# basemap
|
|
if not args.no_basemap:
|
|
add_basemap(ax, styling)
|
|
|
|
# fade mask outside area
|
|
shapeburst_mask(ax, area, extent, styling)
|
|
|
|
# lines
|
|
plot_lines(ax, lines, styling)
|
|
|
|
# active styles (which mode layers are actually drawn)
|
|
active_styles = set(lines["style"].unique())
|
|
|
|
# station markers (symbols, no text)
|
|
if not args.no_stops:
|
|
plot_stations(ax, stations, styling, active_styles)
|
|
rail_names = rail_station_names(stations)
|
|
plot_bus_stops(ax, stops, styling, active_styles, rail_names=rail_names)
|
|
|
|
# station labels (opt-in via --labels)
|
|
if args.labels:
|
|
styling["labels"]["show"] = True
|
|
if styling.get("labels", {}).get("show", False):
|
|
label_stations(ax, stations, styling, active_styles)
|
|
|
|
# title + attribution
|
|
ax.text(0.5, 0.985, "Kurragömma Köpenhamn",
|
|
transform=ax.transAxes, ha="center", va="top",
|
|
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",
|
|
transform=ax.transAxes, ha="left", va="bottom",
|
|
fontsize=max(5, fig_w * 0.5), color="#666")
|
|
|
|
out_dir = OUTPUT
|
|
out_dir.mkdir(parents=True, exist_ok=True)
|
|
if args.out:
|
|
out_path = Path(args.out)
|
|
if not out_path.is_absolute():
|
|
out_path = out_dir / out_path
|
|
else:
|
|
out_path = out_dir / f"copenhagen.{args.format}"
|
|
out_path.parent.mkdir(parents=True, exist_ok=True)
|
|
|
|
plt.savefig(out_path, dpi=args.dpi, bbox_inches="tight",
|
|
pad_inches=0.2, facecolor="white")
|
|
plt.close(fig)
|
|
print(f"wrote {out_path} ({args.format}, {fig_w}x{fig_h}in @ {args.dpi}dpi)")
|
|
|
|
|
|
if __name__ == "__main__":
|
|
main()
|