spa-l3-paris-emergency-routing

Rerun

import osmnx as ox
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point, LineString, MultiLineString, Polygon, MultiPolygon
import networkx as nx
import numpy as np

# 1. Load Incidents
incidents_df = pd.read_csv('incidents.csv')
incidents_gdf = gpd.GeoDataFrame(
    incidents_df, 
    geometry=gpd.points_from_xy(incidents_df.longitude, incidents_df.latitude),
    crs="EPSG:4326"
)

# Target CRS: Paris metric (EPSG:2154 - Lambert-93)
target_crs = "EPSG:2154"
incidents_gdf = incidents_gdf.to_crs(target_crs)

# Bounding box for OSM (from prompt)
bbox = (48.83, 2.30, 48.88, 2.38) # minlat, minlon, maxlat, maxlon

# 2. Fetch OSM Data
print("Fetching road network and hospitals...")
# Fetch graph
G = ox.graph_from_bbox(bbox[0], bbox[2], bbox[1], bbox[3], network_type='drive')

# Fetch hospitals
# Overpass query for hospitals (amenity=hospital)
# Note: ox.geometries (deprecated) or ox.features_from_bbox (new)
hospitals = ox.features_from_bbox(bbox[0], bbox[2], bbox[1], bbox[3], tags={'amenity': 'hospital'})

# Pre-process graph for speed and speed assumptions
G = ox.add_edge_speeds(G)
G = ox.add_edge_travel_times(G)

# Ensure all edges have a travel time based on 30km/h if missing
# (ox.add_edge_speeds handles maxspeed, but if still nan, we set to 30km/h)
for u, v, k, data in G.edges(data=True, keys=True):
    if 'travel_time' not in data or np.isnan(data['travel_time']):
        # 30 km/h = 30000 m / 3600 s = 8.333 m/s
        # length in meters
        length = data.get('length', 1.0) 
        data['travel_time'] = length / 8.333

# Prepare hospital list
# We need the nearest node for each hospital to use in network calculations
# and for distances.
hospitals = hospitals.to_crs(target_crs)
# We'll use the centroid of the hospital if it's a polygon, or the point itself
hospitals['geometry'] = hospitals.geometry.centroid

# 3. Analysis Functions
def get_nearest_node(G, point):
    return ox.nearest_nodes(G, point.x, point.y)

# For distances, we will find the closest hospital to each incident and the top 3.
# First, let's find the nearest nodes for each hospital
hospitals['nearest_node'] = hospitals.geometry.apply(lambda p: get_nearest_node(G, p))

# Distance Matrix & Closest Hospital
results_closest = []
results_matrix = []

for _, incident in incidents_gdf.iterrows():
    inc_id = incident['incident_id']
    inc_node = get_nearest_node(G, incident.geometry)
    
    # Calculate distances to all hospitals via network
    # To be efficient, we find shortest path distances from inc_node to all hospital nodes
    # However, we only have a few hospitals, so it's fine to loop.
    
    h_distances = []
    
    for _, hosp in hospitals.iterrows():
        h_node = hosp['nearest_node']
        h_name = hosp.get('name', 'Unnamed Hospital')
        
        try:
            # Shortest path distance (length)
            dist = nx.shortest_path_length(G, inc_node, h_node, weight='length')
            # Shortest path time
            time = nx.shortest_path_length(G, inc_node, h_node, weight='travel_time')
            
            # Path geometry for the closest one
            path_nodes = nx.shortest_path(G, inc_node, h_node, weight='length')
            path_geom = ox.utils_graph.get_route_edge_attributes(G, path_nodes, 'geometry')
            # This is complex, simpler way to get geometry:
            # nx.shortest_path returns nodes, we need to collect edges.
            
            h_distances.append({
                'hospital_name': h_name,
                'network_distance_m': dist,
                'travel_time_s': time,
                'path_nodes': path_nodes
            })
        except nx.NetworkXNoPath:
            continue
            
    # Sort distances for this incident
    h_distances.sort(key=lambda x: x['network_distance_m'])
    
    # Closest hospital (rank 1)
    if h_distances:
        best = h_distances[0]
        # Create LineString for the route
        # We'll use ox.routing.route_to_gdf or manual reconstruction
        # Manually:
        edges = []
        for i in range(len(best['path_nodes']) - 1):
            u, v = best['path_nodes'][i], best['path_nodes'][i+1]
            # Handle multi-edges (simplification: take first)
            edge_data = list(G.get_edge_data(u, v).values())[0]
            if 'geometry' in edge_data:
                edges.append(edge_data['geometry'])
            else:
                # If no geometry (straight line)
                edges.append(LineString([G.nodes[u]['x'], G.nodes[u]['y'], G.nodes[v]['x'], G.nodes[v]['y']])) # This is wrong, LineString needs 2 points
                # Correct fallback:
                p1 = Point(G.nodes[u]['x'], G.nodes[u]['y'])
                p2 = Point(G.nodes[v]['x'], G.nodes[v]['y'])
                edges.append(LineString([p1, p2]))
        
        # Combine edges into a single LineString or MultiLineString
        # Since we are in target_crs, the x,y from G.nodes are in target_crs? 
        # No, ox.graph_from_bbox uses lat/lon. We need to ensure graph is in target_crs.
        
        # Wait, I should project the graph first.
        pass

