Add road closure cost scripts for opus and sonnet implementations

data/

.venv/

venv/

__pycache__/

*.pyc

cache/

# Generated outputs: keep only the HTML pages in git

output/*

!output/*.html

"""Step 1: traffic counts at the BASt station next to the bridge.

- Download the BASt station list (Jawe2023, with coordinates) once.

- Stream the 1.4 GB hourly file with DuckDB and keep the rows of one station.

- Compute the annual average daily traffic (DTV) per direction, Pkw / Lkw,

  and the average hourly profile.

Output: data/<bridge>/counts_hourly.parquet, data/<bridge>/dtv.json,

        data/stations_2023.csv (all stations with coordinates, used by the BPR step)

"""

import json

import sys

import urllib.request

from pathlib import Path

sys.path.insert(0, str(Path(__file__).resolve().parent))

import duckdb

import pandas as pd

from common import DATA, haversine_km, load_config, run_dir

CLASSES = ["KFZ", "Lkw", "PLZ", "Pkw", "Lfw", "Mot", "PmA", "Bus", "LoA", "Lzg", "Sat", "Son"]

def download(url, path):

    if not path.exists():

        print(f"Downloading {url}")

        req = urllib.request.Request(url, headers={"User-Agent": "Mozilla/5.0 (research script)"})

        with urllib.request.urlopen(req, timeout=120) as r:

            path.write_bytes(r.read())

    return path

def read_station_list(cfg):

    """BASt station list. Semicolon separated, latin-1, German decimal comma."""

    raw = download(cfg["counts"]["station_list_url"], DATA / "stations" / "Jawe2023.csv")

    st = pd.read_csv(raw, sep=";", encoding="latin-1", dtype=str)

    keep = {

        "TK_Nr": "tknr", "DZ_Nr": "zst", "DZ_Name": "name", "Land_Code": "land",

        "Str_Kl": "strklas", "Str_Nr": "strnum", "Fernziel_Ri1": "fernziel_r1",

        "Nahziel_Ri1": "nahziel_r1", "Fernziel_Ri2": "fernziel_r2", "Nahziel_Ri2": "nahziel_r2",

        "Anz_Fs_Q": "lanes", "DTV_Kfz_MobisSo_Q": "dtv_kfz", "DTV_SV_MobisSo_Q": "dtv_sv",

        "Koor_WGS84_N": "lat", "Koor_WGS84_E": "lon",

    }

    st = st[list(keep)].rename(columns=keep)

    for c in ["dtv_kfz", "dtv_sv"]:

        st[c] = pd.to_numeric(st[c].str.replace(".", "", regex=False), errors="coerce")

    for c in ["lat", "lon"]:

        st[c] = pd.to_numeric(st[c].str.replace(",", ".", regex=False), errors="coerce")

    st["tknr"] = pd.to_numeric(st["tknr"], errors="coerce").astype("Int64")

    st.to_csv(DATA / "stations_2023.csv", index=False)

    return st

def read_hourly(cfg, tknr, out_path):

    """Filter one station out of the big hourly file with DuckDB (streaming)."""

    if out_path.exists():

        return pd.read_parquet(out_path)

    src = cfg["counts"]["hourly_file"]

    print(f"Scanning {src} for station {tknr} (takes a minute)")

    con = duckdb.connect()

    df = con.execute(

        """

        SELECT * FROM read_csv(?, delim=';', header=true, all_varchar=true)

        WHERE trim(TKNR) = ?

        """,

        [src, str(tknr)],

    ).df()

    df.columns = [c.strip() for c in df.columns]

    df = df.apply(lambda s: s.str.strip())

    for c in df.columns:

        if not c.startswith("K_") and c not in ("Land", "Strklas", "Fahrtzw"):

            df[c] = pd.to_numeric(df[c], errors="coerce")

    # -1 with flag x means missing

    for cls in CLASSES:

        for r in ("R1", "R2"):

            col = f"{cls}_{r}"

            bad = (df[col] < 0) | (df[f"K_{col}"] == "x")

            df.loc[bad, col] = pd.NA

    df.to_parquet(out_path)

    return df

def daily_totals(df):

    """Daily sums per direction; keep only days with 24 complete hours."""

    df = df.copy()

    for r in ("R1", "R2"):

        df[f"PkwGrp_{r}"] = df[f"KFZ_{r}"] - df[f"Lkw_{r}"]  # all vehicles <= 3.5 t

    cols = [f"{c}_{r}" for c in ["KFZ", "Lkw", "PkwGrp", "Bus", "LoA", "Lzg", "Sat"] for r in ("R1", "R2")]

    g = df.groupby(["Datum", "Fahrtzw"])

    day = g[cols].sum(min_count=24)

    n_ok = g[["KFZ_R1", "KFZ_R2", "Lkw_R1", "Lkw_R2"]].count().min(axis=1)

    day = day[n_ok == 24].dropna(subset=["KFZ_R1", "KFZ_R2", "Lkw_R1", "Lkw_R2"])

    return day.reset_index()

def main():

    cfg = load_config()

    out = run_dir(cfg)

    tknr = cfg["counts"]["station_tknr"]

    stations = read_station_list(cfg)

    st = stations[stations["tknr"] == tknr].iloc[0]

    df = read_hourly(cfg, tknr, DATA / f"counts_hourly_{tknr}.parquet")

    print(f"{len(df)} hourly rows, {df['Datum'].nunique()} days")

    # The "Lkw" column matches the published DTV_SV (heavy vehicles). Its sum does not equal

    # LoA + Lzg + Sat exactly, so the subclasses are only used for the articulated share.

    day = daily_totals(df)

    if cfg["counts"]["days"] == "workdays":

        day = day[day["Fahrtzw"] == "w"]

    print(f"{len(day)} complete days used (days = {cfg['counts']['days']})")

    m = day.mean(numeric_only=True)

    # Hourly profile: mean share of daily traffic by hour (both directions together)

    hours = df.groupby("Stunde")[["KFZ_R1", "KFZ_R2", "Lkw_R1", "Lkw_R2"]].mean()

    profile_kfz = (hours["KFZ_R1"] + hours["KFZ_R2"]) / (hours["KFZ_R1"] + hours["KFZ_R2"]).sum()

    profile_lkw = (hours["Lkw_R1"] + hours["Lkw_R2"]) / (hours["Lkw_R1"] + hours["Lkw_R2"]).sum()

    art = m["Lzg_R1"] + m["Lzg_R2"] + m["Sat_R1"] + m["Sat_R2"]

    trucks = art + m["LoA_R1"] + m["LoA_R2"]

    bridge = cfg["bridge_info"]

    res = {

        "station": {

            "tknr": int(tknr), "zst": st["zst"], "name": st["name"], "lat": float(st["lat"]),

            "lon": float(st["lon"]), "road": f"{st['strklas']}{st['strnum']}",

            "dir1": f"{st['nahziel_r1']} / {st['fernziel_r1']}",

            "dir2": f"{st['nahziel_r2']} / {st['fernziel_r2']}",

            "dtv_kfz_bast": float(st["dtv_kfz"]), "dtv_sv_bast": float(st["dtv_sv"]),

            "distance_to_bridge_km": haversine_km(st["lat"], st["lon"], bridge["lat"], bridge["lon"]),

        },

        "days_used": int(len(day)),

        "days_rule": cfg["counts"]["days"],

        "dtv": {

            "pkw": {"R1": m["PkwGrp_R1"], "R2": m["PkwGrp_R2"]},

            "lkw": {"R1": m["Lkw_R1"], "R2": m["Lkw_R2"]},

        },

        "dtv_kfz_total": m["KFZ_R1"] + m["KFZ_R2"],

        "lkw_articulated_share": art / trucks,

        "hourly_profile_kfz": profile_kfz.round(5).tolist(),

        "hourly_profile_lkw": profile_lkw.round(5).tolist(),

    }

    override = bridge.get("dtv_override")

    if override:

        # Keep the station's direction split, rescale to the published volume at the bridge

        targets = {"pkw": override["kfz"] - override["lkw"], "lkw": override["lkw"]}

        for cls, target in targets.items():

            tot = res["dtv"][cls]["R1"] + res["dtv"][cls]["R2"]

            res["dtv"][cls] = {r: v * target / tot for r, v in res["dtv"][cls].items()}

        res["dtv_kfz_total"] = float(override["kfz"])

        res["dtv_source"] = override["source"]

        print(f"DTV overridden for the bridge: {override['kfz']} Kfz, {override['lkw']} Lkw ({override['source']})")

    else:

        res["dtv_source"] = f"BASt station {tknr}, 2023 hourly counts"

    res["dtv_pkw_total"] = res["dtv"]["pkw"]["R1"] + res["dtv"]["pkw"]["R2"]

    res["dtv_lkw_total"] = res["dtv"]["lkw"]["R1"] + res["dtv"]["lkw"]["R2"]

    (out / "dtv.json").write_text(json.dumps(res, indent=2, default=float))

    print(json.dumps({k: res[k] for k in ["dtv", "dtv_kfz_total", "lkw_articulated_share"]}, indent=2, default=float))

    print(f"BASt published DTV Kfz {st['dtv_kfz']:.0f}, SV {st['dtv_sv']:.0f}")

if __name__ == "__main__":

    main()

"""Step 2: road network and origin-destination points.

- Build the OSM drive network around the bridge with OSMnx (from Geofabrik

  extracts cut with pyosmium, or from Overpass): main roads

  (motorway ... unclassified) up to the largest buffer plus a margin, and all

  drivable streets close to the bridge. Add free-flow speeds and travel times

  for Pkw and Lkw.

- Find the bridge edges in the graph (by OSM way id).

- Read municipality population and centroids (Destatis) and snap the

  municipalities inside the largest buffer to the network.

Output: data/<bridge>/graph.graphml, data/<bridge>/bridge_edges.json,

        data/<bridge>/od_points.csv

"""

import html

import json

import re

import sys

import urllib.request

from pathlib import Path

sys.path.insert(0, str(Path(__file__).resolve().parent))

import networkx as nx

import osmnx as ox

import pandas as pd

from common import DATA, haversine_km, load_config, run_dir

MAIN_ROADS = ('["highway"~"motorway|motorway_link|trunk|trunk_link|primary|primary_link|'

              'secondary|secondary_link|tertiary|tertiary_link|unclassified"]')

ox.settings.use_cache = True

ox.settings.cache_folder = str(DATA / "osmnx_cache")

ox.settings.log_console = False

ox.settings.requests_timeout = 300

# Many small Overpass queries (20 x 20 km) instead of a few huge ones: public servers time out otherwise

ox.settings.max_query_area_size = 20_000 * 20_000

def first(x):

    return x[0] if isinstance(x, list) else x

def parse_maxspeed(v):

    """OSM maxspeed to km/h, or None (e.g. 'none', 'signals', 'DE:rural')."""

    v = first(v)

    if v is None or (isinstance(v, float) and pd.isna(v)):

        return None

    v = str(v)

    if v in ("DE:urban",):

        return 50.0

    if v in ("DE:rural",):

        return 100.0

    m = re.match(r"^\s*(\d+)", v)

    return float(m.group(1)) if m else None

def add_speeds(G, cfg):

    defaults = cfg["network"]["default_speeds_kmh"]

    cap = cfg["network"]["lkw_max_speed_kmh"]

    for u, v, k, d in G.edges(keys=True, data=True):

        hw = first(d.get("highway", "unclassified"))

        speed = parse_maxspeed(d.get("maxspeed"))

        if speed is None:

            speed = defaults.get(hw, 50)

        # 'none' on motorways (no limit): use the default motorway speed

        speed = min(speed, defaults.get("motorway", 120)) if hw.startswith("motorway") else speed

        lkw_cap = cap["motorway"] if hw.startswith("motorway") else cap["other"]

        d["speed_kph"] = speed

        d["travel_time"] = d["length"] / (speed / 3.6)

        d["travel_time_lkw"] = d["length"] / (min(speed, lkw_cap) / 3.6)

    return G

MAIN_TYPES = {"motorway", "motorway_link", "trunk", "trunk_link", "primary", "primary_link", "secondary",

              "secondary_link", "tertiary", "tertiary_link", "unclassified"}

LOCAL_TYPES = MAIN_TYPES | {"residential", "living_street", "road"}

KEEP_TAGS = ("highway", "maxspeed", "oneway", "name", "ref", "bridge", "junction", "lanes")

def extract_osm_xml(cfg, xml_path):

    """Cut the drivable roads around the bridge out of Geofabrik .pbf extracts into a small .osm XML.

    Main road types are kept within dist_main of the bridge, residential streets only within

    local_radius_km. A way is kept when one of its nodes lies inside the radius (all its nodes are kept).

    """

    import osmium

    b = cfg["bridge_info"]

    net = cfg["network"]

    r_main = max(net["sensitivity_buffers_km"] + [net["buffer_km"]]) + net["graph_margin_km"]

    r_loc = net["local_radius_km"]

    pbfs = []

    for region in net["geofabrik_regions"]:

        f = DATA / "osm" / (region.split("/")[-1] + ".osm.pbf")

        if not f.exists():

            f.parent.mkdir(exist_ok=True)

            url = f"https://download.geofabrik.de/{region}-latest.osm.pbf"

            print(f"Downloading {url}")

            urllib.request.urlretrieve(url, f)

        pbfs.append(f)

    def drivable(tags):

        if tags.get("area") == "yes" or tags.get("motor_vehicle") in ("no", "private"):

            return False

        if tags.get("access") in ("no", "private") and tags.get("motor_vehicle") not in ("yes", "designated"):

            return False

        return tags.get("highway") in LOCAL_TYPES

    # One pass per file: osmium stores node locations in C++, only highway ways reach Python

    coords, out_ways, seen = {}, [], set()  # ways on a border are in two extracts

    for f in pbfs:

        print(f"Reading roads from {f.name} (a few minutes)")

        fp = (osmium.FileProcessor(str(f)).with_locations("sparse_mem_array")

              .with_filter(osmium.filter.EntityFilter(osmium.osm.WAY))

              .with_filter(osmium.filter.KeyFilter("highway")))

        for w in fp:

            if w.id in seen or not drivable(w.tags):

                continue

            hw = w.tags["highway"]

            radius = r_main if hw in MAIN_TYPES else r_loc

            try:

                pts = [(n.ref, n.location.lat, n.location.lon) for n in w.nodes]

            except osmium.InvalidLocationError:

                continue

            if not any(haversine_km(la, lo, b["lat"], b["lon"]) <= radius for _, la, lo in pts[::5] + pts[-1:]):

                continue

            seen.add(w.id)

            out_ways.append((w.id, [n for n, _, _ in pts], {k: w.tags[k] for k in KEEP_TAGS if k in w.tags}))

            for n, la, lo in pts:

                coords[n] = (la, lo)

    kept_nodes = coords.keys()

    print(f"Writing {len(out_ways)} ways, {len(kept_nodes)} nodes to {xml_path.name}")

    esc = html.escape

    with open(xml_path, "w", encoding="utf-8") as fh:

        fh.write('<?xml version="1.0" encoding="UTF-8"?>\n<osm version="0.6">\n')

        for nid in kept_nodes:

            la, lo = coords[nid]

            fh.write(f'<node id="{nid}" lat="{la:.7f}" lon="{lo:.7f}"/>\n')

        for wid, refs, tags in out_ways:

            fh.write(f'<way id="{wid}">')

            fh.write("".join(f'<nd ref="{n}"/>' for n in refs))

            fh.write("".join(f'<tag k="{k}" v="{esc(v, quote=True)}"/>' for k, v in tags.items()))

            fh.write("</way>\n")

        fh.write("</osm>\n")

def build_graph(cfg, path):

    if path.exists():

        return ox.load_graphml(path, edge_dtypes={"travel_time": float, "travel_time_lkw": float,

                                                   "speed_kph": float, "length": float})

    b = cfg["bridge_info"]

    net = cfg["network"]

    if net["source"] == "geofabrik":

        xml = path.with_name("roads.osm")

        if not xml.exists():

            extract_osm_xml(cfg, xml)

        G = ox.graph_from_xml(xml, simplify=False, retain_all=True)

    else:

        if net.get("overpass_url"):

            ox.settings.overpass_url = net["overpass_url"]

        dist_main = (max(net["sensitivity_buffers_km"] + [net["buffer_km"]]) + net["graph_margin_km"]) * 1000

        print(f"Downloading main roads within {dist_main / 1000:.0f} km from Overpass")

        G_main = ox.graph_from_point((b["lat"], b["lon"]), dist=dist_main, custom_filter=MAIN_ROADS,

                                     simplify=False, retain_all=True)

        G_loc = ox.graph_from_point((b["lat"], b["lon"]), dist=net["local_radius_km"] * 1000,

                                    network_type="drive", simplify=False)

        G = nx.compose(G_main, G_loc)

    # Speeds on the raw OSM segments first, then merge segments and sum their travel times

    # (otherwise a merged edge with maxspeed [60, 100] would get one speed for its whole length)

    G = add_speeds(G, cfg)

    G = ox.simplify_graph(G, edge_attr_aggs={"length": sum, "travel_time": sum, "travel_time_lkw": sum})

    G = ox.truncate.largest_component(G, strongly=True)

    for u, v, k, d in G.edges(keys=True, data=True):

        d["speed_kph"] = 3.6 * d["length"] / d["travel_time"]

    ox.save_graphml(G, path)

    return G

def find_bridge_edges(G, cfg):

    b = cfg["bridge_info"]

    ids = set(b.get("osm_way_ids") or [])

    edges = []

    for u, v, k, d in G.edges(keys=True, data=True):

        osmid = d.get("osmid")

        osmids = osmid if isinstance(osmid, list) else [osmid]

        if ids & {int(x) for x in osmids}:

            edges.append((u, v, k))

    if not edges:

        # Fallback: the nearest edge to the bridge point and its reverse

        u, v, k = ox.distance.nearest_edges(G, b["lon"], b["lat"])

        edges = [(u, v, k)] + [(v, u, kk) for kk in G.get_edge_data(v, u, default={})]

    return edges

def read_municipalities(cfg):

    path = DATA / "destatis" / "gv_31122023.xlsx"

    if not path.exists():

        path.parent.mkdir(exist_ok=True)

        req = urllib.request.Request(cfg["population"]["url"], headers={"User-Agent": "Mozilla/5.0"})

        with urllib.request.urlopen(req, timeout=120) as r:

            path.write_bytes(r.read())

    raw = pd.read_excel(path, sheet_name=1, header=None, dtype=str)

    g = raw[raw[0] == "60"]  # Satzart 60 = municipality

    df = pd.DataFrame({

        "ags": g[2] + g[3] + g[4] + g[6],

        "name": g[7],

        "pop": pd.to_numeric(g[9], errors="coerce"),

        "lon": pd.to_numeric(g[14].str.replace(",", "."), errors="coerce"),

        "lat": pd.to_numeric(g[15].str.replace(",", "."), errors="coerce"),

    }).dropna()

    return df

def main():

    cfg = load_config()

    out = run_dir(cfg)

    b = cfg["bridge_info"]

    G = build_graph(cfg, out / "graph.graphml")

    print(f"Graph: {G.number_of_nodes()} nodes, {G.number_of_edges()} edges")

    edges = find_bridge_edges(G, cfg)

    info = []

    for u, v, k in edges:

        d = G.edges[u, v, k]

        info.append({"u": u, "v": v, "k": k, "length": d["length"], "name": str(d.get("name")),

                     "ref": str(d.get("ref")), "lat_u": G.nodes[u]["y"], "lon_u": G.nodes[u]["x"],