diff --git a/copenhagen/scripts/download_osm.py b/copenhagen/scripts/download_osm.py
index f55b517..559c74c 100644
--- a/copenhagen/scripts/download_osm.py
+++ b/copenhagen/scripts/download_osm.py
@@ -126,14 +126,64 @@ def parse(data, mode):
     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 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)
     for mode, cfg in modes.items():
-        ql = build_query(cfg, (s, n, w, e))
+        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)
@@ -146,6 +196,18 @@ def main():
         if stops:
             sgdf = gpd.GeoDataFrame(stops, crs="EPSG:4326")
             sgdf.to_file(OSM_RAW / f"{mode}.gpkg", driver="GPKG", layer="stops")
+
+    # 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.")
 
 
diff --git a/copenhagen/scripts/prepare.py b/copenhagen/scripts/prepare.py
new file mode 100644
index 0000000..e325fae
--- /dev/null
+++ b/copenhagen/scripts/prepare.py
@@ -0,0 +1,205 @@
+#!/usr/bin/env python3
+"""Merge OSM 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".
+
+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.
+"""
+import re
+from pathlib import Path
+
+import geopandas as gpd
+import pandas as pd
+import yaml
+
+from _common import CONFIG, OSM_RAW, PROCESSED
+
+MIXED_CRS_WARN = "mixed CRS"
+
+
+def load_modes():
+    return yaml.safe_load((CONFIG / "modes.yaml").read_text())["modes"]
+
+
+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 collect_stops(modes):
+    frames = []
+    for mode, cfg in modes.items():
+        p = OSM_RAW / f"{mode}.gpkg"
+        if not p.exists():
+            continue
+        try:
+            g = gpd.read_file(p, layer="stops")
+        except Exception:
+            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)
+
+
+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 categorise_bus(ref, patterns):
+    if ref is None:
+        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):
+    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
+    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")
+
+
+def main():
+    modes = load_modes()
+    styling = load_styling()
+    palette = styling["palette"]
+    bus_patterns = styling["bus_filters"]["categories"]
+
+    area = gpd.read_file(PROCESSED / "area.gpkg")
+    area_geom = area.geometry.union_all()
+
+    lines = collect_lines(modes)
+    print(f"collected {len(lines)} raw lines across modes", flush=True)
+
+    # 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"])
+    ]
+
+    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)
+
+    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)
+
+    # 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", "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"):
+        refs = grp["ref"].dropna().unique()
+        print(f"  {style:10s}: {len(grp):4d} feats, {len(refs):3d} refs")
+
+
+if __name__ == "__main__":
+    main()