spa-l3-paris-emergency-routing

Rerun

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

# 1. Setup Area and Data
bbox = (48.83, 2.30, 48.88, 2.38)  # min_lat, min_lon, max_lat, max_lon
# osmnx bbox is (north, south, east, west)
ox_bbox = (bbox[2], bbox[0], bbox[3], bbox[1])

print("Downloading road network...")
G = ox.graph_from_bbox(ox_bbox[0], ox_bbox[1], ox_bbox[2], ox_bbox[3], network_type='drive')

print("Downloading hospitals...")
# Using ox.geometries_from_bbox or similar. 
# In newer osmnx, use ox.features_from_bbox
hospitals_gdf = ox.features_from_bbox(ox_bbox[0], ox_bbox[1], ox_bbox[2], ox_bbox[3], tags={'amenity': 'hospital'})

# Clean hospital data
hospitals_gdf = hospitals_gdf[hospitals_gdf['amenity'] == 'hospital'].copy()
hospitals_gdf['hospital_name'] = hospitals_gdf['name'].fillna('Unnamed Hospital')
hospitals_gdf = hospitals_gdf.reset_index()

print(f"Found {len(hospitals_gdf)} hospitals.")

# 2. Prepare Road Network for Routing
# Apply speed limits
# osmnx handles this, but let's ensure speed/travel_time are set
G = ox.add_edge_speeds(G)
G = ox.add_edge_travel_times(G)

# Verify/Set default speed if missing
for u, v, k, data in G.edges(data=True, keys=True):
    if 'speed_kph' not in data or pd.isna(data['speed_kph']):
        data['speed_kph'] = 30.0
    if 'travel_time' not in data or pd.isna(data['travel_time']):
        # travel_time is in seconds. time = distance / speed
        # speed is in km/h. speed_mps = speed_kph / 3.6
        dist = data.get('length', 0)
        speed_mps = data['speed_kph'] / 3.6
        if speed_mps > 0:
            data['travel_time'] = dist / speed_mps
        else:
            data['travel_time'] = 0

# 3. 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"
)

# 4. Coordinate System: Paris metric (EPSG:2154 - Lambert-93)
# We'll project everything to EPSG:2154
G_proj = ox.project_graph(G, to_crs='EPSG:2154')
incidents_proj = incidents_gdf.to_crs('EPSG:2154')
hospitals_proj = hospitals_gdf.to_crs('EPSG:2154')

# Get node coordinates for hospitals
hospitals_nodes = []
for idx, row in hospitals_proj.iterrows():
    # Use centroid for hospital location if it's a polygon
    centroid = row.geometry.centroid
    hospitals_nodes.append({
        'hospital_name': row['hospital_name'],
        'geometry': centroid,
        'id': idx
    })
hospitals_nodes_gdf = gpd.GeoDataFrame(hospitals_nodes, crs='EPSG:2154')

# 5. Routing logic
# We need to find the nearest node in the projected graph for each incident and hospital
incidents_nodes = ox.nearest_nodes(G_proj, incidents_proj.geometry.x, incidents_proj.geometry.y)
hospitals_nodes_idx = ox.nearest_nodes(G_proj, hospitals_nodes_gdf.geometry.x, hospitals_nodes_gdf.geometry.y)

# results containers
closest_hospital_rows = []
distance_matrix_rows = []

print("Calculating routes...")
for i, (idx, incident) in enumerate(incidents_proj.iterrows()):
    orig_node = incidents_nodes[i]
    
    # Distances to all hospitals
    dist_to_hospitals = []
    
    for h_idx, h_node in enumerate(hospitals_nodes_idx):
        try:
            # shortest path by weight 'length'
            path = nx.shortest_path(G_proj, orig_node, h_node, weight='length')
            path_length = nx.shortest_path_length(G_proj, orig_node, h_node, weight='length')
            
            # Build LineString for the route
            nodes_coords = [(G_proj.nodes[n]['x'], G_proj.nodes[n]['y']) for n in path]
            route_geom = LineString(nodes_coords)
            
            dist_to_hospitals.append({
                'hospital_name': hospitals_nodes_gdf.iloc[h_idx]['hospital_name'],
                'network_distance_m': path_length,
                'route_geom': route_geom
            })
        except nx.NetworkXNoPath:
            continue

    if not dist_to_hospitals:
        continue
        
    # Sort by distance
    dist_to_hospitals.sort(key=lambda x: x['network_distance_m'])
    
    # Closest Hospital (Layer 2)
    closest = dist_to_hospitals[0]
    closest_hospital_rows.append({
        'incident_id': incident['incident_id'],
        'hospital_name': closest['hospital_name'],
        'network_distance_m': closest['network_distance_m'],
        'geometry': closest['route_geom']
    })
    
    # Distance Matrix (Layer 3)
    for rank, h_info in enumerate(dist_to_hospitals[:3], 1):
        distance_matrix_rows.append({
            'incident_id': incident['incident_id'],
            'hospital_name': h_info['hospital_name'],
            'rank': rank,
            'network_distance_m': h_info['network_distance_m'],
            'geometry': None
        })

