spa-l3-paris-emergency-routing

Rerun

"""
Emergency routing model for Paris area incidents.
Deliverables: incidents layer, closest-hospital routes, distance matrix, 15-min isochrones.
All outputs in EPSG:2154 (RGF93 / Lambert-93) in emergency_routing.gpkg.
"""
import osmnx as ox
import geopandas as gpd
import pandas as pd
import networkx as nx
import numpy as np
from shapely.geometry import Point, LineString, MultiPolygon
from shapely.ops import unary_union
import warnings

warnings.filterwarnings("ignore")
ox.settings.log_console = False

# ── 0. Config ──────────────────────────────────────────────────────────────
METRIC_CRS = "EPSG:2154"       # RGF93 / Lambert-93 — official metric for Paris
BBOX_N, BBOX_S, BBOX_E, BBOX_W = 48.88, 48.83, 2.38, 2.30
DEFAULT_SPEED_KMPH = 30        # assumed when maxspeed is missing
GPKG_PATH = "/work/emergency_routing.gpkg"
CSV_PATH = "/work/incidents.csv"

# ── 1. Load incidents ──────────────────────────────────────────────────────
inc_raw = pd.read_csv(CSV_PATH)
gdf_inc_wgs84 = gpd.GeoDataFrame(
    inc_raw[["incident_id"]],
    geometry=gpd.points_from_xy(inc_raw.longitude, inc_raw.latitude),
    crs="EPSG:4326",
)
print(f"Loaded {len(gdf_inc_wgs84)} incidents")

# ── 2. Drivable road network from OSM (via OSMnx / Overpass) ───────────────
print("Downloading drivable road network ...")
G = ox.graph_from_bbox(
    north=BBOX_N, south=BBOX_S, east=BBOX_E, west=BBOX_W,
    network_type="drive", simplify=True, retain_all=True,
)

# Add travel time (seconds) to every edge
for u, v, k, data in G.edges(keys=True, data=True):
    raw = data.get("maxspeed")
    if raw is None:
        speed_kmh = DEFAULT_SPEED_KMPH
    else:
        # maxspeed can be a string like "50" or a list ["50", "30"]
        if isinstance(raw, list):
            # take the first (most common) value; parse numeric
            speed_kmh = float(str(raw[0]).replace(" km/h", "").replace(" mph", "").strip())
        else:
            speed_kmh = float(str(raw).replace(" km/h", "").replace(" mph", "").strip())
    speed_ms = speed_kmh / 3.6
    length_m = data.get("length", 0.0)
    data["travel_time_s"] = length_m / speed_ms if speed_ms > 0 else 1.0

# Project graph to metric CRS
G_proj = ox.project_graph(G, to_crs=METRIC_CRS)
print(f"Graph: {len(G_proj.nodes)} nodes, {len(G_proj.edges)} edges")

# ── 3. Hospitals from OSM ─────────────────────────────────────────────────
print("Downloading hospitals ...")
hospitals_gdf = ox.features_from_bbox(
    north=BBOX_N, south=BBOX_S, east=BBOX_E, west=BBOX_W,
    tags={"amenity": "hospital"},
)

# Ensure name column exists
if "name" in hospitals_gdf.columns:
    hospitals_gdf["hospital_name"] = hospitals_gdf["name"].fillna("Unknown Hospital")
else:
    hospitals_gdf["hospital_name"] = "Unknown Hospital"

# Keep valid geometries, get centroids for polygons
hosp_pts = hospitals_gdf[hospitals_gdf.geometry.notna()].copy()
poly_mask = hosp_pts.geometry.type.isin(["Polygon", "MultiPolygon"])
if poly_mask.any():
    hosp_pts.loc[poly_mask, "geometry"] = hosp_pts.loc[poly_mask].centroid
# Keep only Points now
hosp_pts = hosp_pts[hosp_pts.geometry.type == "Point"].copy()
hosp_pts = hosp_pts[["hospital_name", "geometry"]].reset_index(drop=True)
hosp_gdf_wgs84 = gpd.GeoDataFrame(hosp_pts, geometry="geometry", crs="EPSG:4326")
print(f"Found {len(hosp_gdf_wgs84)} hospitals: {hosp_gdf_wgs84['hospital_name'].tolist()}")

# ── 4. Project incidents & hospitals to metric CRS ─────────────────────────
inc_gdf_metric = gdf_inc_wgs84.to_crs(METRIC_CRS)
hosp_gdf_metric = hosp_gdf_wgs84.to_crs(METRIC_CRS)

# Nearest graph node for each point
inc_gdf_metric["nn"] = ox.nearest_nodes(
    G_proj, inc_gdf_metric.geometry.x, inc_gdf_metric.geometry.y,
)
hosp_gdf_metric["nn"] = ox.nearest_nodes(
    G_proj, hosp_gdf_metric.geometry.x, hosp_gdf_metric.geometry.y,
)

# Pre-compute all-pairs shortest distances (fast with a single Dijkstra per incident)
print("Computing shortest paths ...")
inc_nn = inc_gdf_metric["nn"].values
hosp_nn = hosp_gdf_metric["nn"].values

# Build a lookup: node -> (hospital_name, hospital_index)
hosp_lookup = {
    row.nn: (row.hospital_name, idx)
    for idx, row in hosp_gdf_metric.iterrows()
}

# For each incident, find distances to every hospital
all_distances = []  # list of lists of (hosp_name, rank, dist_m)
for i_idx, inc_node in enumerate(inc_nn):
    inc_id = inc_gdf_metric.iloc[i_idx]["incident_id"]
    try:
        lengths = nx.single_source_dijkstra_path_length(
            G_proj, inc_node, cutoff=None, weight="length"
        )
    except Exception:
        lengths = {}

    dists = []
    for h_node, (h_name, h_idx) in hosp_lookup.items():
        d = lengths.get(h_node)
        if d is not None:
            dists.append((h_name, d))
    dists.sort(key=lambda x: x[1])
    all_distances.append((inc_id, dists))
    print(f"  {inc_id}: {len(dists)} reachable hospital(s)")

# ── 5. Closest hospital route (LineString geometry) ───────────────────────
print("Building closest-hospital routes ...")
closest_rows = []
for inc_id, dists in all_distances:
    if not dists:
        continue
    best_hosp_name, best_dist = dists[0]  # closest
    # Find the hospital node
    best_hosp_node = None
    for h_node, (h_name, h_idx) in hosp_lookup.items():
        if h_name == best_hosp_name:
            best_hosp_node = h_node
            break
    inc_node = inc_gdf_metric.loc[inc_gdf_metric["incident_id"] == inc_id, "nn"].iloc[0]
    try:
        path_nodes = nx.shortest_path(G_proj, inc_node, best_hosp_node, weight="length")
    except nx.NetworkXNoPath:
        continue

    # Build LineString geometry from edge geometries
    coords = []
    for i in range(len(path_nodes) - 1):
        u, v = path_nodes[i], path_nodes[i + 1]
        edge_data = G_proj.get_edge_data(u, v)
        if edge_data is None:
            continue
        # Use any key (parallel edges)
        first_key = list(edge_data.keys())[0]
        geom = edge_data[first_key].get("geometry")
        if geom is not None:
            pts = list(geom.coords)
            if i == 0:
                coords.extend(pts)
            else:
                coords.extend(pts[1:])
        else:
            # Straight line between nodes
            x1, y1 = G_proj.nodes[u]["x"], G_proj.nodes[u]["y"]
            x2, y2 = G_proj.nodes[v]["x"], G_proj.nodes[v]["y"]
            if i == 0:
                coords.append((x1, y1))
            coords.append((x2, y2))

    if len(coords) >= 2:
        line = LineString(coords)
    else:
        # Fallback: direct line (should not happen)
        p1 = Point(G_proj.nodes[inc_node]["x"], G_proj.nodes[inc_node]["y"])
        p2 = Point(G_proj.nodes[best_hosp_node]["x"], G_proj.nodes[best_hosp_node]["y"])
        line = LineString([p1, p2])

    closest_rows.append({
        "incident_id": inc_id,
        "hospital_name": best_hosp_name,
        "network_distance_m": round(best_dist, 1),
        "geometry": line,
    })

gdf_closest = gpd.GeoDataFrame(closest_rows, crs=METRIC_CRS)
print(f"  {len(gdf_closest)} routes built")

# ── 6. Distance matrix (3 closest hospitals per incident) ─────────────────
print("Building distance matrix ...")
matrix_rows = []
for inc_id, dists in all_distances:
    for rank, (h_name, d) in enumerate(dists[:3], 1):
        matrix_rows.append({
            "incident_id": inc_id,
            "hospital_name": h_name,
            "rank": rank,
            "network_distance_m": round(d, 1),
        })

gdf_matrix = gpd.GeoDataFrame(matrix_rows, geometry=[None] * len(matrix_rows), crs=METRIC_CRS)
print(f"  {len(gdf_matrix)} matrix entries")

# ── 7. 15-minute drive-time isochrones ────────────────────────────────────
print("Computing 15-min isochrones ...")
iso_rows = []
for h_idx, (_, hosp_row) in enumerate(hosp_gdf_metric.iterrows()):
    h_name = hosp_row["hospital_name"]
    h_node = hosp_row["nn"]

    # All nodes reachable within 900 seconds (15 min) along the road network
    try:
        reachable = nx.single_source_dijkstra_path_length(
            G_proj, h_node, cutoff=900.0, weight="travel_time_s"
        )
    except Exception:
        continue

    reachable_nodes = list(reachable.keys())
    if len(reachable_nodes) < 2:
        continue

    # Collect all edge LineStrings that are fully reachable
    subgraph = G_proj.subgraph(reachable_nodes)
    edge_lines = []
    for u, v, k, data in subgraph.edges(keys=True, data=True):
        geom = data.get("geometry")
        if geom is not None:
            edge_lines.append(geom)
        else:
            edge_lines.append(
                LineString([
                    (G_proj.nodes[u]["x"], G_proj.nodes[u]["y"]),
                    (G_proj.nodes[v]["x"], G_proj.nodes[v]["y"]),
                ])
            )

    if not edge_lines:
        continue

    # Union all edges and buffer by 150 m to create a smooth isochrone polygon
    merged = unary_union(edge_lines)
    isochrone = merged.buffer(150.0, cap_style=2, join_style=2)  # flat caps, round join

    if isochrone.geom_type == "Polygon":
        isochrone = MultiPolygon([isochrone])
    elif isochrone.geom_type != "MultiPolygon":
        # Shouldn't happen, skip
        continue

    iso_rows.append({
        "hospital_name": h_name,
        "travel_time_min": 15,
        "geometry": isochrone,
    })

gdf_iso = gpd.GeoDataFrame(iso_rows, crs=METRIC_CRS)
print(f"  {len(gdf_iso)} isochrones computed")

# ── 8. Write GPKG ─────────────────────────────────────────────────────────
print("Writing GPKG ...")
# Remove existing file if present
import os
if os.path.exists(GPKG_PATH):
    os.remove(GPKG_PATH)

# Layer: incidents — Point geometry, incident_id column
inc_out = inc_gdf_metric[["incident_id", "geometry"]].copy()
inc_out.to_file(GPKG_PATH, layer="incidents", driver="GPKG")
print(f"  Written layer 'incidents': {len(inc_out)} features")

# Layer: closest_hospital — LineString, columns incident_id + hospital_name + network_distance_m
gdf_closest.to_file(GPKG_PATH, layer="closest_hospital", driver="GPKG")
print(f"  Written layer 'closest_hospital': {len(gdf_closest)} features")

# Layer: distance_matrix — tabular (empty/null geometry), columns incident_id + hospital_name + rank + network_distance_m
gdf_matrix.to_file(GPKG_PATH, layer="distance_matrix", driver="GPKG")
print(f"  Written layer 'distance_matrix': {len(gdf_matrix)} features")

# New file mode: must append to existing GPKG; use mode='a' for gpkg
# Actually for subsequent layers we need mode='a' with fiona/geopandas
for layer_name, gdf in [("isochrones_15min", gdf_iso)]:
    gdf.to_file(GPKG_PATH, layer=layer_name, driver="GPKG")

print(f"  Written layer 'isochrones_15min': {len(gdf_iso)} features")

print("✓ All layers written to emergency_routing.gpkg")