print("Wait, I need to re-structure the script to project the graph properly.")

import osmnx as ox
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point, LineString, MultiLineString, Polygon, MultiPolygon
import networkx as nx
import numpy as np

# 1. Load Incidents
incidents_df = pd.read_csv('incidents.csv')
incidents_gdf = gpd.GeoDataFrame(
    incidents_df, 
    geometry=gpd.points_from_xy(incidents_df.longitude, incidents_df.latitude),
    crs="EPSG:4326"
)

# Target CRS: Paris metric (EPSG:2154 - Lambert-93)
target_crs = "EPSG:2154"
incidents_gdf = incidents_gdf.to_crs(target_crs)

# Bounding box for OSM (from prompt) - note: lat/lon
bbox = (48.83, 2.30, 48.88, 2.38) # minlat, minlon, maxlat, maxlon

# 2. Fetch OSM Data
print("Fetching road network and hospitals...")
G = ox.graph_from_bbox(bbox[0], bbox[2], bbox[1], bbox[3], network_type='drive')

# Fetch hospitals
hospitals = ox.features_from_bbox(bbox[0], bbox[2], bbox[1], bbox[3], tags={'amenity': 'hospital'})

# Pre-process graph
G = ox.add_edge_speeds(G)
G = ox.add_edge_travel_times(G)

# Apply default 30 km/h if speed/travel_time is missing
# 30 km/h = 8.333 m/s
for u, v, k, data in G.edges(data=True, keys=True):
    if 'travel_time' not in data or np.isnan(data['travel_time']):
        length = data.get('length', 1.0)
        data['travel_time'] = length / 8.333333333333333

# Project the graph to target CRS
# This makes all node geometries and edge geometries metric
G_proj = ox.project_graph(G, to_crs=target_crs)

# Prepare hospital list
hospitals = hospitals.to_crs(target_crs)
# Use centroid for hospitals
hospitals['geometry'] = hospitals.geometry.centroid
# Ensure hospital names are unique and clean
hospitals['hospital_name'] = hospitals['name'].fillna('Unnamed Hospital')
# If duplicates exist in name, append index
hospitals['hospital_name'] = hospitals['hospital_name'].astype(str) + "_" + hospitals.index.astype(str)

# For distance, we need to find nearest nodes in projected graph
# Get node coords in target CRS
nodes_proj = ox.graph_to_gdfs(G_proj, nodes=True, edges=False)

def get_nearest_node_proj(graph_nodes, point):
    # Using the nodes dataframe to find nearest node
    # A simple approach with scipy or ox.nearest_nodes
    # Since graph_nodes is a GeoDataFrame, we can use it
    return ox.nearest_nodes(G_proj, point.x, point.y)

# 3. Analysis
print("Calculating routes and distances...")

hospitals['nearest_node'] = hospitals.geometry.apply(lambda p: get_nearest_node_proj(nodes_proj, p))

closest_hospital_rows = []
distance_matrix_rows = []

