pbt/scripts/fetch_stops.py

#!/usr/bin/env python3
"""Fetch all Austrian public-transport stops (and the lines serving them) from
OpenStreetMap via the Overpass API and write a deduplicated snapshot to
data/stops_at.csv.gz.

Data © OpenStreetMap contributors, licensed under the ODbL.

Run this only when you want to refresh the snapshot; the import into the
database reads the checked-in CSV and needs no network access. The route pass
downloads a few hundred MB in ~25 tiles and needs roughly 1 GB of RAM.

    python scripts/fetch_stops.py
"""
from __future__ import annotations

import csv
import gzip
import json
import math
import re
import sys
import time
import unicodedata
import urllib.parse
import urllib.request
from pathlib import Path

OUT_PATH = Path(__file__).resolve().parent.parent / "data" / "stops_at.csv.gz"

OVERPASS_ENDPOINTS = [
    "https://overpass-api.de/api/interpreter",
    "https://overpass.kumi.systems/api/interpreter",
    "https://overpass.private.coffee/api/interpreter",
]

# Named nodes in Austria that represent a boardable stop ("Haltestelle").
STOPS_QUERY = """
[out:json][timeout:600];
area["ISO3166-1"="AT"][admin_level=2]->.at;
(
  node(area.at)["highway"="bus_stop"]["name"];
  node(area.at)["railway"="tram_stop"]["name"];
  node(area.at)["railway"="station"]["name"];
  node(area.at)["railway"="halt"]["name"];
  node(area.at)["public_transport"="station"]["name"];
  node(area.at)["amenity"="bus_station"]["name"];
);
out body;
"""

ROUTE_MODES = "bus|trolleybus|tram|light_rail|subway|train|monorail|share_taxi"

# Austria bbox, tiled so each Overpass response stays small enough to parse.
BBOX = (46.35, 9.50, 49.05, 17.20)  # south, west, north, east
LAT_STEP = 0.9
LON_STEP = 1.0

STOP_ROLES = {
    "stop", "platform",
    "stop_entry_only", "stop_exit_only",
    "platform_entry_only", "platform_exit_only",
}

CSV_FIELDS = [
    "osm_type", "osm_id", "name", "stop_type",
    "municipality", "latitude", "longitude", "lines",
]

TYPE_RANK = {"subway": 5, "train": 4, "bus_station": 3, "tram": 2, "bus": 1, "other": 0}

# OSM keeps a separate node per direction for most bus stops (~20-60 m apart);
# nodes with the same name closer than this are treated as one stop.
MERGE_METERS = 250.0


def normalize_name(value: str) -> str:
    """Lowercase, fold accents and ß so search is diacritic-insensitive.

    The same function lives in stops_import.py; keep them in sync.
    """
    value = value.strip().lower().replace("ß", "ss")
    decomposed = unicodedata.normalize("NFKD", value)
    return "".join(ch for ch in decomposed if not unicodedata.combining(ch))


def haversine_m(lat1: float, lon1: float, lat2: float, lon2: float) -> float:
    rlat1, rlon1, rlat2, rlon2 = map(math.radians, (lat1, lon1, lat2, lon2))
    dlat, dlon = rlat2 - rlat1, rlon2 - rlon1
    h = math.sin(dlat / 2) ** 2 + math.cos(rlat1) * math.cos(rlat2) * math.sin(dlon / 2) ** 2
    return 2 * 6_371_000 * math.asin(min(1.0, math.sqrt(h)))


def classify(tags: dict) -> str:
    if tags.get("station") == "subway" or tags.get("subway") == "yes":
        return "subway"
    if tags.get("railway") in {"station", "halt"}:
        return "train"
    if tags.get("railway") == "tram_stop" or tags.get("tram") == "yes":
        return "tram"
    if tags.get("amenity") == "bus_station":
        return "bus_station"
    if tags.get("highway") == "bus_stop" or tags.get("bus") == "yes":
        return "bus"
    return "other"


def municipality_of(tags: dict) -> str:
    for key in ("addr:city", "is_in:municipality", "is_in:city", "is_in"):
        value = tags.get(key)
        if value:
            return value.split(",")[0].strip()
    return ""