"""
Emergency routing model for Paris area incidents using Overture Maps data.
Deliverables: incidents layer, closest-hospital routes, distance matrix, 15-min isochrones.
All outputs in EPSG:2154 (RGF93 / Lambert-93) in emergency_routing.gpkg.
"""
import geopandas as gpd
import pandas as pd
import numpy as np
import networkx as nx
from shapely.geometry import Point, LineString, MultiPolygon
from shapely.ops import unary_union
import overturemaps
import warnings
import os

warnings.filterwarnings("ignore")

# ── 0. Config ──────────────────────────────────────────────────────────────
METRIC_CRS = "EPSG:2154"          # RGF93 / Lambert-93 — official metric for Paris
DEFAULT_SPEED_KMPH = 30           # assumed when maxspeed is missing
GPKG_PATH = "/work/emergency_routing.gpkg"
CSV_PATH = "/work/incidents.csv"
BBOX = (2.30, 48.83, 2.38, 48.88)  # (xmin, ymin, xmax, ymax) for Overture
BBOX_W, BBOX_S, BBOX_E, BBOX_N = BBOX  # unpack for clarity

# ── 1. Load incidents ──────────────────────────────────────────────────────
inc_raw = pd.read_csv(CSV_PATH)
gdf_inc_wgs84 = gpd.GeoDataFrame(
    inc_raw[["incident_id"]],
    geometry=gpd.points_from_xy(inc_raw.longitude, inc_raw.latitude),
    crs="EPSG:4326",
)
print(f"✓ Loaded {len(gdf_inc_wgs84)} incidents")

# ── 2. Drivable road network from Overture ─────────────────────────────────
print("Downloading drivable road segments from Overture ...")
segments = overturemaps.geodataframe("segment", bbox=BBOX)
print(f"  Total segments: {len(segments)}")

# Filter to drivable roads
drivable_classes = [
    "motorway", "trunk", "primary", "secondary", "tertiary",
    "residential", "living_street", "service", "unclassified",
]
roads = segments[(segments["subtype"] == "road") & (segments["class"].isin(drivable_classes))].copy()
print(f"  Drivable road segments: {len(roads)}")

# Ensure metric CRS
roads = roads.to_crs(METRIC_CRS)

# Extract speed limit from each segment
def extract_maxspeed(speed_limits):
    """Extract max speed in km/h from Overture speed_limits field."""
    if speed_limits is None or (isinstance(speed_limits, list) and len(speed_limits) == 0):
        return DEFAULT_SPEED_KMPH
    if isinstance(speed_limits, list):
        # Take the first speed limit entry
        for sl in speed_limits:
            if sl and sl.get("max_speed") and sl["max_speed"].get("value"):
                val = sl["max_speed"]["value"]
                unit = sl["max_speed"].get("unit", "km/h")
                if unit == "mph":
                    return val * 1.60934
                return float(val)
    return DEFAULT_SPEED_KMPH

roads["maxspeed_kmh"] = roads["speed_limits"].apply(extract_maxspeed)
roads["travel_time_s"] = (roads.geometry.length / (roads["maxspeed_kmh"] / 3.6)).round(1)

# ── 3. Build routable graph from road segments ────────────────────────────
print("Building routable graph from road segments ...")
G = nx.MultiDiGraph()

# Track edges by their geometry coordinates to avoid duplicates
edge_key_counter = {}

for idx, row in roads.iterrows():
    geom = row.geometry
    if geom is None or geom.is_empty:
        continue
    coords = list(geom.coords)
    if len(coords) < 2:
        continue

    speed_kmh = row["maxspeed_kmh"]
    length_m = geom.length
    tt_s = row["travel_time_s"]

    # Add nodes
    for i, (x, y) in enumerate(coords):
        node_id = (round(x, 6), round(y, 6))  # approximate matching
        if node_id not in G:
            G.add_node(node_id, x=x, y=y)

    # Add edges between consecutive coordinates
    for i in range(len(coords) - 1):
        u = (round(coords[i][0], 6), round(coords[i][1], 6))
        v = (round(coords[i + 1][0], 6), round(coords[i + 1][1], 6))
        if u == v:
            continue
        # Edge geometry (subsegment)
        edge_geom = LineString([coords[i], coords[i + 1]])
        sub_length = edge_geom.length
        sub_tt = sub_length / (speed_kmh / 3.6) if speed_kmh > 0 else 1.0

        # Use a simple key
        G.add_edge(u, v, key=0, length=sub_length, travel_time_s=sub_tt, geometry=edge_geom)
        G.add_edge(v, u, key=0, length=sub_length, travel_time_s=sub_tt, geometry=edge_geom)

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

# ── 4. Hospitals from Overture ───────────────────────────────────────────
print("Downloading hospitals from Overture ...")
places = overturemaps.geodataframe("place", bbox=BBOX)
# Filter to hospital category
hospitals = places[places["categories"].apply(
    lambda c: isinstance(c, dict) and c.get("primary") == "hospital"
)].copy()
print(f"  Found {len(hospitals)} hospitals in Overture places")

if len(hospitals) == 0:
    # Fallback: try taxonomy
    hospitals = places[places["taxonomy"].apply(
        lambda c: isinstance(c, dict) and c.get("primary") == "hospital"
    )].copy()
    print(f"  (from taxonomy: {len(hospitals)})")

# Get names
def get_name(names_field):
    if isinstance(names_field, dict):
        return names_field.get("primary", "Unknown Hospital")
    return "Unknown Hospital"

hospitals["hospital_name"] = hospitals["names"].apply(get_name)
hospitals = hospitals[["hospital_name", "geometry"]].reset_index(drop=True)
hospitals = hospitals.to_crs("EPSG:4326")

# Use centroids for polygon hospital grounds
poly_mask = hospitals.geometry.type.isin(["Polygon", "MultiPolygon"])
if poly_mask.any():
    hospitals.loc[poly_mask, "geometry"] = hospitals.loc[poly_mask].centroid
hospitals = hospitals[hospitals.geometry.type == "Point"].copy()
hospitals = hospitals.reset_index(drop=True)
print(f"  {len(hospitals)} hospital point locations")

if len(hospitals) == 0:
    print("  WARNING: No hospitals found! Using dummy data.")
    # Fallback: create hospitals at known Paris locations
    hospitals = gpd.GeoDataFrame({
        "hospital_name": ["Hôpital Hôtel-Dieu", "Hôpital Lariboisière", "Hôpital Cochin"],
        "geometry": [Point(2.3481, 48.8542), Point(2.3510, 48.8820), Point(2.3390, 48.8390)]
    }, crs="EPSG:4326")
    print(f"  Using {len(hospitals)} fallback hospitals")

# ── 5. Project to metric CRS & get nearest nodes ──────────────────────────
inc_metric = gdf_inc_wgs84.to_crs(METRIC_CRS)
hosp_metric = hospitals.to_crs(METRIC_CRS)

def nearest_graph_node(gdf_metric, G):
    """Find nearest graph node for each point."""
    nodes_gdf = gpd.GeoDataFrame(
        [{"node": n, "geometry": Point(data["x"], data["y"])}
         for n, data in G.nodes(data=True)],
        crs=METRIC_CRS
    )
    # For each point, find the nearest node
    nearest_nodes_list = []
    for _, pt in gdf_metric.iterrows():
        dists = nodes_gdf.distance(pt.geometry)
        nearest_nodes_list.append(nodes_gdf.iloc[dists.idxmin()]["node"])
    return nearest_nodes_list

inc_metric["nn"] = nearest_graph_node(inc_metric, G)
hosp_metric["nn"] = nearest_graph_node(hosp_metric, G)

# Build hospital node lookup
hosp_lookup = {}
for _, row in hosp_metric.iterrows():
    node = row["nn"]
    if node not in hosp_lookup:
        hosp_lookup[node] = row["hospital_name"]
    # If multiple hospitals at same node, keep both

# ── 6. Shortest paths (distance by length) ────────────────────────────────
print("Computing shortest paths ...")

# For each incident, compute distances to all hospitals
all_distances = []  # list of (inc_id, [(hosp_name, dist_m), ...])

for i_idx, inc_row in inc_metric.iterrows():
    inc_id = inc_row["incident_id"]
    inc_node = inc_row["nn"]

    try:
        lengths = nx.single_source_dijkstra_path_length(
            G, inc_node, cutoff=None, weight="length"
        )
    except nx.NetworkXNoPath:
        lengths = {}

    dists = []
    for h_node, h_name in hosp_lookup.items():
        d = lengths.get(h_node)
        if d is not None:
            dists.append((h_name, d))
    dists.sort(key=lambda x: x[1])
    all_distances.append((inc_id, dists))
    print(f"  {inc_id}: {len(dists)} reachable hospital(s)")

# ── 7. Closest hospital route ─────────────────────────────────────────────
print("Building closest-hospital routes ...")
closest_rows = []

for inc_id, dists in all_distances:
    if not dists:
        continue
    best_hosp, best_dist = dists[0]

    inc_node = inc_metric.loc[inc_metric["incident_id"] == inc_id, "nn"].iloc[0]
    # Find hospital node for the closest hospital
    best_h_node = None
    for h_node, h_name in hosp_lookup.items():
        if h_name == best_hosp:
            best_h_node = h_node
            break
    if best_h_node is None:
        continue

    try:
        path_nodes = nx.shortest_path(G, inc_node, best_h_node, weight="length")
    except nx.NetworkXNoPath:
        continue

    # Build LineString from edge geometries
    coords = []
    for i in range(len(path_nodes) - 1):
        u, v = path_nodes[i], path_nodes[i + 1]
        edge_data = G.get_edge_data(u, v)
        if edge_data is None:
            # Use direct line
            x1, y1 = G.nodes[u]["x"], G.nodes[u]["y"]
            x2, y2 = G.nodes[v]["x"], G.nodes[v]["y"]
            if i == 0:
                coords.append((x1, y1))
            coords.append((x2, y2))
            continue

        # Try first edge key
        first_key = list(edge_data.keys())[0]
        geom = edge_data[first_key].get("geometry")
        if geom is not None:
            pts = list(geom.coords)
            if i == 0:
                coords.extend(pts)
            else:
                coords.extend(pts[1:])
        else:
            x1, y1 = G.nodes[u]["x"], G.nodes[u]["y"]
            x2, y2 = G.nodes[v]["x"], G.nodes[v]["y"]
            if i == 0:
                coords.append((x1, y1))
            coords.append((x2, y2))

    if len(coords) >= 2:
        line = LineString(coords)
    else:
        p1 = Point(G.nodes[inc_node]["x"], G.nodes[inc_node]["y"])
        p2 = Point(G.nodes[best_h_node]["x"], G.nodes[best_h_node]["y"])
        line = LineString([p1, p2])

    closest_rows.append({
        "incident_id": inc_id,
        "hospital_name": best_hosp,
        "network_distance_m": round(best_dist, 1),
        "geometry": line,
    })

gdf_closest = gpd.GeoDataFrame(closest_rows, crs=METRIC_CRS)
print(f"  {len(gdf_closest)} closest routes built")

# ── 8. Distance matrix ────────────────────────────────────────────────────
print("Building distance matrix ...")
matrix_rows = []
for inc_id, dists in all_distances:
    for rank, (h_name, d) in enumerate(dists[:3], 1):
        matrix_rows.append({
            "incident_id": inc_id,
            "hospital_name": h_name,
            "rank": rank,
            "network_distance_m": round(d, 1),
        })

gdf_matrix = gpd.GeoDataFrame(matrix_rows, geometry=[None] * len(matrix_rows), crs=METRIC_CRS)
print(f"  {len(gdf_matrix)} matrix entries")

# ── 9. 15-minute drive-time isochrones ────────────────────────────────────
print("Computing 15-minute isochrones ...")
iso_rows = []

for _, hosp_row in hosp_metric.iterrows():
    h_name = hosp_row["hospital_name"]
    h_node = hosp_row["nn"]

    # All nodes reachable within 900 seconds (15 min) using travel_time_s
    try:
        reachable = nx.single_source_dijkstra_path_length(
            G, h_node, cutoff=900.0, weight="travel_time_s"
        )
    except nx.NetworkXNoPath:
        continue

    reachable_nodes = set(reachable.keys())
    if len(reachable_nodes) < 2:
        continue

    # Collect edge lines that are reachable (both endpoints reachable)
    edge_lines = []
    for u, v, k, data in G.edges(keys=True, data=True):
        if u in reachable_nodes and v in reachable_nodes:
            geom = data.get("geometry")
            if geom is not None:
                edge_lines.append(geom)
            else:
                edge_lines.append(
                    LineString([
                        (G.nodes[u]["x"], G.nodes[u]["y"]),
                        (G.nodes[v]["x"], G.nodes[v]["y"]),
                    ])
                )

    if not edge_lines:
        continue

    merged = unary_union(edge_lines)
    isochrone = merged.buffer(150.0, cap_style=2, join_style=2)

    if isochrone.geom_type == "Polygon":
        isochrone = MultiPolygon([isochrone])
    elif isochrone.geom_type != "MultiPolygon":
        continue

    iso_rows.append({
        "hospital_name": h_name,
        "travel_time_min": 15,
        "geometry": isochrone,
    })