for _, incident in incidents_gdf.iterrows():
    inc_id = incident['incident_id']
    inc_node = get_nearest_node_proj(nodes_proj, incident.geometry)
    
    h_stats = []
    
    for _, hosp in hospitals.iterrows():
        h_node = hosp['nearest_node']
        h_name = hosp['hospital_name']
        
        try:
            # Network distance (length)
            dist = nx.shortest_path_length(G_proj, inc_node, h_node, weight='length')
            # Network time (seconds)
            time = nx.shortest_path_length(G_proj, inc_node, h_node, weight='travel_time')
            # Path for geometry
            path = nx.shortest_path(G_proj, inc_node, h_node, weight='length')
            
            h_stats.append({
                'hospital_name': h_name,
                'network_distance_m': dist,
                'travel_time_s': time,
                'path': path
            })
        except nx.NetworkXNoPath:
            continue
            
    # Sort by distance
    h_stats.sort(key=lambda x: x['network_distance_m'])
    
    # 1. Distance Matrix (up to 3)
    for rank, stat in enumerate(h_stats[:3], 1):
        distance_matrix_rows.append({
            'incident_id': inc_id,
            'hospital_name': stat['hospital_name'],
            'rank': rank,
            'network_distance_m': stat['network_distance_m']
        })
        
    # 2. Closest Hospital
    if h_stats:
        best = h_stats[0]
        # Create route geometry
        path_nodes = best['path']
        
        # Build route from edges
        route_edges = []
        for i in range(len(path_nodes) - 1):
            u, v = path_nodes[i], path_nodes[i+1]
            # Get edge geometry
            edge_data = G_proj.get_edge_data(u, v)
            if edge_data:
                # Pick the first edge if there are multiple (usually 1)
                # edge_data is a dict of dicts (k)
                k = list(edge_data.keys())[0]
                if 'geometry' in edge_data[k]:
                    route_edges.append(edge_data[k]['geometry'])
                else:
                    # Fallback to straight line
                    p1 = Point(G_proj.nodes[u]['x'], G_proj.nodes[u]['y'])
                    p2 = Point(G_proj.nodes[v]['x'], G_proj.nodes[v]['y'])
                    route_edges.append(LineString([p1, p2]))
        
        if route_edges:
            from shapely.ops import linemerge
            merged_route = linemerge(route_edges)
            closest_hospital_rows.append({
                'incident_id': inc_id,
                'hospital_name': best['hospital_name'],
                'network_distance_m': best['network_distance_m'],
                'geometry': merged_route
            })

# 4. Isochrones (15 min)
print("Calculating 15-min isochrones...")
isochrones_rows = []
# 15 minutes = 900 seconds
limit_s = 15 * 60

for _, hosp in hospitals.iterrows():
    h_node = hosp['nearest_node']
    h_name = hosp['hospital_name']
    
    # Get all nodes reachable within limit_s
    subgraph = nx.ego_graph(G_proj, h_node, radius=limit_s, distance='travel_time')
    
    if len(subgraph.nodes) > 0:
        # Get all edges in the subgraph
        sub_edges = subgraph.edges(data=True)
        edge_geoms = []
        for u, v, data in sub_edges:
            if 'geometry' in data:
                edge_geoms.append(data['geometry'])
            else:
                p1 = Point(G_proj.nodes[u]['x'], G_proj.nodes[u]['y'])
                p2 = Point(G_proj.nodes[v]['x'], G_proj.nodes[v]['y'])
                edge_geoms.append(LineString([p1, p2]))
        
        if edge_geoms:
            from shapely.ops import linemerge, unary_union
            merged_lines = linemerge(edge_geoms)
            # To create a polygon, we need to buffer the lines. 
            # A standard way to do isochrones is to buffer the edges, 
            # but that's not quite a "coverage area". 
            # However, for a "travel-time isochrone" in a network context, 
            # a common approach is the union of buffers of the reachable edges.
            # Let's buffer edges by something reasonable or just find the boundary.
            # Wait, a more accurate way is the concave hull or alpha shape of the nodes.
            # But buffering edges is a good approximation for "area reached".
            # Let's buffer edges by a small amount to ensure they connect, then unary_union.
            # Actually, let's buffer edges by say 20m.
            buffered_edges = unary_union([geom.buffer(20) for geom in edge_geoms])
            
            # For a cleaner isochrone, we can take the convex hull or a concave hull.
            # But the prompt asks for a "15-minute drive-time isochrone".
            # Let's use the union of buffered edges.
            isochrones_rows.append({
                'hospital_name': h_name,
                'travel_time_min': 15.0,
                'geometry': buffered_edges
            })
        else:
            # No edges? Just buffer the node
            p = Point(G_proj.nodes[h_node]['x'], G_proj.nodes[h_node]['y'])
            isochrones_rows.append({
                'hospital_name': h_name,
                'travel_time_min': 15.0,
                'geometry': p.buffer(50) # Small default
            })
    else:
        # No nodes reachable?
        p = Point(G_proj.nodes[h_node]['x'], G_proj.nodes[h_node]['y'])
        isochrones_rows.append({
            'hospital_name': h_name,
            'travel_time_min': 15.0,
            'geometry': p.buffer(50)
        })