closest_hospital_gdf = gpd.GeoDataFrame(closest_hospital_rows, crs='EPSG:2154')
distance_matrix_gdf = gpd.GeoDataFrame(distance_matrix_rows, crs='EPSG:2154')

# 6. Isochrones (Layer 4)
# 15 minutes = 900 seconds
print("Calculating isochrones...")
isochrones_rows = []
for idx, h_row in hospitals_nodes_gdf.iterrows():
    h_node = hospitals_nodes_idx[idx]
    
    # Use Dijkstra to find all nodes reachable within 900s
    try:
        subgraph = nx.single_source_dijkstra_path_length(G_proj, h_node, cutoff=900, weight='travel_time')
        
        # Get the nodes and edges
        reachable_nodes = list(subgraph.keys())
        if not reachable_nodes:
            continue
            
        # Collect all edges in the subgraph
        edges_list = []
        for u, v, k, data in G_proj.edges(data=True, keys=True):
            if u in subgraph and v in subgraph:
                # Check if both ends are within the time limit (travel_time is one-way, but for simplicity...)
                # Actually, single_source_dijkstra finds nodes reachable from h_node.
                # The edge (u,v) is part of a path from h_node to v if v is in subgraph and path(h,u)+travel(u,v) <= 900
                # But for simplicity, if both u and v are in reachable nodes, the edge is likely part of the reachable set.
                # To be precise, we should check the time to the further node.
                time_to_u = subgraph[u]
                time_to_v = subgraph[v]
                # Since it's a directed graph, we check if u is reachable from h_node
                # and the edge (u,v) helps reach v.
                if time_to_u + data['travel_time'] <= 900:
                    edges_list.append((u, v, data))
                elif time_to_v + data['travel_time'] <= 900: # (Wait, this is for reversed edges, but let's stick to directed)
                     # We'll only use edges that are part of the paths from the source.
                     pass
        
        # A better way: extract the subgraph of edges where the target node is within cutoff
        # We'll use the property that if v is in subgraph, the edge (u,v) that was used to reach it is valid.
        valid_edges = []
        for v in subgraph:
            # We need to find which edge brought us here. Dijkstra doesn't give edges directly.
            # Let's use nx.ego_graph
            pass
        
        # Re-try:
        subgraph_nodes = list(subgraph.keys())
        # To get the edges, we use nx.ego_graph or filter G_proj
        # Ego graph uses distance, we need travel_time.
        # Let's use the nodes found and filter edges where travel_time from source is <= 900
        
        # Simplified approach: use nx.ego_graph with weight='travel_time'
        ego_subgraph = nx.ego_graph(G_proj, h_node, radius=900, distance='travel_time')
        
        # Create geometry from edges in ego_subgraph
        edge_geoms = []
        for u, v, k, data in ego_subgraph.edges(data=True, keys=True):
            # Get coords
            x1, y1 = G_proj.nodes[u]['x'], G_proj.nodes[u]['y']
            x2, y2 = G_proj.nodes[v]['x'], G_proj.nodes[v]['y']
            edge_geoms.append(LineString([(x1, y1), (x2, y2)]))
        
        if edge_geoms:
            # Combine edges into a MultiLineString and then buffer it to create a polygon
            # or just unary_union them and buffer.
            combined_lines = MultiLineString(edge_geoms).unary_union
            # Buffer creates a polygon around the lines. 
            # However, for isochrones, the "area" is often approximated by buffering the lines.
            # Since this is a road network, buffering the lines (the roads) is a common way.
            # But the instruction asks for "isochrone" which is the area.
            # A better way: concave hull of the nodes or buffer the lines and union them.
            # Let's buffer the lines by a small amount and union, or use the nodes' buffer.
            # Actually, the most common way in OSMnx is to buffer the lines.
            # But we want the "service area". Let's buffer the lines with a certain width 
            # to represent the road area and then buffer the result.
            # Or, even better, just use the union of buffers around each node.
            # Let's try union of buffers around nodes.
            node_geoms = [Point(G_proj.nodes[n]['x'], G_proj.nodes[n]['y']).buffer(50) for n in subgraph_nodes]
            isochrone_poly = gpd.GeoSeries(node_geoms).unary_union.convex_hull # convex_hull is too simple
            # Let's use the union of buffers of the lines.
            poly = combined_lines.buffer(50) # 50m buffer around roads
            # To make it more like an area, we can buffer it more or use the union of node buffers.
            # Let's use the union of node buffers and then a small buffer to close gaps.
            isochrone_poly = gpd.GeoSeries(node_geoms).unary_union.buffer(100) # buffer the union of node points
            
            isochrones_rows.append({
                'hospital_name': h_row['hospital_name'],
                'travel_time_min': 15.0,
                'geometry': isochrone_poly
            })
        else:
            # If no edges, maybe just a point buffer?
            isochrones_rows.append({
                'hospital_name': h_row['hospital_name'],
                'travel_time_min': 15.0,
                'geometry': h_row['geometry'].buffer(100)
            })
            
    except Exception as e:
        print(f"Error calculating isochrone for {h_row['hospital_name']}: {e}")
        continue