def overpass(query: str, *, tries: int = 4) -> dict:
    body = urllib.parse.urlencode({"data": query}).encode()
    last_error: Exception | None = None
    for attempt in range(tries):
        endpoint = OVERPASS_ENDPOINTS[attempt % len(OVERPASS_ENDPOINTS)]
        try:
            req = urllib.request.Request(
                endpoint, data=body, headers={"User-Agent": "pbt-stop-import/1.0"}
            )
            with urllib.request.urlopen(req, timeout=600) as resp:
                return json.load(resp)
        except Exception as exc:  # noqa: BLE001 - retry on another mirror
            last_error = exc
            print(f"  {endpoint} failed ({exc}); retrying", file=sys.stderr)
            time.sleep(10 * (attempt + 1))
    raise SystemExit(f"Overpass failed after {tries} tries; last error: {last_error}")


# --------------------------------------------------------------------------- #
# Stops
# --------------------------------------------------------------------------- #
def build_stop_rows(payload: dict) -> list[dict]:
    by_name: dict[str, list[dict]] = {}
    for element in payload.get("elements", []):
        if element.get("type") != "node":
            continue
        tags = element.get("tags", {})
        name = (tags.get("name") or "").strip()
        lat, lon = element.get("lat"), element.get("lon")
        if not name or lat is None or lon is None:
            continue
        node = {
            "osm_id": int(element["id"]),
            "name": name,
            "stop_type": classify(tags),
            "municipality": municipality_of(tags),
            "latitude": round(float(lat), 6),
            "longitude": round(float(lon), 6),
        }
        by_name.setdefault(normalize_name(name), []).append(node)

    rows: list[dict] = []
    for group in by_name.values():
        clusters: list[list[dict]] = []
        for node in sorted(group, key=lambda n: n["osm_id"]):
            for cluster in clusters:
                if any(
                    haversine_m(
                        node["latitude"], node["longitude"],
                        other["latitude"], other["longitude"],
                    ) <= MERGE_METERS
                    for other in cluster
                ):
                    cluster.append(node)
                    break
            else:
                clusters.append([node])

        for cluster in clusters:
            lead = max(
                cluster,
                key=lambda n: (TYPE_RANK[n["stop_type"]], n["municipality"] != "", -n["osm_id"]),
            )
            rows.append({
                "osm_type": "node",
                "osm_id": min(n["osm_id"] for n in cluster),
                "name": lead["name"],
                "stop_type": lead["stop_type"],
                "municipality": lead["municipality"]
                or next((n["municipality"] for n in cluster if n["municipality"]), ""),
                "latitude": round(sum(n["latitude"] for n in cluster) / len(cluster), 6),
                "longitude": round(sum(n["longitude"] for n in cluster) / len(cluster), 6),
                "lines": set(),
            })

    rows.sort(key=lambda r: (normalize_name(r["name"]), r["osm_id"]))
    return rows


# --------------------------------------------------------------------------- #
# Lines (route relations -> stops)
# --------------------------------------------------------------------------- #
def fetch_route_members() -> tuple[list[tuple[str, int]], dict[int, tuple[str, float, float]]]:
    """Return (ref, node_id) pairs plus node_id -> (name, lat, lon)."""
    seen_rel: set[int] = set()
    ref_members: list[tuple[str, int]] = []
    node_pos: dict[int, tuple[str, float, float]] = {}

    south, west, north, east = BBOX
    lat = south
    tiles: list[tuple[float, float, float, float]] = []
    while lat < north:
        lon = west
        while lon < east:
            tiles.append((lat, lon, min(lat + LAT_STEP, north), min(lon + LON_STEP, east)))
            lon += LON_STEP
        lat += LAT_STEP

    for i, (s, w, n, e) in enumerate(tiles, 1):
        print(f"  route tile {i}/{len(tiles)} ({s:.1f},{w:.1f})", file=sys.stderr)
        query = (
            f"[out:json][timeout:400];"
            f'relation["type"="route"]["route"~"^({ROUTE_MODES})$"]["ref"]'
            f"({s},{w},{n},{e})->.r;"
            f".r out body;"
            f"node(r.r);out body;"
        )
        payload = overpass(query)
        for el in payload.get("elements", []):
            if el["type"] == "node":
                tags = el.get("tags", {})
                node_pos[el["id"]] = (tags.get("name", ""), el["lat"], el["lon"])
            elif el["type"] == "relation":
                if el["id"] in seen_rel:
                    continue
                seen_rel.add(el["id"])
                # A few relations pack several numbers into one ref ("407, 413").
                refs = [r.strip() for r in re.split(r"[;,]", el["tags"]["ref"]) if r.strip()]
                if not refs:
                    continue
                members = [
                    m["ref"] for m in el.get("members", [])
                    if m["type"] == "node" and m["role"] in STOP_ROLES
                ]
                for ref in refs:
                    for node_id in members:
                        ref_members.append((ref, node_id))
        del payload

    print(f"  {len(seen_rel)} routes, {len(ref_members)} stop memberships", file=sys.stderr)
    return ref_members, node_pos