# 5. Prepare Outputs
print("Saving outputs to GPKG...")

# incidents layer
closest_hospital_gdf = gpd.GeoDataFrame(closest_hospital_rows, crs=target_crs)
# Note: closest_hospital_gdf might have MultiLineString or LineString
# If merged_route is a LineString, it's fine.

# distance_matrix layer
dist_matrix_gdf = gpd.GeoDataFrame(distance_matrix_rows, crs=target_crs)
# The prompt says: distance_matrix layer is tabular (geometry may be empty or null)
# So we set geometry to None.
dist_matrix_gdf.set_geometry(None, inplace=True)

# isochrones layer
isochrones_gdf = gpd.GeoDataFrame(isochrones_rows, crs=target_crs)
# Ensure MultiPolygon/Polygon
isochrones_gdf['geometry'] = isochrones_gdf['geometry'].apply(lambda x: x if x.geom_type in ['Polygon', 'MultiPolygon'] else x.convex_hull)

# Final save
with pd.io.formats.format.DataFrameFormatter.precision_context(6):
    with gpd.GeoDataFrame(incidents_gdf, crs=target_crs).to_crs(target_crs).explode(index_parts=False) as inc_exploded:
        # Actually the prompt says incidents layer has Point (or MultiPoint)
        # We'll just use incidents_gdf directly.
        pass

# Let's do a more robust way to ensure geometry type for isochrones
# Some union results might be weird.
isochrones_gdf['geometry'] = isochrones_gdf['geometry'].buffer(0) # Clean up

# Create the GPKG
# Layers: incidents, closest_hospital, distance_matrix, isochrones_15min
incidents_gdf.to_file('emergency_routing.gpkg', layer='incidents', driver='GPKG')
closest_hospital_gdf.to_file('emergency_routing.gpkg', layer='closest_hospital', driver='GPKG')
dist_matrix_gdf.to_file('emergency_routing.gpkg', layer='distance_matrix', driver='GPKG')
isochrones_gdf.to_file('emergency_routing.gpkg', layer='isochrones_15min', driver='GPKG')

print("Done!")

import osmnx as ox
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point, LineString, MultiLineString, Polygon, MultiPolygon
import networkx as nx
import numpy as np
from shapely.ops import linemerge, unary_union

# 1. Load Incidents
incidents_df = pd.read_csv('incidents.csv')
incidents_gdf = gpd.GeoDataFrame(
    incidents_df, 
    geometry=gpd.points_from_xy(incidents_df.longitude, incidents_df.latitude),
    crs="EPSG:4326"
)

# Target CRS: Paris metric (EPSG:2154 - Lambert-93)
target_crs = "EPSG:2154"
incidents_gdf = incidents_gdf.to_crs(target_crs)

# Bounding box for OSM (left, bottom, right, top)
bbox = (2.30, 48.83, 2.38, 48.88) 

# 2. Fetch OSM Data
print("Fetching road network and hospitals...")
G = ox.graph_from_bbox(bbox, network_type='drive')
hospitals = ox.features_from_bbox(bbox, tags={'amenity': 'hospital'})

# Pre-process graph
G = ox.add_edge_speeds(G)
G = ox.add_edge_travel_times(G)

# Apply default 30 km/h if speed/travel_time is missing
for u, v, k, data in G.edges(data=True, keys=True):
    if 'travel_time' not in data or np.isnan(data['travel_time']):
        length = data.get('length', 1.0)
        data['travel_time'] = length / 8.333333333333333