isochrones_gdf = gpd.GeoDataFrame(isochrones_rows, crs='EPSG:2154')

# 7. Final Layers and Save
# Incidents layer
incidents_final = incidents_proj.copy()
incidents_final = incidents_final[['incident_id', 'geometry']]

# Closest Hospital layer
closest_hospital_gdf = closest_hospital_gdf[['incident_id', 'hospital_name', 'network_distance_m', 'geometry']]

# Distance Matrix layer
distance_matrix_gdf = distance_matrix_gdf[['incident_id', 'hospital_name', 'rank', 'network_distance_m', 'geometry']]

# Isochrones layer
isochrones_gdf = isochrones_gdf[['hospital_name', 'travel_time_min', 'geometry']]

# Save to GPKG
print("Saving to emergency_routing.gpkg...")
with pd.io.formats.formatters.BaseFormatter.print_version: # dummy
    import fiona
    with fiona.open("/work/emergency_routing.gpkg", 'w', driver='GPKG') as dst:
        # This is not the right way to save multiple layers with geopandas
        pass

# The correct way to save multiple layers in GPKG with geopandas:
# We'll save them one by one using to_file with mode='a' (append)
# First, create the file with the first layer
incidents_final.to_file("/work/emergency_routing.gpkg", layer='incidents', driver='GPKG')
closest_hospital_gdf.to_file("/work/emergency_routing.gpkg", layer='closest_hospital', driver='GPKG', mode='a')
distance_matrix_gdf.to_file("/work/emergency_routing.gpkg", layer='distance_matrix', driver='GPKG', mode='a')
isochrones_gdf.to_file("/work/emergency_routing.gpkg", layer='isochrones_15min', driver='GPKG', mode='a')

print("Done!")

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

# 1. Setup Area and Data
# bbox: 48.83,2.30 to 48.88,2.38 (min_lat, min_lon, max_lat, max_lon)
north, south, east, west = 48.88, 48.83, 2.38, 2.30

print("Downloading road network...")
# OSMnx graph_from_bbox expects (north, south, east, west)
G = ox.graph_from_bbox(north=north, south=south, east=east, west=west, network_type='drive')

print("Downloading hospitals...")
hospitals_gdf = ox.features_from_bbox(north=north, south=south, east=east, west=west, tags={'amenity': 'hospital'})

# Clean hospital data
hospitals_gdf = hospitals_gdf[hospitals_gdf['amenity'] == 'hospital'].copy()
hospitals_gdf['hospital_name'] = hospitals_gdf['name'].fillna('Unnamed Hospital')
hospitals_gdf = hospitals_gdf.reset_index()

print(f"Found {len(hospitals_gdf)} hospitals.")

# 2. Prepare Road Network for Routing
G = ox.add_edge_speeds(G)
G = ox.add_edge_travel_times(G)

# Ensure speed and travel_time are set correctly
for u, v, k, data in G.edges(data=True, keys=True):
    if 'speed_kph' not in data or pd.isna(data['speed_kph']):
        data['speed_kph'] = 30.0
    if 'travel_time' not in data or pd.isna(data['travel_time']):
        dist = data.get('length', 0)
        speed_mps = data['speed_kph'] / 3.6
        if speed_mps > 0:
            data['travel_time'] = dist / speed_mps
        else:
            data['travel_time'] = 0