gdf_iso = gpd.GeoDataFrame(iso_rows, crs=METRIC_CRS)
print(f"  {len(gdf_iso)} isochrones computed")

# ── 10. Write GPKG ────────────────────────────────────────────────────────
print("Writing GPKG ...")
if os.path.exists(GPKG_PATH):
    os.remove(GPKG_PATH)

# Layer: incidents
inc_out = inc_metric[["incident_id", "geometry"]].copy()
inc_out.to_file(GPKG_PATH, layer="incidents", driver="GPKG")
print(f"  Layer 'incidents': {len(inc_out)} features")

# Layer: closest_hospital
gdf_closest.to_file(GPKG_PATH, layer="closest_hospital", driver="GPKG")
print(f"  Layer 'closest_hospital': {len(gdf_closest)} features")

# Layer: distance_matrix
gdf_matrix.to_file(GPKG_PATH, layer="distance_matrix", driver="GPKG")
print(f"  Layer 'distance_matrix': {len(gdf_matrix)} features")

# Layer: isochrones_15min
gdf_iso.to_file(GPKG_PATH, layer="isochrones_15min", driver="GPKG")
print(f"  Layer 'isochrones_15min': {len(gdf_iso)} features")

print("\n✓ All layers written to emergency_routing.gpkg")

"""
Emergency routing model for Paris area incidents using Overture Maps data.
Deliverables: incidents layer, closest-hospital routes, distance matrix, 15-min isochrones.
All outputs in EPSG:2154 (RGF93 / Lambert-93) in emergency_routing.gpkg.
"""
import geopandas as gpd
import pandas as pd
import numpy as np
import networkx as nx
from shapely.geometry import Point, LineString, MultiPolygon, box
from shapely.ops import unary_union
from shapely import STRtree
import overturemaps
import warnings
import os

warnings.filterwarnings("ignore")

# ── 0. Config ──────────────────────────────────────────────────────────────
METRIC_CRS = "EPSG:2154"
DEFAULT_SPEED_KMPH = 30
GPKG_PATH = "/work/emergency_routing.gpkg"
CSV_PATH = "/work/incidents.csv"
BBOX = (2.30, 48.83, 2.38, 48.88)   # (xmin, ymin, xmax, ymax)

# ── 1. Load incidents ──────────────────────────────────────────────────────
inc_raw = pd.read_csv(CSV_PATH)
gdf_inc_wgs84 = gpd.GeoDataFrame(
    inc_raw[["incident_id"]],
    geometry=gpd.points_from_xy(inc_raw.longitude, inc_raw.latitude),
    crs="EPSG:4326",
)
print(f"✓ Loaded {len(gdf_inc_wgs84)} incidents")
inc_metric = gdf_inc_wgs84.to_crs(METRIC_CRS)

# ── 2. Drivable road network from Overture ─────────────────────────────────
print("Downloading drivable road segments from Overture ...")
segments = overturemaps.geodataframe("segment", bbox=BBOX)
print(f"  Total segments: {len(segments)}")

drivable_classes = [
    "motorway", "trunk", "primary", "secondary", "tertiary",
    "residential", "living_street", "service", "unclassified",
]
roads = segments[(segments["subtype"] == "road") & (segments["class"].isin(drivable_classes))].copy()
print(f"  Drivable road segments: {len(roads)}")

if roads.crs is None:
    roads.set_crs("EPSG:4326", inplace=True)
roads = roads.to_crs(METRIC_CRS)

# Extract max speed
def extract_maxspeed(speed_limits):
    if speed_limits is None or not isinstance(speed_limits, list) or len(speed_limits) == 0:
        return DEFAULT_SPEED_KMPH
    for sl in speed_limits:
        if sl and sl.get("max_speed") and sl["max_speed"].get("value"):
            val = sl["max_speed"]["value"]
            unit = sl["max_speed"].get("unit", "km/h")
            if unit == "mph":
                return val * 1.60934
            return float(val)
    return DEFAULT_SPEED_KMPH

roads["maxspeed_kmh"] = roads["speed_limits"].apply(extract_maxspeed)
roads["travel_time_s"] = (roads.geometry.length / (roads["maxspeed_kmh"] / 3.6)).round(1)

# ── 3. Build routable graph (simplified: endpoints as nodes) ──────────────
# Build a graph where each road segment is an edge between its two endpoints
# Use connector IDs as node identifiers for connecting segments at junctions
print("Building routable graph ...")
G = nx.MultiDiGraph()

# Map: connector_id -> (x, y) from first occurrence
connector_coords = {}
# For edges: (u_connector_id, v_connector_id) -> list of edge data
edge_list = []

for idx, row in roads.iterrows():
    geom = row.geometry
    if geom is None or geom.is_empty:
        continue
    coords = list(geom.coords)
    if len(coords) < 2:
        continue

    # Get connectors
    connectors = row.get("connectors")
    if connectors and isinstance(connectors, list) and len(connectors) >= 2:
        # Use connector IDs as nodes
        try:
            c_start = connectors[0]  # at=0.0
            c_end = connectors[1]    # at=1.0
            if isinstance(c_start, dict) and isinstance(c_end, dict):
                u_id = c_start["connector_id"]
                v_id = c_end["connector_id"]
            elif isinstance(c_start, tuple) or isinstance(c_start, list):
                u_id = c_start[0]
                v_id = c_start[1]
            else:
                continue
        except (IndexError, KeyError, TypeError):
            continue
    else:
        # Fall back to coordinate-based nodes at endpoints
        u_id = (round(coords[0][0], 3), round(coords[0][1], 3))
        v_id = (round(coords[-1][0], 3), round(coords[-1][1], 3))

    # Store coordinates for nodes
    geom_start = Point(coords[0])
    geom_end = Point(coords[-1])
    connector_coords[u_id] = (geom_start.x, geom_start.y)
    connector_coords[v_id] = (geom_end.x, geom_end.y)

    # Add directed edges (both directions for drivable roads)
    length_m = geom.length
    speed_kmh = row["maxspeed_kmh"]
    tt_s = row["travel_time_s"]

    edge_list.append((u_id, v_id, length_m, tt_s, geom))

print(f"  {len(connector_coords)} unique nodes, {len(edge_list)} segments")

# Add all nodes and edges to the graph
for node_id, (x, y) in connector_coords.items():
    G.add_node(node_id, x=x, y=y)

for u_id, v_id, length_m, tt_s, geom in edge_list:
    if u_id in G and v_id in G:
        G.add_edge(u_id, v_id, key=0, length=length_m, travel_time_s=tt_s, geometry=geom)
        G.add_edge(v_id, u_id, key=0, length=length_m, travel_time_s=tt_s, geometry=geom)

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

# ── 4. Hospitals from Overture ───────────────────────────────────────────
print("Downloading hospitals from Overture ...")
places = overturemaps.geodataframe("place", bbox=BBOX)

# Filter to hospitals
hospitals = places[places["categories"].apply(
    lambda c: isinstance(c, dict) and c.get("primary") == "hospital"
)].copy()

if len(hospitals) == 0:
    hospitals = places[places["taxonomy"].apply(
        lambda c: isinstance(c, dict) and c.get("primary") == "hospital"
    )].copy()

def get_name(names_field):
    if isinstance(names_field, dict):
        return names_field.get("primary", "Unknown Hospital")
    return "Unknown Hospital"

hospitals["hospital_name"] = hospitals["names"].apply(get_name)
hospitals = hospitals[["hospital_name", "geometry"]].reset_index(drop=True)
if hospitals.crs is None:
    hospitals.set_crs("EPSG:4326", inplace=True)
hospitals = hospitals.to_crs(METRIC_CRS)

# Use centroids for polygons
poly_mask = hospitals.geometry.type.isin(["Polygon", "MultiPolygon"])
if poly_mask.any():
    hospitals.loc[poly_mask, "geometry"] = hospitals.loc[poly_mask].centroid
hospitals = hospitals[hospitals.geometry.type == "Point"].copy()
hospitals = hospitals.reset_index(drop=True)
print(f"  {len(hospitals)} hospitals")

# ── 5. Find nearest graph node for each POI using STRtree ────────────────
print("Finding nearest graph nodes ...")
node_ids = list(G.nodes(data=False))
node_coords = np.array([(G.nodes[n]["x"], G.nodes[n]["y"]) for n in node_ids])
node_points = [Point(xy) for xy in node_coords]

# Build spatial index
tree = STRtree(node_points)

def nearest_node(point_geom):
    """Return nearest graph node id for a point."""
    # Query the tree
    idx = tree.nearest(point_geom)
    return node_ids[idx]

inc_metric["nn"] = inc_metric.geometry.apply(nearest_node)
hospitals["nn"] = hospitals.geometry.apply(nearest_node)

# Build hospital lookup: mapping from node -> list of (hospital_name, index)
hosp_lookup = {}
for idx, row in hospitals.iterrows():
    node = row["nn"]
    if node not in hosp_lookup:
        hosp_lookup[node] = []
    hosp_lookup[node].append(row["hospital_name"])

# ── 6. Shortest paths ─────────────────────────────────────────────────────
print("Computing shortest paths ...")
all_distances = []

for _, inc_row in inc_metric.iterrows():
    inc_id = inc_row["incident_id"]
    inc_node = inc_row["nn"]

    try:
        lengths = nx.single_source_dijkstra_path_length(
            G, inc_node, cutoff=None, weight="length"
        )
    except Exception:
        lengths = {}

    dists = []
    for h_node, h_names in hosp_lookup.items():
        d = lengths.get(h_node)
        if d is not None:
            for h_name in h_names:
                dists.append((h_name, d))
    dists.sort(key=lambda x: x[1])
    all_distances.append((inc_id, dists))
    print(f"  {inc_id}: {len(dists)} reachable hospital(s)")

# ── 7. Closest hospital routes ────────────────────────────────────────────
print("Building closest-hospital routes ...")
closest_rows = []

for inc_id, dists in all_distances:
    if not dists:
        continue
    best_hosp, best_dist = dists[0]

    inc_node = inc_metric.loc[inc_metric["incident_id"] == inc_id, "nn"].iloc[0]
    # Find hospital node for the closest hospital
    best_h_node = None
    for h_node, h_names in hosp_lookup.items():
        if best_hosp in h_names:
            best_h_node = h_node
            break
    if best_h_node is None:
        continue

    try:
        path_nodes = nx.shortest_path(G, inc_node, best_h_node, weight="length")
    except nx.NetworkXNoPath:
        continue

    # Build LineString from edge geometries
    coords = []
    for i in range(len(path_nodes) - 1):
        u, v = path_nodes[i], path_nodes[i + 1]
        edge_data = G.get_edge_data(u, v)
        if edge_data is None:
            x1, y1 = G.nodes[u]["x"], G.nodes[u]["y"]
            x2, y2 = G.nodes[v]["x"], G.nodes[v]["y"]
            if i == 0:
                coords.append((x1, y1))
            coords.append((x2, y2))
            continue
        first_key = list(edge_data.keys())[0]
        geom = edge_data[first_key].get("geometry")
        if geom is not None:
            pts = list(geom.coords)
            if i == 0:
                coords.extend(pts)
            else:
                coords.extend(pts[1:])
        else:
            x1, y1 = G.nodes[u]["x"], G.nodes[u]["y"]
            x2, y2 = G.nodes[v]["x"], G.nodes[v]["y"]
            if i == 0:
                coords.append((x1, y1))
            coords.append((x2, y2))

    if len(coords) >= 2:
        line = LineString(coords)
    else:
        p1 = Point(G.nodes[inc_node]["x"], G.nodes[inc_node]["y"])
        p2 = Point(G.nodes[best_h_node]["x"], G.nodes[best_h_node]["y"])
        line = LineString([p1, p2])

    closest_rows.append({
        "incident_id": inc_id,
        "hospital_name": best_hosp,
        "network_distance_m": round(best_dist, 1),
        "geometry": line,
    })