# Project the graph to target CRS
G_proj = ox.project_graph(G, to_crs=target_crs)

# Prepare hospital list
hospitals = hospitals.to_crs(target_crs)
hospitals['geometry'] = hospitals.geometry.centroid
hospitals['hospital_name'] = hospitals['name'].fillna('Unnamed Hospital').astype(str) + "_" + hospitals.index.astype(str)

# Get node coords in target CRS
nodes_proj = ox.graph_to_gdfs(G_proj, nodes=True, edges=False)

def get_nearest_node_proj(graph_nodes, point):
    return ox.nearest_nodes(G_proj, point.x, point.y)

# 3. Analysis
print("Calculating routes and distances...")
hospitals['nearest_node'] = hospitals.geometry.apply(lambda p: get_nearest_node_proj(nodes_proj, p))

closest_hospital_rows = []
distance_matrix_rows = []

for _, incident in incidents_gdf.iterrows():
    inc_id = incident['incident_id']
    inc_node = get_nearest_node_proj(nodes_proj, incident.geometry)
    
    h_stats = []
    for _, hosp in hospitals.iterrows():
        h_node = hosp['nearest_node']
        h_name = hosp['hospital_name']
        try:
            dist = nx.shortest_path_length(G_proj, inc_node, h_node, weight='length')
            time = nx.shortest_path_length(G_proj, inc_node, h_node, weight='travel_time')
            path = nx.shortest_path(G_proj, inc_node, h_node, weight='length')
            h_stats.append({'hospital_name': h_name, 'network_distance_m': dist, 'travel_time_s': time, 'path': path})
        except nx.NetworkXNoPath:
            continue
            
    h_stats.sort(key=lambda x: x['network_distance_m'])
    
    # 1. Distance Matrix (up to 3)
    for rank, stat in enumerate(h_stats[:3], 1):
        distance_matrix_rows.append({
            'incident_id': inc_id,
            'hospital_name': stat['hospital_name'],
            'rank': rank,
            'network_distance_m': stat['network_distance_m']
        })
        
    # 2. Closest Hospital
    if h_stats:
        best = h_stats[0]
        path_nodes = best['path']
        route_edges = []
        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:
                k = list(edge_data.keys())[0]
                if 'geometry' in edge_data[k]:
                    route_edges.append(edge_data[k]['geometry'])
                else:
                    p1 = Point(G_proj.nodes[u]['x'], G_proj.nodes[u]['y'])
                    p2 = Point(G_proj.nodes[v]['x'], G_proj.nodes[v]['y'])
                    route_edges.append(LineString([p1, p2]))
        if route_edges:
            merged_route = linemerge(route_edges)
            closest_hospital_rows.append({
                'incident_id': inc_id,
                'hospital_name': best['hospital_name'],
                'network_distance_m': best['network_distance_m'],
                'geometry': merged_route
            })

# 4. Isochrones (15 min)
print("Calculating 15-min isochrones...")
isochrones_rows = []
limit_s = 15 * 60

for _, hosp in hospitals.iterrows():
    h_node = hosp['nearest_node']
    h_name = hosp['hospital_name']
    subgraph = nx.ego_graph(G_proj, h_node, radius=limit_s, distance='travel_time')
    
    if len(subgraph.nodes) > 0:
        sub_edges = subgraph.edges(data=True)
        edge_geoms = []
        for u, v, data in sub_edges:
            if 'geometry' in data:
                edge_geoms.append(data['geometry'])
            else:
                p1 = Point(G_proj.nodes[u]['x'], G_proj.nodes[u]['y'])
                p2 = Point(G_proj.nodes[v]['x'], G_proj.nodes[v]['y'])
                edge_geoms.append(LineString([p1, p2]))
        if edge_geoms:
            buffered_edges = unary_union([geom.buffer(20) for geom in edge_geoms])
            isochrones_rows.append({'hospital_name': h_name, 'travel_time_min': 15.0, 'geometry': buffered_edges})
        else:
            p = Point(G_proj.nodes[h_node]['x'], G_proj.nodes[h_node]['y'])
            isochrones_rows.append({'hospital_name': h_name, 'travel_time_min': 15.0, 'geometry': p.buffer(50)})
    else:
        p = Point(G_proj.nodes[h_node]['x'], G_proj.nodes[h_node]['y'])
        isochrones_rows.append({'hospital_name': h_name, 'travel_time_min': 15.0, 'geometry': p.buffer(50)})