# 3. 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"
)

# 4. Coordinate System: Paris metric (EPSG:2154 - Lambert-93)
G_proj = ox.project_graph(G, to_crs='EPSG:2154')
incidents_proj = incidents_gdf.to_crs('EPSG:2154')
hospitals_proj = hospitals_gdf.to_crs('EPSG:2154')

# Get centroid for hospital locations
hospitals_nodes_gdf = hospitals_proj.copy()
hospitals_nodes_gdf['geometry'] = hospitals_nodes_gdf.geometry.centroid

print(f"Routing with {len(hospitals_nodes_gdf)} hospitals and {len(incidents_proj)} incidents...")

# 5. Routing logic
# Get nearest nodes for incidents and hospitals
# For incidents
incidents_nodes = []
for _, row in incidents_proj.iterrows():
    node = ox.nearest_nodes(G_proj, row.geometry.x, row.geometry.y)
    incidents_nodes.append(node)

# For hospitals
hospitals_nodes_idx = []
for _, row in hospitals_nodes_gdf.iterrows():
    node = ox.nearest_nodes(G_proj, row.geometry.x, row.geometry.y)
    hospitals_nodes_idx.append(node)

closest_hospital_rows = []
distance_matrix_rows = []

print("Calculating routes...")
for i, (idx, incident) in enumerate(incidents_proj.iterrows()):
    orig_node = incidents_nodes[i]
    
    dist_to_hospitals = []
    
    for h_idx, h_node in enumerate(hospitals_nodes_idx):
        try:
            # shortest path by weight 'length'
            path_length = nx.shortest_path_length(G_proj, orig_node, h_node, weight='length')
            path = nx.shortest_path(G_proj, orig_node, h_node, weight='length')
            
            # Build LineString for the route
            nodes_coords = [(G_proj.nodes[n]['x'], G_proj.nodes[n]['y']) for n in path]
            route_geom = LineString(nodes_coords)
            
            dist_to_hospitals.append({
                'hospital_name': hospitals_nodes_gdf.iloc[h_idx]['hospital_name'],
                'network_distance_m': path_length,
                'route_geom': route_geom
            })
        except (nx.NetworkXNoPath, nx.NodeNotFound):
            continue

    if not dist_to_hospitals:
        continue
        
    # Sort by distance
    dist_to_hospitals.sort(key=lambda x: x['network_distance_m'])
    
    # Closest Hospital (Layer 2)
    closest = dist_to_hospitals[0]
    closest_hospital_rows.append({
        'incident_id': incident['incident_id'],
        'hospital_name': closest['hospital_name'],
        'network_distance_m': closest['network_distance_m'],
        'geometry': closest['route_geom']
    })
    
    # Distance Matrix (Layer 3)
    for rank, h_info in enumerate(dist_to_hospitals[:3], 1):
        distance_matrix_rows.append({
            'incident_id': incident['incident_id'],
            'hospital_name': h_info['hospital_name'],
            'rank': rank,
            'network_distance_m': h_info['network_distance_m'],
            'geometry': None
        })

closest_hospital_gdf = gpd.GeoDataFrame(closest_hospital_rows, crs='EPSG:2154')
distance_matrix_gdf = gpd.GeoDataFrame(distance_matrix_rows, crs='EPSG:2154')

# 6. Isochrones (Layer 4)
# 15 minutes = 900 seconds
print("Calculating isochrones...")
isochrones_rows = []
for idx, h_row in hospitals_nodes_gdf.iterrows():
    h_node = hospitals_nodes_idx[idx]
    
    try:
        # Find nodes reachable within 900s
        subgraph_nodes = nx.single_source_dijkstra_path_length(G_proj, h_node, cutoff=900, weight='travel_time')
        
        if not subgraph_nodes:
            continue
            
        # Extract reachable nodes' geometry
        node_geoms = []
        for n in subgraph_nodes.keys():
            node_geoms.append(Point(G_proj.nodes[n]['x'], G_proj.nodes[n]['y']))
        
        if node_geoms:
            # Create the isochrone polygon by buffering the union of reachable nodes
            # A small buffer on the union of nodes is a common approximation
            isochrone_poly = gpd.GeoSeries(node_geoms).unary_union.buffer(150)
            
            isochrones_rows.append({
                'hospital_name': h_row['hospital_name'],
                'travel_time_min': 15.0,
                'geometry': isochrone_poly
            })
            
    except Exception as e:
        print(f"Error calculating isochrone for {h_row['hospital_name']}: {e}")
        continue