gdf_closest = gpd.GeoDataFrame(closest_rows, crs=METRIC_CRS)
print(f"  {len(gdf_closest)} routes")

# ── 8. Distance matrix ────────────────────────────────────────────────────
print("Building distance matrix ...")
matrix_rows = []
for inc_id, dists in all_distances:
    for rank, (h_name, d) in enumerate(dists[:3], 1):
        matrix_rows.append({
            "incident_id": inc_id,
            "hospital_name": h_name,
            "rank": rank,
            "network_distance_m": round(d, 1),
        })

gdf_matrix = gpd.GeoDataFrame(matrix_rows, geometry=[None] * len(matrix_rows), crs=METRIC_CRS)
print(f"  {len(gdf_matrix)} entries")

# ── 9. 15-minute isochrones ───────────────────────────────────────────────
print("Computing 15-min isochrones ...")
iso_rows = []

for _, hosp_row in hospitals.iterrows():
    h_name = hosp_row["hospital_name"]
    h_node = hosp_row["nn"]

    try:
        reachable = nx.single_source_dijkstra_path_length(
            G, h_node, cutoff=900.0, weight="travel_time_s"
        )
    except Exception:
        continue

    reachable_nodes = set(reachable.keys())
    if len(reachable_nodes) < 2:
        continue

    edge_lines = []
    for u, v, k, data in G.edges(keys=True, data=True):
        if u in reachable_nodes and v in reachable_nodes:
            geom = data.get("geometry")
            if geom is not None:
                edge_lines.append(geom)
            else:
                edge_lines.append(
                    LineString([
                        (G.nodes[u]["x"], G.nodes[u]["y"]),
                        (G.nodes[v]["x"], G.nodes[v]["y"]),
                    ])
                )

    if not edge_lines:
        continue

    merged = unary_union(edge_lines)
    isochrone = merged.buffer(150.0, cap_style=2, join_style=2)

    if isochrone.geom_type == "Polygon":
        isochrone = MultiPolygon([isochrone])
    elif isochrone.geom_type != "MultiPolygon":
        continue

    iso_rows.append({
        "hospital_name": h_name,
        "travel_time_min": 15,
        "geometry": isochrone,
    })

gdf_iso = gpd.GeoDataFrame(iso_rows, crs=METRIC_CRS)
print(f"  {len(gdf_iso)} isochrones")

# ── 10. Write GPKG ────────────────────────────────────────────────────────
print("Writing GPKG ...")
if os.path.exists(GPKG_PATH):
    os.remove(GPKG_PATH)

inc_out = inc_metric[["incident_id", "geometry"]].copy()
inc_out.to_file(GPKG_PATH, layer="incidents", driver="GPKG")
print(f"  Layer 'incidents': {len(inc_out)} features")

gdf_closest.to_file(GPKG_PATH, layer="closest_hospital", driver="GPKG")
print(f"  Layer 'closest_hospital': {len(gdf_closest)} features")

gdf_matrix.to_file(GPKG_PATH, layer="distance_matrix", driver="GPKG")
print(f"  Layer 'distance_matrix': {len(gdf_matrix)} features")

gdf_iso.to_file(GPKG_PATH, layer="isochrones_15min", driver="GPKG")
print(f"  Layer 'isochrones_15min': {len(gdf_iso)} features")

print("\n✓ All layers written to emergency_routing.gpkg")

"""
Emergency routing model for Paris area incidents using Overture Maps data.
Deliverables: incidents layer, closest-hospital routes, distance matrix, 15-min isochrones.
All outputs in EPSG:2154 (RGF93 / Lambert-93) in emergency_routing.gpkg.
"""
import geopandas as gpd
import pandas as pd
import numpy as np
import networkx as nx
from shapely.geometry import Point, LineString, MultiPolygon
from shapely.ops import unary_union
from shapely import STRtree
import overturemaps
import warnings
import os

warnings.filterwarnings("ignore")

# ── 0. Config ──────────────────────────────────────────────────────────────
METRIC_CRS = "EPSG:2154"
DEFAULT_SPEED_KMPH = 30
GPKG_PATH = "/work/emergency_routing.gpkg"
CSV_PATH = "/work/incidents.csv"
BBOX = (2.30, 48.83, 2.38, 48.88)   # (xmin, ymin, xmax, ymax)

# ── 1. Load incidents ──────────────────────────────────────────────────────
inc_raw = pd.read_csv(CSV_PATH)
gdf_inc_wgs84 = gpd.GeoDataFrame(
    inc_raw[["incident_id"]],
    geometry=gpd.points_from_xy(inc_raw.longitude, inc_raw.latitude),
    crs="EPSG:4326",
)
print(f"✓ Loaded {len(gdf_inc_wgs84)} incidents")
inc_metric = gdf_inc_wgs84.to_crs(METRIC_CRS)

# ── 2. Drivable road network from Overture ─────────────────────────────────
print("Downloading drivable road segments from Overture ...")
segments = overturemaps.geodataframe("segment", bbox=BBOX)
print(f"  Total segments: {len(segments)}")

drivable_classes = [
    "motorway", "trunk", "primary", "secondary", "tertiary",
    "residential", "living_street", "service", "unclassified",
]
roads = segments[(segments["subtype"] == "road") & (segments["class"].isin(drivable_classes))].copy()
print(f"  Drivable road segments: {len(roads)}")

if roads.crs is None:
    roads.set_crs("EPSG:4326", inplace=True)
roads = roads.to_crs(METRIC_CRS)

# Extract max speed
def extract_maxspeed(speed_limits):
    if speed_limits is None or not isinstance(speed_limits, list) or len(speed_limits) == 0:
        return DEFAULT_SPEED_KMPH
    for sl in speed_limits:
        if sl and sl.get("max_speed") and sl["max_speed"].get("value"):
            val = sl["max_speed"]["value"]
            unit = sl["max_speed"].get("unit", "km/h")
            if unit == "mph":
                return val * 1.60934
            return float(val)
    return DEFAULT_SPEED_KMPH

roads["maxspeed_kmh"] = roads["speed_limits"].apply(extract_maxspeed)
roads["travel_time_s"] = (roads.geometry.length / (roads["maxspeed_kmh"] / 3.6)).round(1)

# ── 3. Build routable graph ──────────────────────────────────────────────
# Use segment endpoint coordinates (hashed) as node IDs for robustness.
# When connectors are available, use connector_id as node identifier.
print("Building routable graph ...")
G = nx.MultiDiGraph()

# Node registry: node_id -> (x, y)
node_coords = {}
# For nearest-node lookup: list of (node_id, Point geometry)
node_geom_list = []

for idx, row in roads.iterrows():
    geom = row.geometry
    if geom is None or geom.is_empty:
        continue
    coords = list(geom.coords)
    if len(coords) < 2:
        continue

    # Get connectors for node identification
    connectors = row.get("connectors")
    has_good_connectors = (
        connectors is not None
        and not isinstance(connectors, (float, np.floating))
        and len(connectors) >= 2
    )

    if has_good_connectors:
        try:
            c0 = connectors[0]
            c1 = connectors[-1]
            if isinstance(c0, dict):
                u_id = c0.get("connector_id", str(coords[0]))
                v_id = c1.get("connector_id", str(coords[-1]))
            elif isinstance(c0, (list, tuple)):
                u_id = c0[0]
                v_id = c1[0]
            else:
                u_id = f"c_{coords[0][0]:.3f}_{coords[0][1]:.3f}"
                v_id = f"c_{coords[-1][0]:.3f}_{coords[-1][1]:.3f}"
        except (IndexError, KeyError, TypeError):
            u_id = f"c_{coords[0][0]:.3f}_{coords[0][1]:.3f}"
            v_id = f"c_{coords[-1][0]:.3f}_{coords[-1][1]:.3f}"
    else:
        u_id = f"c_{coords[0][0]:.3f}_{coords[0][1]:.3f}"
        v_id = f"c_{coords[-1][0]:.3f}_{coords[-1][1]:.3f}"

    # Store/retrieve coordinates
    p_start = Point(coords[0])
    p_end = Point(coords[-1])
    if u_id not in node_coords:
        node_coords[u_id] = (p_start.x, p_start.y)
        node_geom_list.append((u_id, p_start))
    if v_id not in node_coords:
        node_coords[v_id] = (p_end.x, p_end.y)
        node_geom_list.append((v_id, p_end))

    # Add nodes
    for nid, (x, y) in [(u_id, (p_start.x, p_start.y)), (v_id, (p_end.x, p_end.y))]:
        if nid not in G:
            G.add_node(nid, x=x, y=y)

    length_m = geom.length
    speed_kmh = row["maxspeed_kmh"]
    tt_s = row["travel_time_s"]

    # Directed edges
    G.add_edge(u_id, v_id, key=0, length=length_m, travel_time_s=tt_s, geometry=geom)
    G.add_edge(v_id, u_id, key=0, length=length_m, travel_time_s=tt_s, geometry=geom)

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

# Build spatial index for nearest-node queries
node_ids_arr = [item[0] for item in node_geom_list]
node_pts_arr = [item[1] for item in node_geom_list]
node_tree = STRtree(node_pts_arr)

# Map tree index back to node id
# The tree returns the geometry; we need to map back. Let's use a simpler approach.
# STRtree.nearest(geom) returns the index of the geometry in the input list
# Wait, the test showed it returns the geometry itself? Let me check.
# Actually from the test: tree = STRtree(pts) where pts=[Point(0,0), Point(1,1), Point(2,2)]
# tree.nearest(Point(0.1, 0.1)) returned 0 (the index). So it returns an integer index.

def nearest_node_id(point_geom):
    """Find the graph node id closest to the given point."""
    idx = int(node_tree.nearest(point_geom))
    return node_ids_arr[idx]

# ── 4. Hospitals from Overture ───────────────────────────────────────────
print("Downloading hospitals from Overture ...")
places = overturemaps.geodataframe("place", bbox=BBOX)

hospitals = places[places["categories"].apply(
    lambda c: isinstance(c, dict) and c.get("primary") == "hospital"
)].copy()

if len(hospitals) == 0:
    hospitals = places[places["taxonomy"].apply(
        lambda c: isinstance(c, dict) and c.get("primary") == "hospital"
    )].copy()

def get_name(names_field):
    if isinstance(names_field, dict):
        return names_field.get("primary", "Unknown Hospital")
    return "Unknown Hospital"

hospitals["hospital_name"] = hospitals["names"].apply(get_name)
hospitals = hospitals[["hospital_name", "geometry"]].reset_index(drop=True)
if hospitals.crs is None:
    hospitals.set_crs("EPSG:4326", inplace=True)
hospitals = hospitals.to_crs(METRIC_CRS)

poly_mask = hospitals.geometry.type.isin(["Polygon", "MultiPolygon"])
if poly_mask.any():
    hospitals.loc[poly_mask, "geometry"] = hospitals.loc[poly_mask].centroid
hospitals = hospitals[hospitals.geometry.type == "Point"].copy()
hospitals = hospitals.reset_index(drop=True)
print(f"  {len(hospitals)} hospitals")

# ── 5. Nearest graph nodes ────────────────────────────────────────────────
print("Finding nearest graph nodes ...")
inc_metric["nn"] = inc_metric.geometry.apply(nearest_node_id)
hospitals["nn"] = hospitals.geometry.apply(nearest_node_id)

# Build hospital lookup
hosp_lookup = {}
for _, row in hospitals.iterrows():
    node = row["nn"]
    if node not in hosp_lookup:
        hosp_lookup[node] = []
    hosp_lookup[node].append(row["hospital_name"])

print(f"  Incidents linked to {inc_metric['nn'].nunique()} unique graph nodes")
print(f"  Hospitals linked to {hospitals['nn'].nunique()} unique graph nodes")

# Verify nodes exist in graph
for inc_node in inc_metric["nn"]:
    if inc_node not in G:
        print(f"  WARNING: incident node {inc_node} not in graph!")
for h_node in hospitals["nn"]:
    if h_node not in G:
        print(f"  WARNING: hospital node {h_node} not in graph!")

# ── 6. Shortest paths ─────────────────────────────────────────────────────
print("Computing shortest paths ...")
all_distances = []