# 5. Prepare Outputs
print("Saving outputs to GPKG...")
closest_hospital_gdf = gpd.GeoDataFrame(closest_hospital_rows, crs=target_crs)
dist_matrix_gdf = gpd.GeoDataFrame(distance_matrix_rows, crs=target_crs)
dist_matrix_gdf.set_geometry(None, inplace=True)
isochrones_gdf = gpd.GeoDataFrame(isochrones_rows, crs=target_crs)
isochrones_gdf['geometry'] = isochrones_gdf['geometry'].apply(lambda x: x if x.geom_type in ['Polygon', 'MultiPolygon'] else x.convex_hull)
isochrones_gdf['geometry'] = isochrones_gdf['geometry'].buffer(0)

incidents_gdf.to_file('emergency_routing.gpkg', layer='incidents', driver='GPKG')
closest_hospital_gdf.to_file('emergency_routing.gpkg', layer='closest_hospital', driver='GPKG')
dist_matrix_gdf.to_file('emergency_routing.gpkg', layer='distance_matrix', driver='GPKG')
isochrones_gdf.to_file('emergency_routing.gpkg', layer='isochrones_15min', driver='GPKG')

print("Done!")

import osmnx as ox
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point, LineString, MultiLineString, Polygon, MultiPolygon
import networkx as nx
import numpy as np
from shapely.ops import linemerge, unary_union

# 1. Load Incidents
incidents_df = pd.read_csv('incidents.csv')
incidents_gdf = gpd.GeoDataFrame(
    incidents_df, 
    geometry=gpd.points_from_xy(incidents_df.longitude, incidents_df.latitude),
    crs="EPSG:4326"
)

# Target CRS: Paris metric (EPSG:2154 - Lambert-93)
target_crs = "EPSG:2154"
incidents_gdf = incidents_gdf.to_crs(target_crs)

# Bounding box for OSM (left, bottom, right, top)
bbox = (2.30, 48.83, 2.38, 48.88) 

# 2. Fetch OSM Data
print("Fetching road network and hospitals...")
G = ox.graph_from_bbox(bbox=bbox, network_type='drive')
hospitals = ox.features_from_bbox(bbox=bbox, tags={'amenity': 'hospital'})

# Pre-process graph
G = ox.add_edge_speeds(G)
G = ox.add_edge_travel_times(G)

# Apply default 30 km/h if speed/travel_time is missing
for u, v, k, data in G.edges(data=True, keys=True):
    if 'travel_time' not in data or np.isnan(data['travel_time']):
        length = data.get('length', 1.0)
        data['travel_time'] = length / 8.333333333333333

# Project the graph to target CRS
G_proj = ox.project_graph(G, to_crs=target_crs)

# Prepare hospital list
hospitals = hospitals.to_crs(target_crs)
hospitals['geometry'] = hospitals.geometry.centroid
# Ensure hospital names are unique and clean
hospitals['hospital_name'] = hospitals['name'].fillna('Unnamed Hospital').astype(str) + "_" + hospitals.index.astype(str)

# Get node coords in target CRS
nodes_proj = ox.graph_to_gdfs(G_proj, nodes=True, edges=False)

def get_nearest_node_proj(graph_nodes, point):
    return ox.nearest_nodes(G_proj, point.x, point.y)

# 3. Analysis
print("Calculating routes and distances...")
hospitals['nearest_node'] = hospitals.geometry.apply(lambda p: get_nearest_node_proj(nodes_proj, p))

closest_hospital_rows = []
distance_matrix_rows = []