isochrones_gdf = gpd.GeoDataFrame(isochrones_rows, crs='EPSG:2154')

# 7. Final Layers and Save
# Incidents layer
incidents_final = incidents_proj[['incident_id', 'geometry']].copy()

# Closest Hospital layer
closest_hospital_final = closest_hospital_gdf[['incident_id', 'hospital_name', 'network_distance_m', 'geometry']]

# Distance Matrix layer
distance_matrix_final = distance_matrix_gdf[['incident_id', 'hospital_name', 'rank', 'network_distance_m', 'geometry']]

# Isochrones layer
isochrones_final = isochrones_gdf[['hospital_name', 'travel_time_min', 'geometry']]

# Save to GPKG
print("Saving to emergency_routing.gpkg...")
incidents_final.to_file("/work/emergency_routing.gpkg", layer='incidents', driver='GPKG')
closest_hospital_final.to_file("/work/emergency_routing.gpkg", layer='closest_hospital', driver='GPKG', mode='a')
distance_matrix_final.to_file("/work/emergency_routing.gpkg", layer='distance_matrix', driver='GPKG', mode='a')
isochrones_final.to_file("/work/emergency_routing.gpkg", layer='isochrones_15min', driver='GPKG', mode='a')

print("Done!")

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

# 1. Setup Area and Data
# bbox: 48.83,2.30 to 48.88,2.38 (min_lat, min_lon, max_lat, max_lon)
# help(ox.graph_from_bbox) says bbox is (left, bottom, right, top) -> (west, south, east, north)
west, south, east, north = 2.30, 48.83, 2.38, 48.88
bbox_tuple = (west, south, east, north)

print(f"Downloading road network for bbox {bbox_tuple}...")
G = ox.graph_from_bbox(bbox=bbox_tuple, network_type='drive')

print("Downloading hospitals...")
# features_from_bbox also takes bbox tuple (left, bottom, right, top)
hospitals_gdf = ox.features_from_bbox(bbox=bbox_tuple, tags={'amenity': 'hospital'})

# Clean hospital data
hospitals_gdf = hospitals_gdf[hospitals_gdf['amenity'] == 'hospital'].copy()
hospitals_gdf['hospital_name'] = hospitals_gdf['name'].fillna('Unnamed Hospital')
hospitals_gdf = hospitals_gdf.reset_index()

print(f"Found {len(hospitals_gdf)} hospitals.")

# 2. Prepare Road Network for Routing
G = ox.add_edge_speeds(G)
G = ox.add_edge_travel_times(G)

# Ensure speed and travel_time are set correctly
for u, v, k, data in G.edges(data=True, keys=True):
    if 'speed_kph' not in data or pd.isna(data['speed_kph']):
        data['speed_kph'] = 30.0
    if 'travel_time' not in data or pd.isna(data['travel_time']):
        dist = data.get('length', 0)
        speed_mps = data['speed_kph'] / 3.6
        if speed_mps > 0:
            data['travel_time'] = dist / speed_mps
        else:
            data['travel_time'] = 0

# 3. 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"
)

# 4. Coordinate System: Paris metric (EPSG:2154 - Lambert-93)
G_proj = ox.project_graph(G, to_crs='EPSG:2154')
incidents_proj = incidents_gdf.to_crs('EPSG:2154')
hospitals_proj = hospitals_gdf.to_crs('EPSG:2154')

# Get centroid for hospital locations
hospitals_nodes_gdf = hospitals_proj.copy()
hospitals_nodes_gdf['geometry'] = hospitals_nodes_gdf.geometry.centroid

print(f"Routing with {len(hospitals_nodes_gdf)} hospitals and {len(incidents_proj)} incidents...")

# 5. Routing logic
# Get nearest nodes for incidents and hospitals
incidents_nodes = []
for _, row in incidents_proj.iterrows():
    node = ox.nearest_nodes(G_proj, row.geometry.x, row.geometry.y)
    incidents_nodes.append(node)

hospitals_nodes_idx = []
for _, row in hospitals_nodes_gdf.iterrows():
    node = ox.nearest_nodes(G_proj, row.geometry.x, row.geometry.y)
    hospitals_nodes_idx.append(node)

closest_hospital_rows = []
distance_matrix_rows = []