for _, inc_row in inc_metric.iterrows():
    inc_id = inc_row["incident_id"]
    inc_node = inc_row["nn"]

    try:
        lengths = nx.single_source_dijkstra_path_length(
            G, inc_node, cutoff=None, weight="length"
        )
    except Exception:
        lengths = {}

    dists = []
    for h_node, h_names in hosp_lookup.items():
        d = lengths.get(h_node)
        if d is not None:
            for h_name in h_names:
                dists.append((h_name, d))
    dists.sort(key=lambda x: x[1])
    all_distances.append((inc_id, dists))
    n_reach = len(dists)
    print(f"  {inc_id}: {n_reach} reachable hospital(s)")
    if n_reach == 0:
        # Debug: check if the incident node is isolated
        print(f"    degree of node {inc_node}: {G.degree(inc_node)}")

# ── 7. Closest hospital routes ────────────────────────────────────────────
print("Building closest-hospital routes ...")
closest_rows = []

for inc_id, dists in all_distances:
    if not dists:
        continue
    best_hosp, best_dist = dists[0]
    inc_node = inc_metric.loc[inc_metric["incident_id"] == inc_id, "nn"].iloc[0]

    best_h_node = None
    for h_node, h_names in hosp_lookup.items():
        if best_hosp in h_names:
            best_h_node = h_node
            break
    if best_h_node is None:
        continue

    try:
        path_nodes = nx.shortest_path(G, inc_node, best_h_node, weight="length")
    except nx.NetworkXNoPath:
        continue

    coords = []
    for i in range(len(path_nodes) - 1):
        u, v = path_nodes[i], path_nodes[i + 1]
        edge_data = G.get_edge_data(u, v)
        if edge_data is None:
            x1, y1 = G.nodes[u]["x"], G.nodes[u]["y"]
            x2, y2 = G.nodes[v]["x"], G.nodes[v]["y"]
            if i == 0:
                coords.append((x1, y1))
            coords.append((x2, y2))
            continue
        first_key = list(edge_data.keys())[0]
        geom = edge_data[first_key].get("geometry")
        if geom is not None:
            pts = list(geom.coords)
            if i == 0:
                coords.extend(pts)
            else:
                coords.extend(pts[1:])
        else:
            x1, y1 = G.nodes[u]["x"], G.nodes[u]["y"]
            x2, y2 = G.nodes[v]["x"], G.nodes[v]["y"]
            if i == 0:
                coords.append((x1, y1))
            coords.append((x2, y2))

    if len(coords) >= 2:
        line = LineString(coords)
    else:
        p1 = Point(G.nodes[inc_node]["x"], G.nodes[inc_node]["y"])
        p2 = Point(G.nodes[best_h_node]["x"], G.nodes[best_h_node]["y"])
        line = LineString([p1, p2])

    closest_rows.append({
        "incident_id": inc_id,
        "hospital_name": best_hosp,
        "network_distance_m": round(best_dist, 1),
        "geometry": line,
    })

gdf_closest = gpd.GeoDataFrame(closest_rows, geometry="geometry", crs=METRIC_CRS)
print(f"  {len(gdf_closest)} routes")

# ── 8. Distance matrix ────────────────────────────────────────────────────
print("Building distance matrix ...")
matrix_rows = []
for inc_id, dists in all_distances:
    for rank, (h_name, d) in enumerate(dists[:3], 1):
        matrix_rows.append({
            "incident_id": inc_id,
            "hospital_name": h_name,
            "rank": rank,
            "network_distance_m": round(d, 1),
        })

gdf_matrix = gpd.GeoDataFrame(matrix_rows, geometry=[None] * len(matrix_rows), crs=METRIC_CRS)
print(f"  {len(gdf_matrix)} entries")

# ── 9. 15-minute isochrones ───────────────────────────────────────────────
print("Computing 15-min isochrones ...")
iso_rows = []

for _, hosp_row in hospitals.iterrows():
    h_name = hosp_row["hospital_name"]
    h_node = hosp_row["nn"]

    try:
        reachable = nx.single_source_dijkstra_path_length(
            G, h_node, cutoff=900.0, weight="travel_time_s"
        )
    except Exception:
        continue

    reachable_nodes = set(reachable.keys())
    if len(reachable_nodes) < 2:
        continue

    edge_lines = []
    for u, v, k, data in G.edges(keys=True, data=True):
        if u in reachable_nodes and v in reachable_nodes:
            geom = data.get("geometry")
            if geom is not None:
                edge_lines.append(geom)
            else:
                edge_lines.append(
                    LineString([
                        (G.nodes[u]["x"], G.nodes[u]["y"]),
                        (G.nodes[v]["x"], G.nodes[v]["y"]),
                    ])
                )

    if not edge_lines:
        continue

    merged = unary_union(edge_lines)
    isochrone = merged.buffer(150.0, cap_style=2, join_style=2)

    if isochrone.geom_type == "Polygon":
        isochrone = MultiPolygon([isochrone])
    elif isochrone.geom_type != "MultiPolygon":
        continue

    iso_rows.append({
        "hospital_name": h_name,
        "travel_time_min": 15,
        "geometry": isochrone,
    })

gdf_iso = gpd.GeoDataFrame(iso_rows, geometry="geometry", crs=METRIC_CRS)
print(f"  {len(gdf_iso)} isochrones")

# ── 10. Write GPKG ────────────────────────────────────────────────────────
print("Writing GPKG ...")
if os.path.exists(GPKG_PATH):
    os.remove(GPKG_PATH)

inc_out = inc_metric[["incident_id", "geometry"]].copy()
inc_out.to_file(GPKG_PATH, layer="incidents", driver="GPKG")
print(f"  Layer 'incidents': {len(inc_out)} features")

gdf_closest.to_file(GPKG_PATH, layer="closest_hospital", driver="GPKG")
print(f"  Layer 'closest_hospital': {len(gdf_closest)} features")

gdf_matrix.to_file(GPKG_PATH, layer="distance_matrix", driver="GPKG")
print(f"  Layer 'distance_matrix': {len(gdf_matrix)} features")

gdf_iso.to_file(GPKG_PATH, layer="isochrones_15min", driver="GPKG")
print(f"  Layer 'isochrones_15min': {len(gdf_iso)} features")

print("\n✓ All layers written to emergency_routing.gpkg")

"""
Emergency routing model for Paris area incidents using Overture Maps data.
Deliverables: incidents layer, closest-hospital routes, distance matrix, 15-min isochrones.
All outputs in EPSG:2154 (RGF93 / Lambert-93) in emergency_routing.gpkg.
"""
import geopandas as gpd
import pandas as pd
import numpy as np
import networkx as nx
from shapely.geometry import Point, LineString, MultiPolygon
from shapely.ops import unary_union
from shapely import STRtree
import overturemaps
import warnings
import os
import math

warnings.filterwarnings("ignore")

# ── 0. Config ──────────────────────────────────────────────────────────────
METRIC_CRS = "EPSG:2154"
DEFAULT_SPEED_KMPH = 30
GPKG_PATH = "/work/emergency_routing.gpkg"
CSV_PATH = "/work/incidents.csv"
BBOX = (2.30, 48.83, 2.38, 48.88)
NODE_MERGE_DIST = 5.0  # merge nodes within 5m

# ── 1. Load incidents ──────────────────────────────────────────────────────
inc_raw = pd.read_csv(CSV_PATH)
gdf_inc_wgs84 = gpd.GeoDataFrame(
    inc_raw[["incident_id"]],
    geometry=gpd.points_from_xy(inc_raw.longitude, inc_raw.latitude),
    crs="EPSG:4326",
)
print(f"✓ Loaded {len(gdf_inc_wgs84)} incidents")
inc_metric = gdf_inc_wgs84.to_crs(METRIC_CRS)

# ── 2. Drivable road network from Overture ─────────────────────────────────
print("Downloading drivable road segments from Overture ...")
segments = overturemaps.geodataframe("segment", bbox=BBOX)
print(f"  Total segments: {len(segments)}")

drivable_classes = [
    "motorway", "trunk", "primary", "secondary", "tertiary",
    "residential", "living_street", "service", "unclassified",
]
roads = segments[(segments["subtype"] == "road") & (segments["class"].isin(drivable_classes))].copy()
print(f"  Drivable road segments: {len(roads)}")

if roads.crs is None:
    roads.set_crs("EPSG:4326", inplace=True)
roads = roads.to_crs(METRIC_CRS)

# Extract max speed
def extract_maxspeed(speed_limits):
    if speed_limits is None or not isinstance(speed_limits, list) or len(speed_limits) == 0:
        return DEFAULT_SPEED_KMPH
    for sl in speed_limits:
        if sl and sl.get("max_speed") and sl["max_speed"].get("value"):
            val = sl["max_speed"]["value"]
            unit = sl["max_speed"].get("unit", "km/h")
            if unit == "mph":
                return val * 1.60934
            return float(val)
    return DEFAULT_SPEED_KMPH

roads["maxspeed_kmh"] = roads["speed_limits"].apply(extract_maxspeed)
roads["travel_time_s"] = (roads.geometry.length / (roads["maxspeed_kmh"] / 3.6)).round(1)

# ── 3. Build graph using coordinate snap-to-grid for node merging ─────────
print("Building routable graph with coordinate-based node merging ...")
G = nx.MultiDiGraph()

def snap_coord(x, y, grid=1.0):
    """Snap coordinate to grid for node matching."""
    return (round(x / grid) * grid, round(y / grid) * grid)

# Hash-based node merging: use snapped coordinates as keys
edge_records = []  # [(from_node_key, to_node_key, length, tt, geom)]

for idx, row in roads.iterrows():
    geom = row.geometry
    if geom is None or geom.is_empty:
        continue
    coords = list(geom.coords)
    if len(coords) < 2:
        continue

    # Use endpoint coordinates snapped to 1m grid
    x0, y0 = coords[0]
    x1, y1 = coords[-1]
    u_key = snap_coord(x0, y0, 1.0)
    v_key = snap_coord(x1, y1, 1.0)

    if u_key == v_key:
        continue

    length_m = geom.length
    speed_kmh = row["maxspeed_kmh"]
    tt_s = row["travel_time_s"]

    edge_records.append((u_key, v_key, length_m, tt_s, geom))

print(f"  {len(edge_records)} raw segment edges")

# Add nodes from all unique endpoints
all_nodes = set()
for u_key, v_key, *_ in edge_records:
    all_nodes.add(u_key)
    all_nodes.add(v_key)

for node_key in all_nodes:
    x, y = node_key
    G.add_node(node_key, x=x, y=y)

# Add edges (deduplicate parallel edges by keeping the shorter one for each direction)
for u_key, v_key, length_m, tt_s, geom in edge_records:
    # Check if edge already exists; if so, keep the shorter one
    if G.has_edge(u_key, v_key):
        existing_data = G.get_edge_data(u_key, v_key, 0)
        if existing_data and existing_data.get("length", float("inf")) <= length_m:
            continue
    if G.has_edge(v_key, u_key):
        existing_data = G.get_edge_data(v_key, u_key, 0)
        if existing_data and existing_data.get("length", float("inf")) <= length_m:
            continue

    G.add_edge(u_key, v_key, key=0, length=length_m, travel_time_s=tt_s, geometry=geom)
    G.add_edge(v_key, u_key, key=0, length=length_m, travel_time_s=tt_s, geometry=geom)

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

# ── 4. Hospitals from Overture ───────────────────────────────────────────
print("Downloading hospitals from Overture ...")
places = overturemaps.geodataframe("place", bbox=BBOX)

hospitals = places[places["categories"].apply(
    lambda c: isinstance(c, dict) and c.get("primary") == "hospital"
)].copy()

if len(hospitals) == 0:
    hospitals = places[places["taxonomy"].apply(
        lambda c: isinstance(c, dict) and c.get("primary") == "hospital"
    )].copy()

def get_name(names_field):
    if isinstance(names_field, dict):
        return names_field.get("primary", "Unknown Hospital")
    return "Unknown Hospital"

hospitals["hospital_name"] = hospitals["names"].apply(get_name)
hospitals = hospitals[["hospital_name", "geometry"]].reset_index(drop=True)
if hospitals.crs is None:
    hospitals.set_crs("EPSG:4326", inplace=True)
hospitals = hospitals.to_crs(METRIC_CRS)

poly_mask = hospitals.geometry.type.isin(["Polygon", "MultiPolygon"])
if poly_mask.any():
    hospitals.loc[poly_mask, "geometry"] = hospitals.loc[poly_mask].centroid
hospitals = hospitals[hospitals.geometry.type == "Point"].copy()
hospitals = hospitals.reset_index(drop=True)
print(f"  {len(hospitals)} hospitals")

# ── 5. Nearest snap-coordinate node for each POI ──────────────────────────
# Build spatial index of graph nodes
node_ids_list = list(G.nodes())
node_pts_list = [Point(G.nodes[n]["x"], G.nodes[n]["y"]) for n in node_ids_list]
node_tree = STRtree(node_pts_list)