def assign_lines(rows: list[dict], ref_members, node_pos) -> tuple[int, int]:
    by_name: dict[str, list[int]] = {}
    grid: dict[tuple[int, int], list[int]] = {}
    for idx, row in enumerate(rows):
        by_name.setdefault(normalize_name(row["name"]), []).append(idx)
        cell = (round(row["latitude"] / 0.01), round(row["longitude"] / 0.01))
        grid.setdefault(cell, []).append(idx)

    def nearest(lat: float, lon: float, max_m: float) -> int | None:
        best, best_d = None, max_m
        base = (round(lat / 0.01), round(lon / 0.01))
        for dr in (-1, 0, 1):
            for dc in (-1, 0, 1):
                for idx in grid.get((base[0] + dr, base[1] + dc), ()):
                    d = haversine_m(lat, lon, rows[idx]["latitude"], rows[idx]["longitude"])
                    if d < best_d:
                        best, best_d = idx, d
        return best

    matched = unmatched = 0
    for ref, node_id in ref_members:
        pos = node_pos.get(node_id)
        if pos is None:
            unmatched += 1
            continue
        name, lat, lon = pos
        idx = None
        key = normalize_name(name) if name else ""
        if key and key in by_name:
            idx = min(
                by_name[key],
                key=lambda i: haversine_m(lat, lon, rows[i]["latitude"], rows[i]["longitude"]),
            )
            if haversine_m(lat, lon, rows[idx]["latitude"], rows[idx]["longitude"]) > 600:
                idx = None
        if idx is None:
            idx = nearest(lat, lon, 120)
        if idx is None:
            unmatched += 1
        else:
            matched += 1
            rows[idx]["lines"].add(ref)
    return matched, unmatched


def line_sort_key(ref: str):
    head, i = "", 0
    while i < len(ref) and not ref[i].isdigit():
        head, i = head + ref[i], i + 1
    digits = ""
    while i < len(ref) and ref[i].isdigit():
        digits, i = digits + ref[i], i + 1
    return (head.lower(), int(digits) if digits else -1, ref[i:].lower())


# --------------------------------------------------------------------------- #
def main() -> None:
    print("Fetching stops ...", file=sys.stderr)
    rows = build_stop_rows(overpass(STOPS_QUERY))
    if len(rows) < 10_000:
        raise SystemExit(
            f"Only {len(rows)} stops parsed - refusing to overwrite the snapshot."
        )

    print("Fetching routes ...", file=sys.stderr)
    ref_members, node_pos = fetch_route_members()
    matched, unmatched = assign_lines(rows, ref_members, node_pos)

    OUT_PATH.parent.mkdir(parents=True, exist_ok=True)
    with gzip.open(OUT_PATH, "wt", newline="", encoding="utf-8") as fh:
        writer = csv.DictWriter(fh, fieldnames=CSV_FIELDS)
        writer.writeheader()
        for row in rows:
            row["lines"] = ";".join(sorted(row["lines"], key=line_sort_key))
            writer.writerow(row)

    by_type: dict[str, int] = {}
    for row in rows:
        by_type[row["stop_type"]] = by_type.get(row["stop_type"], 0) + 1
    with_lines = sum(1 for r in rows if r["lines"])
    print(f"Wrote {len(rows)} stops to {OUT_PATH}", file=sys.stderr)
    print(f"  by type: {by_type}", file=sys.stderr)
    print(
        f"  lines: {with_lines} stops have >=1 line "
        f"(memberships matched {matched}, unmatched {unmatched})",
        file=sys.stderr,
    )


if __name__ == "__main__":
    main()