print("Calculating routes...")
for i, (idx, incident) in enumerate(incidents_proj.iterrows()):
    orig_node = incidents_nodes[i]
    
    dist_to_hospitals = []
    
    for h_idx, h_node in enumerate(hospitals_nodes_idx):
        try:
            # shortest path by weight 'length'
            path_length = nx.shortest_path_length(G_proj, orig_node, h_node, weight='length')
            path = nx.shortest_path(G_proj, orig_node, h_node, weight='length')
            
            # Build LineString for the route
            nodes_coords = [(G_proj.nodes[n]['x'], G_proj.nodes[n]['y']) for n in path]
            route_geom = LineString(nodes_coords)
            
            dist_to_hospitals.append({
                'hospital_name': hospitals_nodes_gdf.iloc[h_idx]['hospital_name'],
                'network_distance_m': path_length,
                'route_geom': route_geom
            })
        except (nx.NetworkXNoPath, nx.NodeNotFound):
            continue

    if not dist_to_hospitals:
        continue
        
    # Sort by distance
    dist_to_hospitals.sort(key=lambda x: x['network_distance_m'])
    
    # Closest Hospital (Layer 2)
    closest = dist_to_hospitals[0]
    closest_hospital_rows.append({
        'incident_id': incident['incident_id'],
        'hospital_name': closest['hospital_name'],
        'network_distance_m': closest['network_distance_m'],
        'geometry': closest['route_geom']
    })
    
    # Distance Matrix (Layer 3)
    for rank, h_info in enumerate(dist_to_hospitals[:3], 1):
        distance_matrix_rows.append({
            'incident_id': incident['incident_id'],
            'hospital_name': h_info['hospital_name'],
            'rank': rank,
            'network_distance_m': h_info['network_distance_m'],
            'geometry': None
        })

closest_hospital_gdf = gpd.GeoDataFrame(closest_hospital_rows, crs='EPSG:2154')
distance_matrix_gdf = gpd.GeoDataFrame(distance_matrix_rows, crs='EPSG:2154')

# 6. Isochrones (Layer 4)
# 15 minutes = 900 seconds
print("Calculating isochrones...")
isochrones_rows = []
for idx, h_row in hospitals_nodes_gdf.iterrows():
    h_node = hospitals_nodes_idx[idx]
    
    try:
        # Find nodes reachable within 900s
        subgraph_nodes = nx.single_source_dijkstra_path_length(G_proj, h_node, cutoff=900, weight='travel_time')
        
        if not subgraph_nodes:
            continue
            
        # Extract reachable nodes' geometry
        node_geoms = []
        for n in subgraph_nodes.keys():
            node_geoms.append(Point(G_proj.nodes[n]['x'], G_proj.nodes[n]['y']))
        
        if node_geoms:
            # Create the isochrone polygon by buffering the union of reachable nodes
            # A small buffer on the union of nodes is a common approximation
            # Using a slightly larger buffer to create a cohesive area
            isochrone_poly = gpd.GeoSeries(node_geoms).unary_union.buffer(200)
            
            isochrones_rows.append({
                'hospital_name': h_row['hospital_name'],
                'travel_time_min': 15.0,
                'geometry': isochrone_poly
            })
            
    except Exception as e:
        print(f"Error calculating isochrone for {h_row['hospital_name']}: {e}")
        continue

isochrones_gdf = gpd.GeoDataFrame(isochrones_rows, crs='EPSG:2154')

# 7. Final Layers and Save
# Incidents layer
incidents_final = incidents_proj[['incident_id', 'geometry']].copy()

# Closest Hospital layer
closest_hospital_final = closest_hospital_gdf[['incident_id', 'hospital_name', 'network_distance_m', 'geometry']]

# Distance Matrix layer
distance_matrix_final = distance_matrix_gdf[['incident_id', 'hospital_name', 'rank', 'network_distance_m', 'geometry']]

# Isochrones layer
isochrones_final = isochrones_gdf[['hospital_name', 'travel_time_min', 'geometry']]

# Save to GPKG
print("Saving to emergency_routing.gpkg...")
incidents_final.to_file("/work/emergency_routing.gpkg", layer='incidents', driver='GPKG')
closest_hospital_final.to_file("/work/emergency_routing.gpkg", layer='closest_hospital', driver='GPKG', mode='a')
distance_matrix_final.to_file("/work/emergency_routing.gpkg", layer='distance_matrix', driver='GPKG', mode='a')
isochrones_final.to_file("/work/emergency_routing.gpkg", layer='isochrones_15min', driver='GPKG', mode='a')

print("Done!")