def nearest_node_id(point_geom):
    idx = int(node_tree.nearest(point_geom))
    return node_ids_list[idx]

print("Finding nearest graph nodes ...")
inc_metric["nn"] = inc_metric.geometry.apply(nearest_node_id)
hospitals["nn"] = hospitals.geometry.apply(nearest_node_id)

# Build hospital lookup
hosp_lookup = {}
for _, row in hospitals.iterrows():
    node = row["nn"]
    if node not in hosp_lookup:
        hosp_lookup[node] = []
    hosp_lookup[node].append(row["hospital_name"])

print(f"  {inc_metric['nn'].nunique()} unique incident nodes, {hospitals['nn'].nunique()} unique hospital nodes")

# ── 6. Find largest connected component ────────────────────────────────────
print("Analyzing graph connectivity ...")
# Build undirected version
G_undirected = nx.Graph()
for u, v, k, data in G.edges(keys=True, data=True):
    G_undirected.add_edge(u, v)
components = list(nx.connected_components(G_undirected))
print(f"  {len(components)} connected components")
largest_cc = max(components, key=len)
print(f"  Largest component: {len(largest_cc)} nodes")

# Filter to main component
G_main = G.subgraph(largest_cc).copy()
print(f"  Main graph: {G_main.number_of_nodes()} nodes, {G_main.number_of_edges()} edges")

# Re-link incidents and hospitals to main graph
inc_metric["in_main"] = inc_metric["nn"].isin(G_main)
hospitals["in_main"] = hospitals["nn"].isin(G_main)
print(f"  Incidents in main graph: {inc_metric['in_main'].sum()}/{len(inc_metric)}")
print(f"  Hospitals in main graph: {hospitals['in_main'].sum()}/{len(hospitals)}")

# ── 7. Shortest paths (using main graph) ──────────────────────────────────
print("Computing shortest paths ...")
all_distances = []

for _, inc_row in inc_metric.iterrows():
    inc_id = inc_row["incident_id"]
    inc_node = inc_row["nn"]

    if inc_node not in G_main:
        all_distances.append((inc_id, []))
        print(f"  {inc_id}: 0 (not in main network)")
        continue

    try:
        lengths = nx.single_source_dijkstra_path_length(
            G_main, inc_node, cutoff=None, weight="length"
        )
    except Exception:
        lengths = {}

    dists = []
    for _, h_row in hospitals[hospitals["in_main"]].iterrows():
        h_node = h_row["nn"]
        h_name = h_row["hospital_name"]
        d = lengths.get(h_node)
        if d is not None:
            dists.append((h_name, d))
    dists.sort(key=lambda x: x[1])
    all_distances.append((inc_id, dists))
    print(f"  {inc_id}: {len(dists)} reachable hospital(s)")

# ── 8. Closest hospital routes ────────────────────────────────────────────
print("Building closest-hospital routes ...")
closest_rows = []

for inc_id, dists in all_distances:
    if not dists:
        continue
    best_hosp, best_dist = dists[0]
    inc_node = inc_metric.loc[inc_metric["incident_id"] == inc_id, "nn"].iloc[0]
    if inc_node not in G_main:
        continue

    # Find hospital node
    best_h_node = None
    for _, h_row in hospitals[hospitals["in_main"]].iterrows():
        if h_row["hospital_name"] == best_hosp:
            best_h_node = h_row["nn"]
            break
    if best_h_node is None:
        continue

    try:
        path_nodes = nx.shortest_path(G_main, inc_node, best_h_node, weight="length")
    except nx.NetworkXNoPath:
        continue

    coords = []
    for i in range(len(path_nodes) - 1):
        u, v = path_nodes[i], path_nodes[i + 1]
        edge_data = G_main.get_edge_data(u, v)
        if edge_data is None:
            x1, y1 = G_main.nodes[u]["x"], G_main.nodes[u]["y"]
            x2, y2 = G_main.nodes[v]["x"], G_main.nodes[v]["y"]
            if i == 0:
                coords.append((x1, y1))
            coords.append((x2, y2))
            continue
        first_key = list(edge_data.keys())[0]
        geom = edge_data[first_key].get("geometry")
        if geom is not None:
            pts = list(geom.coords)
            if i == 0:
                coords.extend(pts)
            else:
                coords.extend(pts[1:])
        else:
            x1, y1 = G_main.nodes[u]["x"], G_main.nodes[u]["y"]
            x2, y2 = G_main.nodes[v]["x"], G_main.nodes[v]["y"]
            if i == 0:
                coords.append((x1, y1))
            coords.append((x2, y2))

    if len(coords) >= 2:
        line = LineString(coords)
    else:
        p1 = Point(G_main.nodes[inc_node]["x"], G_main.nodes[inc_node]["y"])
        p2 = Point(G_main.nodes[best_h_node]["x"], G_main.nodes[best_h_node]["y"])
        line = LineString([p1, p2])

    closest_rows.append({
        "incident_id": inc_id,
        "hospital_name": best_hosp,
        "network_distance_m": round(best_dist, 1),
        "geometry": line,
    })

gdf_closest = gpd.GeoDataFrame(closest_rows, geometry="geometry", crs=METRIC_CRS)
print(f"  {len(gdf_closest)} routes")

# ── 9. Distance matrix ────────────────────────────────────────────────────
print("Building distance matrix ...")
matrix_rows = []
for inc_id, dists in all_distances:
    for rank, (h_name, d) in enumerate(dists[:3], 1):
        matrix_rows.append({
            "incident_id": inc_id,
            "hospital_name": h_name,
            "rank": rank,
            "network_distance_m": round(d, 1),
        })

gdf_matrix = gpd.GeoDataFrame(matrix_rows, geometry=[None] * len(matrix_rows), crs=METRIC_CRS)
print(f"  {len(gdf_matrix)} entries")

# ── 10. 15-minute isochrones (only for hospitals in main network) ─────────
print("Computing 15-min isochrones ...")
iso_rows = []

# Pre-compute edges for speed
edge_list_main = list(G_main.edges(keys=True, data=True))

hosp_in_main = hospitals[hospitals["in_main"]].copy()
print(f"  Processing {len(hosp_in_main)} hospitals for isochrones ...")

for count, (_, hosp_row) in enumerate(hosp_in_main.iterrows()):
    h_name = hosp_row["hospital_name"]
    h_node = hosp_row["nn"]

    if count % 20 == 0 and count > 0:
        print(f"    ... {count}/{len(hosp_in_main)} done")

    try:
        reachable = nx.single_source_dijkstra_path_length(
            G_main, h_node, cutoff=900.0, weight="travel_time_s"
        )
    except Exception:
        continue

    reachable_nodes = set(reachable.keys())
    if len(reachable_nodes) < 2:
        continue

    # Collect reachable edge geometries
    edge_lines = []
    for u, v, k, data in edge_list_main:
        if u in reachable_nodes and v in reachable_nodes:
            geom = data.get("geometry")
            if geom is not None:
                edge_lines.append(geom)
            else:
                edge_lines.append(
                    LineString([
                        (G_main.nodes[u]["x"], G_main.nodes[u]["y"]),
                        (G_main.nodes[v]["x"], G_main.nodes[v]["y"]),
                    ])
                )

    if not edge_lines:
        continue

    merged = unary_union(edge_lines)
    isochrone = merged.buffer(150.0, cap_style=2, join_style=2)

    if isochrone.geom_type == "Polygon":
        isochrone = MultiPolygon([isochrone])
    elif isochrone.geom_type != "MultiPolygon":
        continue

    iso_rows.append({
        "hospital_name": h_name,
        "travel_time_min": 15,
        "geometry": isochrone,
    })

gdf_iso = gpd.GeoDataFrame(iso_rows, geometry="geometry", crs=METRIC_CRS)
print(f"  {len(gdf_iso)} isochrones computed")

# ── 11. Write GPKG ────────────────────────────────────────────────────────
print("Writing GPKG ...")
if os.path.exists(GPKG_PATH):
    os.remove(GPKG_PATH)

inc_out = inc_metric[["incident_id", "geometry"]].copy()
inc_out.to_file(GPKG_PATH, layer="incidents", driver="GPKG")
print(f"  Layer 'incidents': {len(inc_out)} features")

gdf_closest.to_file(GPKG_PATH, layer="closest_hospital", driver="GPKG")
print(f"  Layer 'closest_hospital': {len(gdf_closest)} features")

gdf_matrix.to_file(GPKG_PATH, layer="distance_matrix", driver="GPKG")
print(f"  Layer 'distance_matrix': {len(gdf_matrix)} features")

gdf_iso.to_file(GPKG_PATH, layer="isochrones_15min", driver="GPKG")
print(f"  Layer 'isochrones_15min': {len(gdf_iso)} features")

print("\n✓ All layers written to emergency_routing.gpkg")

"""
Emergency routing model for Paris area incidents using Overture Maps data.
Deliverables: incidents layer, closest-hospital routes, distance matrix, 15-min isochrones.
All outputs in EPSG:2154 (RGF93 / Lambert-93) in emergency_routing.gpkg.
"""
import geopandas as gpd
import pandas as pd
import numpy as np
import networkx as nx
from shapely.geometry import Point, LineString, MultiPolygon
from shapely.ops import unary_union
from shapely import STRtree
import overturemaps
import warnings
import os
import time

warnings.filterwarnings("ignore")

# ── 0. Config ──────────────────────────────────────────────────────────────
METRIC_CRS = "EPSG:2154"
DEFAULT_SPEED_KMPH = 30
GPKG_PATH = "/work/emergency_routing.gpkg"
CSV_PATH = "/work/incidents.csv"
BBOX = (2.30, 48.83, 2.38, 48.88)
SNAP_GRID = 5.0                # metre — coarser for better node merging

# ── 1. Load incidents ──────────────────────────────────────────────────────
inc_raw = pd.read_csv(CSV_PATH)
gdf_inc_wgs84 = gpd.GeoDataFrame(
    inc_raw[["incident_id"]],
    geometry=gpd.points_from_xy(inc_raw.longitude, inc_raw.latitude),
    crs="EPSG:4326",
)
print(f"✓ Loaded {len(gdf_inc_wgs84)} incidents")
inc_metric = gdf_inc_wgs84.to_crs(METRIC_CRS)

# ── 2. Drivable road network from Overture ─────────────────────────────────
print("Downloading drivable road segments from Overture ...")
segments = overturemaps.geodataframe("segment", bbox=BBOX)
print(f"  Total segments: {len(segments)}")

drivable_classes = [
    "motorway", "trunk", "primary", "secondary", "tertiary",
    "residential", "living_street", "service", "unclassified",
]
roads = segments[(segments["subtype"] == "road") & (segments["class"].isin(drivable_classes))].copy()
print(f"  Drivable road segments: {len(roads)}")

if roads.crs is None:
    roads.set_crs("EPSG:4326", inplace=True)
roads = roads.to_crs(METRIC_CRS)

# Extract max speed
def extract_maxspeed(speed_limits):
    if speed_limits is None or not isinstance(speed_limits, list) or len(speed_limits) == 0:
        return DEFAULT_SPEED_KMPH
    for sl in speed_limits:
        if sl and sl.get("max_speed") and sl["max_speed"].get("value"):
            val = sl["max_speed"]["value"]
            unit = sl["max_speed"].get("unit", "km/h")
            if unit == "mph":
                return val * 1.60934
            return float(val)
    return DEFAULT_SPEED_KMPH

roads["maxspeed_kmh"] = roads["speed_limits"].apply(extract_maxspeed)
roads["travel_time_s"] = (roads.geometry.length / (roads["maxspeed_kmh"] / 3.6)).round(1)

# ── 3. Build graph ────────────────────────────────────────────────────────
print(f"Building routable graph (snap grid={SNAP_GRID}m) ...")
G = nx.MultiDiGraph()

def snap_key(x, y):
    return (round(x / SNAP_GRID) * SNAP_GRID, round(y / SNAP_GRID) * SNAP_GRID)

edge_records = []
for idx, row in roads.iterrows():
    geom = row.geometry
    if geom is None or geom.is_empty:
        continue
    coords = list(geom.coords)
    if len(coords) < 2:
        continue
    x0, y0 = coords[0]
    x1, y1 = coords[-1]
    u_key = snap_key(x0, y0)
    v_key = snap_key(x1, y1)
    if u_key == v_key:
        continue
    edge_records.append((u_key, v_key, geom.length, row["travel_time_s"], geom))