for _, incident in incidents_gdf.iterrows():
    inc_id = incident['incident_id']
    inc_node = get_nearest_node_proj(nodes_proj, incident.geometry)
    
    h_stats = []
    for _, hosp in hospitals.iterrows():
        h_node = hosp['nearest_node']
        h_name = hosp['hospital_name']
        try:
            dist = nx.shortest_path_length(G_proj, inc_node, h_node, weight='length')
            time = nx.shortest_path_length(G_proj, inc_node, h_node, weight='travel_time')
            path = nx.shortest_path(G_proj, inc_node, h_node, weight='length')
            h_stats.append({'hospital_name': h_name, 'network_distance_m': dist, 'travel_time_s': time, 'path': path})
        except nx.NetworkXNoPath:
            continue
            
    h_stats.sort(key=lambda x: x['network_distance_m'])
    
    # 1. Distance Matrix (up to 3)
    for rank, stat in enumerate(h_stats[:3], 1):
        distance_matrix_rows.append({
            'incident_id': inc_id,
            'hospital_name': stat['hospital_name'],
            'rank': rank,
            'network_distance_m': stat['network_distance_m']
        })
        
    # 2. Closest Hospital
    if h_stats:
        best = h_stats[0]
        path_nodes = best['path']
        route_edges = []
        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:
                k = list(edge_data.keys())[0]
                if 'geometry' in edge_data[k]:
                    route_edges.append(edge_data[k]['geometry'])
                else:
                    p1 = Point(G_proj.nodes[u]['x'], G_proj.nodes[u]['y'])
                    p2 = Point(G_proj.nodes[v]['x'], G_proj.nodes[v]['y'])
                    route_edges.append(LineString([p1, p2]))
        if route_edges:
            merged_route = linemerge(route_edges)
            closest_hospital_rows.append({
                'incident_id': inc_id,
                'hospital_name': best['hospital_name'],
                'network_distance_m': best['network_distance_m'],
                'geometry': merged_route
            })

# 4. Isochrones (15 min)
print("Calculating 15-min isochrones...")
isochrones_rows = []
limit_s = 15 * 60

for _, hosp in hospitals.iterrows():
    h_node = hosp['nearest_node']
    h_name = hosp['hospital_name']
    subgraph = nx.ego_graph(G_proj, h_node, radius=limit_s, distance='travel_time')
    
    if len(subgraph.nodes) > 0:
        sub_edges = subgraph.edges(data=True)
        edge_geoms = []
        for u, v, data in sub_edges:
            if 'geometry' in data:
                edge_geoms.append(data['geometry'])
            else:
                p1 = Point(G_proj.nodes[u]['x'], G_proj.nodes[u]['y'])
                p2 = Point(G_proj.nodes[v]['x'], G_proj.nodes[v]['y'])
                edge_geoms.append(LineString([p1, p2]))
        if edge_geoms:
            buffered_edges = unary_union([geom.buffer(20) for geom in edge_geoms])
            isochrones_rows.append({'hospital_name': h_name, 'travel_time_min': 15.0, 'geometry': buffered_edges})
        else:
            p = Point(G_proj.nodes[h_node]['x'], G_proj.nodes[h_node]['y'])
            isochrones_rows.append({'hospital_name': h_name, 'travel_time_min': 15.0, 'geometry': p.buffer(50)})
    else:
        p = Point(G_proj.nodes[h_node]['x'], G_proj.nodes[h_node]['y'])
        isochrones_rows.append({'hospital_name': h_name, 'travel_time_min': 15.0, 'geometry': p.buffer(50)})

# 5. Prepare Outputs
print("Saving outputs to GPKG...")
closest_hospital_gdf = gpd.GeoDataFrame(closest_hospital_rows, crs=target_crs)
dist_matrix_gdf = gpd.GeoDataFrame(distance_matrix_rows, crs=target_crs)
dist_matrix_gdf.set_geometry(None, inplace=True)
isochrones_gdf = gpd.GeoDataFrame(isochrones_rows, crs=target_crs)
isochrones_gdf['geometry'] = isochrones_gdf['geometry'].apply(lambda x: x if x.geom_type in ['Polygon', 'MultiPolygon'] else x.convex_hull)
isochrones_gdf['geometry'] = isochrones_gdf['geometry'].buffer(0)

incidents_gdf.to_file('emergency_routing.gpkg', layer='incidents', driver='GPKG')
closest_hospital_gdf.to_file('emergency_routing.gpkg', layer='closest_hospital', driver='GPKG')
dist_matrix_gdf.to_file('emergency_routing.gpkg', layer='distance_matrix', driver='GPKG')
isochrones_gdf.to_file('emergency_routing.gpkg', layer='isochrones_15min', driver='GPKG')

print("Done!")