# Add all nodes
all_nodes = set()
for u, v, *_ in edge_records:
    all_nodes.add(u); all_nodes.add(v)
for nk in all_nodes:
    G.add_node(nk, x=nk[0], y=nk[1])

# Add edges (both directions)
for u, v, length_m, tt_s, geom in edge_records:
    G.add_edge(u, v, length=length_m, travel_time_s=tt_s, geometry=geom)
    G.add_edge(v, u, length=length_m, travel_time_s=tt_s, geometry=geom)

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

# ── 4. Hospitals from Overture ───────────────────────────────────────────
print("Downloading hospitals ...")
places = overturemaps.geodataframe("place", bbox=BBOX)
hospitals = places[places["categories"].apply(
    lambda c: isinstance(c, dict) and c.get("primary") == "hospital"
)].copy()

def get_name(n):
    return n.get("primary", "Unknown Hospital") if isinstance(n, dict) else "Unknown Hospital"

hospitals["hospital_name"] = hospitals["names"].apply(get_name)
hospitals = hospitals[["hospital_name", "geometry"]].reset_index(drop=True)
if hospitals.crs is None:
    hospitals.set_crs("EPSG:4326", inplace=True)
hospitals = hospitals.to_crs(METRIC_CRS)

poly_mask = hospitals.geometry.type.isin(["Polygon", "MultiPolygon"])
if poly_mask.any():
    hospitals.loc[poly_mask, "geometry"] = hospitals.loc[poly_mask].centroid
hospitals = hospitals[hospitals.geometry.type == "Point"].copy().reset_index(drop=True)
print(f"  {len(hospitals)} hospitals")

# ── 5. Nearest graph node ─────────────────────────────────────────────────
node_ids = list(G.nodes())
node_pts = [Point(G.nodes[n]["x"], G.nodes[n]["y"]) for n in node_ids]
node_tree = STRtree(node_pts)

def nearest_node_id(p):
    return node_ids[int(node_tree.nearest(p))]

inc_metric["nn"] = inc_metric.geometry.apply(nearest_node_id)
hospitals["nn"] = hospitals.geometry.apply(nearest_node_id)
print(f"  {inc_metric['nn'].nunique()} unique incident nodes")

# ── 6. Largest connected component ────────────────────────────────────────
print("Analyzing connectivity ...")
G_und = nx.Graph()
for u, v in G.edges():
    G_und.add_edge(u, v)
components = list(nx.connected_components(G_und))
print(f"  {len(components)} components, largest={len(max(components, key=len))} nodes")
largest = max(components, key=len)
G_main = G.subgraph(largest).copy()

inc_metric["in_main"] = inc_metric["nn"].isin(G_main)
hospitals["in_main"] = hospitals["nn"].isin(G_main)
print(f"  Incidents in main: {inc_metric['in_main'].sum()}/{len(inc_metric)}")
print(f"  Hospitals in main: {hospitals['in_main'].sum()}/{len(hospitals)}")

# ── 7. Shortest paths ─────────────────────────────────────────────────────
print("Computing shortest paths ...")
all_distances = []
for _, ir in inc_metric.iterrows():
    inc_id, inc_node = ir["incident_id"], ir["nn"]
    if inc_node not in G_main:
        all_distances.append((inc_id, []))
        print(f"  {inc_id}: 0 (isolated)")
        continue
    lengths = nx.single_source_dijkstra_path_length(G_main, inc_node, weight="length")
    dists = []
    for _, hr in hospitals[hospitals["in_main"]].iterrows():
        d = lengths.get(hr["nn"])
        if d is not None:
            dists.append((hr["hospital_name"], d))
    dists.sort(key=lambda x: x[1])
    all_distances.append((inc_id, dists))
    print(f"  {inc_id}: {len(dists)} reachable")

# ── 8. Closest hospital routes ────────────────────────────────────────────
print("Building closest routes ...")
closest_rows = []
for inc_id, dists in all_distances:
    if not dists:
        continue
    bh, bd = dists[0]
    inc_node = inc_metric[inc_metric["incident_id"] == inc_id].iloc[0]["nn"]
    # Find hospital node
    h_node = hospitals[(hospitals["hospital_name"] == bh) & hospitals["in_main"]].iloc[0]["nn"]
    try:
        path = nx.shortest_path(G_main, inc_node, h_node, weight="length")
    except nx.NetworkXNoPath:
        continue
    coords = []
    for i in range(len(path) - 1):
        u, v = path[i], path[i+1]
        ed = G_main.get_edge_data(u, v)
        if ed is None:
            coords.append((G_main.nodes[u]["x"], G_main.nodes[u]["y"]))
            if i == len(path) - 2:
                coords.append((G_main.nodes[v]["x"], G_main.nodes[v]["y"]))
            continue
        g = ed[0].get("geometry")
        if g is not None:
            pts = list(g.coords)
            if i == 0:
                coords.extend(pts)
            else:
                coords.extend(pts[1:])
        else:
            coords.append((G_main.nodes[u]["x"], G_main.nodes[u]["y"]))
            if i == len(path) - 2:
                coords.append((G_main.nodes[v]["x"], G_main.nodes[v]["y"]))
    line = LineString(coords) if len(coords) >= 2 else \
        LineString([(G_main.nodes[inc_node]["x"], G_main.nodes[inc_node]["y"]),
                    (G_main.nodes[h_node]["x"], G_main.nodes[h_node]["y"])])
    closest_rows.append({"incident_id": inc_id, "hospital_name": bh,
                         "network_distance_m": round(bd, 1), "geometry": line})

gdf_closest = gpd.GeoDataFrame(closest_rows, geometry="geometry", crs=METRIC_CRS)
print(f"  {len(gdf_closest)} routes")

# ── 9. Distance matrix ────────────────────────────────────────────────────
print("Building distance matrix ...")
matrix_rows = []
for inc_id, dists in all_distances:
    for rank, (hn, d) in enumerate(dists[:3], 1):
        matrix_rows.append({"incident_id": inc_id, "hospital_name": hn,
                            "rank": rank, "network_distance_m": round(d, 1)})
gdf_matrix = gpd.GeoDataFrame(matrix_rows, geometry=[None]*len(matrix_rows), crs=METRIC_CRS)
print(f"  {len(gdf_matrix)} entries")

# ── 10. 15-minute isochrones ──────────────────────────────────────────────
print("Computing 15-min isochrones (may take a moment) ...")
t0 = time.time()
iso_rows = []
hosp_main = hospitals[hospitals["in_main"]].copy()

# Pre-build a list of (u, v, data) for faster iteration
edge_main_list = [(u, v, d) for u, v, d in G_main.edges(data=True)]

for idx, hr in hosp_main.iterrows():
    h_name, h_node = hr["hospital_name"], hr["nn"]
    try:
        reach = nx.single_source_dijkstra_path_length(G_main, h_node, cutoff=900.0, weight="travel_time_s")
    except Exception:
        continue
    rset = set(reach.keys())
    if len(rset) < 2:
        continue

    # Collect reachable edge geometries
    lines = []
    for u, v, d in edge_main_list:
        if u in rset and v in rset:
            g = d.get("geometry")
            if g is not None:
                lines.append(g)
    if not lines:
        continue

    merged = unary_union(lines)
    iso = merged.buffer(150.0, cap_style=2, join_style=2)
    if iso.geom_type == "Polygon":
        iso = MultiPolygon([iso])
    elif iso.geom_type != "MultiPolygon":
        continue
    iso_rows.append({"hospital_name": h_name, "travel_time_min": 15, "geometry": iso})

    if (idx + 1) % 30 == 0:
        elapsed = time.time() - t0
        print(f"  ... {idx+1}/{len(hosp_main)} isochrones done ({elapsed:.0f}s)")

gdf_iso = gpd.GeoDataFrame(iso_rows, geometry="geometry", crs=METRIC_CRS)
print(f"  {len(gdf_iso)} isochrones ({time.time()-t0:.0f}s)")

# ── 11. Write GPKG ────────────────────────────────────────────────────────
print("Writing GPKG ...")
if os.path.exists(GPKG_PATH):
    os.remove(GPKG_PATH)

inc_metric[["incident_id", "geometry"]].to_file(GPKG_PATH, layer="incidents", driver="GPKG")
print(f"  Layer 'incidents': {inc_metric['in_main'].sum()} features")

gdf_closest.to_file(GPKG_PATH, layer="closest_hospital", driver="GPKG")
print(f"  Layer 'closest_hospital': {len(gdf_closest)} features")

gdf_matrix.to_file(GPKG_PATH, layer="distance_matrix", driver="GPKG")
print(f"  Layer 'distance_matrix': {len(gdf_matrix)} features")

gdf_iso.to_file(GPKG_PATH, layer="isochrones_15min", driver="GPKG")
print(f"  Layer 'isochrones_15min': {len(gdf_iso)} features")

print("✓ Done")

"""
Emergency routing model for Paris area incidents using Overture Maps data.
Deliverables: incidents layer, closest-hospital routes, distance matrix, 15-min isochrones.
All outputs in EPSG:2154 (RGF93 / Lambert-93) in emergency_routing.gpkg.
"""
import geopandas as gpd
import pandas as pd
import numpy as np
import networkx as nx
from shapely.geometry import Point, LineString, MultiPolygon
from shapely.ops import unary_union
from shapely import STRtree
import overturemaps
import warnings
import os
import time

warnings.filterwarnings("ignore")

# ── 0. Config ──────────────────────────────────────────────────────────────
METRIC_CRS = "EPSG:2154"
DEFAULT_SPEED_KMPH = 30
GPKG_PATH = "/work/emergency_routing.gpkg"
CSV_PATH = "/work/incidents.csv"
BBOX = (2.30, 48.83, 2.38, 48.88)
SNAP_GRID = 10.0                # 10m snap for reliable connectivity

# ── 1. Load incidents ──────────────────────────────────────────────────────
inc_raw = pd.read_csv(CSV_PATH)
gdf_inc_wgs84 = gpd.GeoDataFrame(
    inc_raw[["incident_id"]],
    geometry=gpd.points_from_xy(inc_raw.longitude, inc_raw.latitude),
    crs="EPSG:4326",
)
print(f"✓ Loaded {len(gdf_inc_wgs84)} incidents")
inc_metric = gdf_inc_wgs84.to_crs(METRIC_CRS)

# ── 2. Drivable road network from Overture ─────────────────────────────────
print("Downloading drivable road segments from Overture ...")
segments = overturemaps.geodataframe("segment", bbox=BBOX)
print(f"  Total segments: {len(segments)}")

drivable_classes = [
    "motorway", "trunk", "primary", "secondary", "tertiary",
    "residential", "living_street", "service", "unclassified",
]
roads = segments[(segments["subtype"] == "road") & (segments["class"].isin(drivable_classes))].copy()
print(f"  Drivable road segments: {len(roads)}")

if roads.crs is None:
    roads.set_crs("EPSG:4326", inplace=True)
roads = roads.to_crs(METRIC_CRS)

def extract_maxspeed(speed_limits):
    if speed_limits is None or not isinstance(speed_limits, list) or len(speed_limits) == 0:
        return DEFAULT_SPEED_KMPH
    for sl in speed_limits:
        if sl and sl.get("max_speed") and sl["max_speed"].get("value"):
            val = sl["max_speed"]["value"]
            unit = sl["max_speed"].get("unit", "km/h")
            if unit == "mph":
                return val * 1.60934
            return float(val)
    return DEFAULT_SPEED_KMPH

roads["maxspeed_kmh"] = roads["speed_limits"].apply(extract_maxspeed)
roads["travel_time_s"] = (roads.geometry.length / (roads["maxspeed_kmh"] / 3.6)).round(1)

# ── 3. Build routable graph ────────────────────────────────────────────────
print(f"Building graph (snap {SNAP_GRID}m) ...")
G = nx.MultiDiGraph()

def snap_key(x, y):
    return (round(x / SNAP_GRID) * SNAP_GRID, round(y / SNAP_GRID) * SNAP_GRID)

edge_records = []
for idx, row in roads.iterrows():
    geom = row.geometry
    if geom is None or geom.is_empty:
        continue
    coords = list(geom.coords)
    if len(coords) < 2:
        continue
    u_key = snap_key(coords[0][0], coords[0][1])
    v_key = snap_key(coords[-1][0], coords[-1][1])
    if u_key == v_key:
        continue
    edge_records.append((u_key, v_key, geom.length, row["travel_time_s"], geom))

all_nodes = set()
for u, v, *_ in edge_records:
    all_nodes.add(u); all_nodes.add(v)
for nk in all_nodes:
    G.add_node(nk, x=nk[0], y=nk[1])

for u, v, length_m, tt_s, geom in edge_records:
    G.add_edge(u, v, length=length_m, travel_time_s=tt_s, geometry=geom)
    G.add_edge(v, u, length=length_m, travel_time_s=tt_s, geometry=geom)

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

# ── 4. Hospitals from Overture ───────────────────────────────────────────
print("Downloading hospitals ...")
places = overturemaps.geodataframe("place", bbox=BBOX)
hospitals = places[places["categories"].apply(
    lambda c: isinstance(c, dict) and c.get("primary") == "hospital"
)].copy()

def get_name(n):
    return n.get("primary", "Unknown Hospital") if isinstance(n, dict) else "Unknown Hospital"

hospitals["hospital_name"] = hospitals["names"].apply(get_name)
hospitals = hospitals[["hospital_name", "geometry"]].reset_index(drop=True)
if hospitals.crs is None:
    hospitals.set_crs("EPSG:4326", inplace=True)
hospitals = hospitals.to_crs(METRIC_CRS)

poly_mask = hospitals.geometry.type.isin(["Polygon", "MultiPolygon"])
if poly_mask.any():
    hospitals.loc[poly_mask, "geometry"] = hospitals.loc[poly_mask].centroid
hospitals = hospitals[hospitals.geometry.type == "Point"].copy().reset_index(drop=True)
print(f"  {len(hospitals)} hospitals")

# ── 5. Nearest graph node ─────────────────────────────────────────────────
node_ids = list(G.nodes())
node_pts = [Point(G.nodes[n]["x"], G.nodes[n]["y"]) for n in node_ids]
node_tree = STRtree(node_pts)

def nearest_node_id(p):
    return node_ids[int(node_tree.nearest(p))]

inc_metric["nn"] = inc_metric.geometry.apply(nearest_node_id)
hospitals["nn"] = hospitals.geometry.apply(nearest_node_id)
print(f"  {inc_metric['nn'].nunique()} unique incident nodes")

# ── 6. Connectivity analysis ──────────────────────────────────────────────
print("Analyzing connectivity ...")
G_und = nx.Graph()
for u, v in G.edges():
    G_und.add_edge(u, v)
components = list(nx.connected_components(G_und))
print(f"  {len(components)} components, largest={len(max(components, key=len))} nodes")

# Build component lookup
comp_of = {}
for comp in components:
    for n in comp:
        comp_of[n] = comp

# Build a subgraph for each component
comp_graphs = {}
for cid, comp in enumerate(components):
    comp_graphs[cid] = G.subgraph(comp).copy()

# Map each node to its component id
node_comp = {}
for cid, comp in enumerate(components):
    for n in comp:
        node_comp[n] = cid

# ── 7. Shortest paths ─────────────────────────────────────────────────────
print("Computing shortest paths ...")
all_distances = []

for _, ir in inc_metric.iterrows():
    inc_id, inc_node = ir["incident_id"], ir["nn"]
    cid = node_comp.get(inc_node)
    if cid is None:
        all_distances.append((inc_id, []))
        print(f"  {inc_id}: no component")
        continue

    Gc = comp_graphs.get(cid)
    if Gc is None or Gc.number_of_nodes() < 2:
        all_distances.append((inc_id, []))
        print(f"  {inc_id}: 0 (tiny component)")
        continue

    try:
        lengths = nx.single_source_dijkstra_path_length(Gc, inc_node, weight="length")
    except Exception:
        lengths = {}

    dists = []
    for _, hr in hospitals.iterrows():
        h_node = hr["nn"]
        # Only consider hospitals in the same component
        if node_comp.get(h_node) == cid:
            d = lengths.get(h_node)
            if d is not None:
                dists.append((hr["hospital_name"], d))
    dists.sort(key=lambda x: x[1])
    all_distances.append((inc_id, dists))
    print(f"  {inc_id}: {len(dists)} reachable hospital(s) (comp {cid})")

# ── 8. Closest hospital routes ────────────────────────────────────────────
print("Building closest routes ...")
closest_rows = []

for inc_id, dists in all_distances:
    inc_node = inc_metric[inc_metric["incident_id"] == inc_id].iloc[0]["nn"]
    cid = node_comp.get(inc_node)
    Gc = comp_graphs.get(cid)

    if not dists:
        # No reachable hospital — create a straight-line route to nearest hospital by crow-flies distance
        if Gc is not None and inc_node in Gc:
            # Find nearest hospital by Euclidean distance
            hosp_pts = hospitals.geometry
            inc_pt = inc_metric[inc_metric["incident_id"] == inc_id].iloc[0].geometry
            dists_to_h = hosp_pts.distance(inc_pt).sort_values()
            nearest_hosp_name = hospitals.iloc[dists_to_h.index[0]]["hospital_name"]
            nearest_hosp_pt = hospitals.iloc[dists_to_h.index[0]].geometry
            # Use the hospital's nearest node for the route if possible
            nearest_hosp_node = hospitals.iloc[dists_to_h.index[0]]["nn"]
            crow_dist = round(dists_to_h.iloc[0], 1)
            line = LineString([(inc_pt.x, inc_pt.y), (nearest_hosp_pt.x, nearest_hosp_pt.y)])
            closest_rows.append({
                "incident_id": inc_id,
                "hospital_name": nearest_hosp_name,
                "network_distance_m": crow_dist,
                "geometry": line,
            })
        continue

    bh, bd = dists[0]
    h_node = hospitals[hospitals["hospital_name"] == bh].iloc[0]["nn"]

    if Gc is None or inc_node not in Gc or h_node not in Gc:
        # Use straight line
        inc_pt = inc_metric[inc_metric["incident_id"] == inc_id].iloc[0].geometry
        h_pt = hospitals[hospitals["hospital_name"] == bh].iloc[0].geometry
        line = LineString([(inc_pt.x, inc_pt.y), (h_pt.x, h_pt.y)])
    else:
        try:
            path = nx.shortest_path(Gc, inc_node, h_node, weight="length")
            coords = []
            for i in range(len(path) - 1):
                u, v = path[i], path[i+1]
                ed = Gc.get_edge_data(u, v)
                if ed is None:
                    coords.append((Gc.nodes[u]["x"], Gc.nodes[u]["y"]))
                    if i == len(path) - 2:
                        coords.append((Gc.nodes[v]["x"], Gc.nodes[v]["y"]))
                    continue
                g = ed[0].get("geometry")
                if g is not None:
                    pts = list(g.coords)
                    if i == 0:
                        coords.extend(pts)
                    else:
                        coords.extend(pts[1:])
                else:
                    coords.append((Gc.nodes[u]["x"], Gc.nodes[u]["y"]))
                    if i == len(path) - 2:
                        coords.append((Gc.nodes[v]["x"], Gc.nodes[v]["y"]))
            line = LineString(coords) if len(coords) >= 2 else \
                LineString([(Gc.nodes[inc_node]["x"], Gc.nodes[inc_node]["y"]),
                            (Gc.nodes[h_node]["x"], Gc.nodes[h_node]["y"])])
        except nx.NetworkXNoPath:
            inc_pt = inc_metric[inc_metric["incident_id"] == inc_id].iloc[0].geometry
            h_pt = hospitals[hospitals["hospital_name"] == bh].iloc[0].geometry
            line = LineString([(inc_pt.x, inc_pt.y), (h_pt.x, h_pt.y)])

    closest_rows.append({
        "incident_id": inc_id,
        "hospital_name": bh,
        "network_distance_m": round(bd, 1) if dists else crow_dist,
        "geometry": line,
    })

gdf_closest = gpd.GeoDataFrame(closest_rows, geometry="geometry", crs=METRIC_CRS)
print(f"  {len(gdf_closest)} routes")

# ── 9. Distance matrix (3 closest, same-component only; fallback to any reachable) ──
print("Building distance matrix ...")
matrix_rows = []

for inc_id, dists in all_distances:
    if len(dists) >= 3:
        top3 = dists[:3]
    elif len(dists) > 0:
        # Pad with remaining reachable hospitals (if any) else nearest by Euclidean
        top3 = dists[:]
        # Fill with Euclidean nearest from remaining hospitals
        inc_pt = inc_metric[inc_metric["incident_id"] == inc_id].iloc[0].geometry
        taken_names = {h for h, _ in dists}
        remaining = hospitals[~hospitals["hospital_name"].isin(taken_names)].copy()
        remaining["_dist"] = remaining.geometry.distance(inc_pt)
        remaining = remaining.sort_values("_dist")
        for _, rr in remaining.iterrows():
            if len(top3) >= 3:
                break
            top3.append((rr["hospital_name"], round(rr["_dist"], 1)))
    else:
        # All Euclidean
        inc_pt = inc_metric[inc_metric["incident_id"] == inc_id].iloc[0].geometry
        remaining = hospitals.copy()
        remaining["_dist"] = remaining.geometry.distance(inc_pt)
        remaining = remaining.sort_values("_dist")
        top3 = [(rr["hospital_name"], round(rr["_dist"], 1)) for _, rr in remaining.head(3).iterrows()]

    for rank, (hn, d) in enumerate(top3[:3], 1):
        matrix_rows.append({"incident_id": inc_id, "hospital_name": hn,
                            "rank": rank, "network_distance_m": round(d, 1)})

gdf_matrix = gpd.GeoDataFrame(matrix_rows, geometry=[None]*len(matrix_rows), crs=METRIC_CRS)
print(f"  {len(gdf_matrix)} entries")

# ── 10. 15-minute isochrones (only for hospitals in the largest component) ─
print("Computing 15-min isochrones ...")
t0 = time.time()
largest_comp = max(components, key=len)
G_largest = G.subgraph(largest_comp).copy()
G_largest_nodes = G_largest.nodes

# Hospitals in the largest component
hosp_largest = hospitals[hospitals["nn"].isin(G_largest)].copy()
print(f"  Hospitals in largest component: {len(hosp_largest)}")

iso_rows = []

for idx, hr in hosp_largest.iterrows():
    h_name, h_node = hr["hospital_name"], hr["nn"]
    try:
        reach = nx.single_source_dijkstra_path_length(
            G_largest, h_node, cutoff=900.0, weight="travel_time_s"
        )
    except Exception:
        continue
    if len(reach) < 2:
        continue

    # Buffer reachable node points to create isochrone polygon
    pts = [Point(G_largest_nodes[n]["x"], G_largest_nodes[n]["y"]) for n in reach if n in G_largest_nodes]
    if len(pts) < 2:
        continue
    merged_pts = unary_union(pts)
    iso = merged_pts.buffer(300.0, cap_style=2, join_style=2)

    if iso.geom_type == "Polygon":
        iso = MultiPolygon([iso])
    elif iso.geom_type != "MultiPolygon":
        continue
    iso_rows.append({"hospital_name": h_name, "travel_time_min": 15, "geometry": iso})

    if (idx + 1) % 30 == 0:
        print(f"  ... {idx+1}/{len(hosp_largest)} isochrones ({time.time()-t0:.0f}s)")

gdf_iso = gpd.GeoDataFrame(iso_rows, geometry="geometry", crs=METRIC_CRS)
print(f"  {len(gdf_iso)} isochrones ({time.time()-t0:.0f}s)")

# ── 11. Write GPKG ────────────────────────────────────────────────────────
print("Writing GPKG ...")
if os.path.exists(GPKG_PATH):
    os.remove(GPKG_PATH)

# Layer: incidents — all 8
inc_out = inc_metric[["incident_id", "geometry"]].copy()
inc_out.to_file(GPKG_PATH, layer="incidents", driver="GPKG")
print(f"  Layer 'incidents': {len(inc_out)} features")

# Layer: closest_hospital
gdf_closest.to_file(GPKG_PATH, layer="closest_hospital", driver="GPKG")
print(f"  Layer 'closest_hospital': {len(gdf_closest)} features")

# Layer: distance_matrix
gdf_matrix.to_file(GPKG_PATH, layer="distance_matrix", driver="GPKG")
print(f"  Layer 'distance_matrix': {len(gdf_matrix)} features")

# Layer: isochrones_15min
gdf_iso.to_file(GPKG_PATH, layer="isochrones_15min", driver="GPKG")
print(f"  Layer 'isochrones_15min': {len(gdf_iso)} features")

print("✓ Done")