spa-l3-paris-emergency-routing

Rerun

import pandas as pd
import geopandas as gpd
from shapely.geometry import Point

# Load incidents
df = pd.read_csv('incidents.csv')

# Convert to GeoDataFrame
# Coordinates are lat, lon. GeoPandas expects lon, lat (x, y)
geometry = [Point(xy) for xy in zip(df['longitude'], df['latitude'])]
incidents_gdf = gpd.GeoDataFrame(df, geometry=geometry, crs="EPSG:4326")

print("Incidents loaded:")
print(incidents_gdf)
print(incidents_gdf.crs)

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

def solve():
    bbox = (48.83, 2.30, 48.88, 2.38) # min_lat, min_lon, max_lat, max_lon
    
    print("Downloading road network...")
    # osmnx uses (north, south, east, west) or bbox
    # Let's use the bbox as (north, south, east, west)
    # bbox = (48.88, 48.83, 2.38, 2.30)
    G = ox.graph_from_bbox(48.88, 48.83, 2.38, 2.30, network_type='drive')
    
    # Add speed and travel time to edges
    # osmnx's add_edge_speeds and add_edge_travel_times use the maxspeed tag 
    # and default to 30 km/h if missing? Actually we should check defaults.
    # The prompt says "Assume 30 km/h for edges without a maxspeed tag; use posted maxspeed otherwise."
    G = ox.add_edge_speeds(G, fallback=30)
    G = ox.add_edge_travel_times(G)

    print("Downloading hospitals...")
    # Query hospitals via Overpass (using osmnx or similar)
    # osmnx can query tags
    tags = {'amenity': 'hospital'}
    # Since we are constrained by bbox, let's fetch amenities in the area
    # We can use ox.features_from_bbox
    hospitals_gdf = ox.features_from_bbox(48.88, 48.83, 2.38, 2.30, tags=tags)
    
    # Filter for hospital points or polygons and get a name
    # Some might not have a name, let's provide a fallback
    hospitals_gdf['hospital_name'] = hospitals_gdf.get('name', 'Unnamed Hospital')
    # Ensure we have points for hospital locations (use centroid if polygon)
    hospitals_gdf['geometry'] = hospitals_gdf.geometry.centroid
    
    print(f"Found {len(hospitals_gdf)} hospitals.")

    # Load incidents
    df_incidents = pd.read_csv('incidents.csv')
    incidents_gdf = gpd.GeoDataFrame(
        df_incidents, 
        geometry=[Point(xy) for xy in zip(df_incidents['longitude'], df_incidents['latitude'])],
        crs="EPSG:4326"
    )

    # Project to Paris metric system (EPSG:2154 - Lambert-93)
    target_crs = "EPSG:2154"
    incidents_gdf = incidents_gdf.to_crs(target_crs)
    hospitals_gdf = hospitals_gdf.to_crs(target_crs)
    G_proj = ox.project_graph(G, to_crs=target_crs)

    # Pre-calculate node positions for routing
    # Note: G might be in 4326, but ox.project_graph converts it.
    # Let's make sure G and its nodes are in target_crs.
    # Wait, the prompt implies we use the network. If we use ox.project_graph, 
    # the nodes are in target_crs.
    
    # Let's re-verify G is projected
    G = ox.project_graph(G, to_crs=target_crs)

    # Prepare outputs
    closest_hospital_rows = []
    distance_matrix_rows = []

    print("Calculating routes...")
    for _, incident in incidents_gdf.iterrows():
        inc_id = incident['incident_id']
        inc_geom = incident.geometry
        
        # Find nearest node in the network
        inc_node = ox.distance.nearest_nodes(G, inc_geom.x, inc_geom.y)
        
        # Find all hospitals and calculate distances to incident
        # We'll find the distances from the incident to ALL hospitals first 
        # to get the top 3.
        
        # For each hospital, find the nearest node
        hosp_data = []
        for _, hosp in hospitals_gdf.iterrows():
            h_name = hosp['hospital_name']
            h_geom = hosp.geometry
            h_node = ox.distance.nearest_nodes(G, h_geom.x, h_geom.y)
            
            try:
                # Shortest path distance (network distance)
                dist = nx.shortest_path_length(G, inc_node, h_node, weight='length')
                hosp_data.append({
                    'hospital_name': h_name,
                    'node': h_node,
                    'network_distance_m': dist
                })
            except nx.NetworkXNoPath:
                continue
        
        # Sort hospitals by distance
        hosp_data = sorted(hosp_data, key=lambda x: x['network_distance_m'])
        
        # 1. Closest hospital route
        if hosp_data:
            closest = hosp_data[0]
            try:
                path = nx.shortest_path(G, inc_node, closest['node'], weight='length')
                # Convert path to LineString
                route_nodes = [G.nodes[n]['x'], G.nodes[n]['y']] # This is wrong, need geom
                # Better way to get geometry of path
                route_edges = list(zip(path[:-1], path[1:]))
                edge_geoms = []
                for u, v in route_edges:
                    # There might be multiple edges between u and v, take the shortest
                    edge_data = G.get_edge_data(u, v)
                    # edge_data is a dict of dicts (key is key)
                    # We want the shortest one for simplicity or just one.
                    # Usually, for 'drive', edges are unique-ish.
                    best_key = min(edge_data.keys(), key=lambda k: edge_data[k].get('length', float('inf')))
                    edge_geoms.append(edge_data[best_key]['geometry'])
                
                from shapely.ops import linemerge
                route_line = linemerge(edge_geoms)
                
                closest_hospital_rows.append({
                    'incident_id': inc_id,
                    'hospital_name': closest['hospital_name'],
                    'network_distance_m': closest['network_distance_m'],
                    'geometry': route_line
                })
            except Exception as e:
                print(f"Error calculating route for {inc_id}: {e}")

            # 2. Distance Matrix (Top 3)
            for rank, h in enumerate(hosp_data[:3], 1):
                distance_matrix_rows.append({
                    'incident_id': inc_id,
                    'hospital_name': h['hospital_name'],
                    'rank': rank,
                    'network_distance_m': h['network_distance_m']
                })

    # Create Layers
    # Layer 1: incidents
    incidents_out = incidents_gdf[['incident_id', 'geometry']]

    # Layer 2: closest_hospital
    closest_hospital_gdf = gpd.GeoDataFrame(closest_hospital_rows, crs=target_crs)
    # Note: the prompt says "LineString (or MultiLineString) geometry"
    # linemerge returns LineString or MultiLineString.

    # Layer 3: distance_matrix
    distance_matrix_gdf = gpd.GeoDataFrame(distance_matrix_rows, crs=target_crs)
    # Tabular, geometry is empty/null.
    
    # Layer 4: isochrones_15min
    print("Calculating isochrones...")
    isochrone_rows = []
    for _, hosp in hospitals_gdf.iterrows():
        h_name = hosp['hospital_name']
        h_node = ox.distance.nearest_nodes(G, hosp.geometry.x, hosp.geometry.y)
        
        # 15 minutes at 30 km/h (avg) is not quite right because we should use travel time.
        # Prompt: "15-minute drive-time isochrone around every hospital"
        # We need to find all nodes reachable within 15 minutes.
        target_time_seconds = 15 * 60
        
        # Use nx.single_source_dijkstra to find nodes within target_time_seconds
        try:
            subgraph_nodes = nx.single_source_dijkstra_path_length(G, h_node, cutoff=target_time_seconds, weight='travel_time')
            
            if subgraph_nodes:
                # Get geometries of these nodes
                node_points = [Point(G.nodes[n]['x'], G.nodes[n]['y']) for n in subgraph_nodes.keys()]
                # Create a convex hull or alpha shape? 
                # Usually, an isochrone is the union of geometries of edges/nodes.
                # A common way: buffer the nodes and union them, or find the boundary of the reachable subgraph.
                # Let's use the union of edges in the subgraph for a more accurate isochrone.
                
                # To be more robust and faster, let's find the edges in the subgraph
                subgraph_edges = []
                for u, v, k, data in G.edges(subgraph_nodes.keys(), data=True, keys=True):
                    if u in subgraph_nodes and v in subgraph_nodes:
                        # Check if travel time from h_node to u + edge time u->v <= cutoff
                        # But single_source_dijkstra already handles that.
                        # We need to verify if the edge itself is part of the "reachable" set.
                        # The easiest way is to check if both ends are in the subgraph_nodes.
                        # Wait, we need to be careful about the time to reach the middle of the edge.
                        # For simplicity, if both ends are in, the edge is reachable.
                        subgraph_edges.append(data['geometry'])
                
                if subgraph_edges:
                    isochrone_geom = MultiPolygon([Polygon(subgraph_edges[0].exterior).buffer(10)]) # Placeholder logic
                    # Let's do it better:
                    from shapely.ops import unary_union
                    # Buffer the nodes slightly to create a polygon
                    node_geoms = [Point(G.nodes[n]['x'], G.nodes[n]['y']).buffer(20) for n in subgraph_nodes.keys()]
                    isochrone_geom = unary_union(node_geoms).convex_hull # Convex hull is too simple.
                    # Actually, let's buffer the edges.
                    edge_geoms = [data['geometry'].buffer(25) for u, v, k, data in G.edges(subgraph_nodes.keys(), data=True, keys=True) 
                                  if u in subgraph_nodes and v in subgraph_nodes]
                    if edge_geoms:
                        isochrone_geom = unary_union(edge_geoms).convex_hull # Still convex hull? No.
                        # Let's use unary_union of buffered edges.
                        isochrone_geom = unary_union(edge_geoms)
                        # A proper isochrone is usually the boundary. 
                        # But the prompt asks for MultiPolygon.
                        # unary_union of buffered edges should give a good shape.
                        
                        isochrone_rows.append({
                            'hospital_name': h_name,
                            'travel_time_min': 15.0,
                            'geometry': isochrone_geom
                        })
        except Exception as e:
            print(f"Error calculating isochrone for {h_name}: {e}")

    isochrones_gdf = gpd.GeoDataFrame(isochrone_rows, crs=target_crs)

    # Save to GPKG
    print("Saving to GPKG...")
    with gpd.io.file.FileHandler("emergency_routing.gpkg") as fh:
        incidents_out.to_file("emergency_routing.gpkg", layer='incidents', driver="GPKG")
        closest_hospital_gdf.to_file("emergency_routing.gpkg", layer='closest_hospital', driver="GPKG")
        distance_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!")

if __name__ == "__main__":
    solve()

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

def solve():
    # Bounding box: min_lat, min_lon, max_lat, max_lon
    # Prompt: 48.83, 2.30 to 48.88, 2.38
    north, south, east, west = 48.88, 48.83, 2.38, 2.30
    
    print("Downloading road network...")
    G = ox.graph_from_bbox(north, south, east, west, network_type='drive')
    
    # Add speeds and travel times
    # Prompt: "Assume 30 km/h for edges without a maxspeed tag; use posted maxspeed otherwise."
    # ox.add_edge_speeds uses 45 km/h by default in some versions, so we set fallback=30.
    G = ox.add_edge_speeds(G, fallback=30)
    G = ox.add_edge_travel_times(G)

    print("Downloading hospitals...")
    tags = {'amenity': 'hospital'}
    hospitals_gdf = ox.features_from_bbox(north, south, east, west, tags=tags)
    
    # Prepare hospital data
    hospitals_gdf['hospital_name'] = hospitals_gdf.get('name', 'Unnamed Hospital')
    # Use centroid for point location
    hospitals_gdf['geometry'] = hospitals_gdf.geometry.centroid
    
    print(f"Found {len(hospitals_gdf)} hospitals.")

    # Load incidents
    df_incidents = pd.read_csv('incidents.csv')
    incidents_gdf = gpd.GeoDataFrame(
        df_incidents, 
        geometry=[Point(xy) for xy in zip(df_incidents['longitude'], df_incidents['latitude'])],
        crs="EPSG:4326"
    )

    # Project to Paris metric system (EPSG:2154 - Lambert-93)
    target_crs = "EPSG:2154"
    incidents_gdf = incidents_gdf.to_crs(target_crs)
    hospitals_gdf = hospitals_gdf.to_crs(target_crs)
    G = ox.project_graph(G, to_crs=target_crs)

    # Prepare outputs
    closest_hospital_rows = []
    distance_matrix_rows = []
    isochrone_rows = []

    print("Calculating routes and distance matrix...")
    for _, incident in incidents_gdf.iterrows():
        inc_id = incident['incident_id']
        inc_geom = incident.geometry
        
        # Find nearest node in the network
        inc_node = ox.distance.nearest_nodes(G, inc_geom.x, inc_geom.y)
        
        # Find all hospitals and calculate distances to incident
        hosp_data = []
        for _, hosp in hospitals_gdf.iterrows():
            h_name = hosp['hospital_name']
            h_geom = hosp.geometry
            h_node = ox.distance.nearest_nodes(G, h_geom.x, h_geom.y)
            
            try:
                # Shortest path distance (network distance)
                dist = nx.shortest_path_length(G, inc_node, h_node, weight='length')
                hosp_data.append({
                    'hospital_name': h_name,
                    'node': h_node,
                    'network_distance_m': dist
                })
            except nx.NetworkXNoPath:
                continue
        
        # Sort hospitals by distance
        hosp_data = sorted(hosp_data, key=lambda x: x['network_distance_m'])
        
        # 1. Closest hospital route
        if hosp_data:
            closest = hosp_data[0]
            try:
                path = nx.shortest_path(G, inc_node, closest['node'], weight='length')
                route_edges = list(zip(path[:-1], path[1:]))
                edge_geoms = []
                for u, v in route_edges:
                    edge_data = G.get_edge_data(u, v)
                    # Get the geometry of the shortest edge between u and v
                    best_key = min(edge_data.keys(), key=lambda k: edge_data[k].get('length', float('inf')))
                    edge_geoms.append(edge_data[best_key]['geometry'])
                
                route_line = linemerge(edge_geoms)
                
                closest_hospital_rows.append({
                    'incident_id': inc_id,
                    'hospital_name': closest['hospital_name'],
                    'network_distance_m': float(closest['network_distance_m']),
                    'geometry': route_line
                })
            except Exception as e:
                print(f"Error calculating route for {inc_id}: {e}")

            # 2. Distance Matrix (Top 3)
            for rank, h in enumerate(hosp_data[:3], 1):
                distance_matrix_rows.append({
                    'incident_id': inc_id,
                    'hospital_name': h['hospital_name'],
                    'rank': rank,
                    'network_distance_m': float(h['network_distance_m'])
                })

    # Create Layers
    # Layer 1: incidents
    incidents_out = incidents_gdf[['incident_id', 'geometry']]

    # Layer 2: closest_hospital
    closest_hospital_gdf = gpd.GeoDataFrame(closest_hospital_rows, crs=target_crs)

    # Layer 3: distance_matrix
    distance_matrix_gdf = gpd.GeoDataFrame(distance_matrix_rows, crs=target_crs)
    # Ensure it has the correct columns and is tabular
    distance_matrix_gdf = distance_matrix_gdf[['incident_id', 'hospital_name', 'rank', 'network_distance_m']]
    
    # Layer 4: isochrones_15min
    print("Calculating isochrones...")
    target_time_seconds = 15 * 60
    for _, hosp in hospitals_gdf.iterrows():
        h_name = hosp['hospital_name']
        h_node = ox.distance.nearest_nodes(G, hosp.geometry.x, hosp.geometry.y)
        
        try:
            # Find nodes reachable within 15 minutes
            subgraph_nodes = nx.single_source_dijkstra_path_length(G, h_node, cutoff=target_time_seconds, weight='travel_time')
            
            if subgraph_nodes:
                # Collect all edges in the reachable subgraph
                edge_geoms = []
                for u, v, k, data in G.edges(subgraph_nodes.keys(), data=True, keys=True):
                    if u in subgraph_nodes and v in subgraph_nodes:
                        edge_geoms.append(data['geometry'].buffer(25))
                
                if edge_geoms:
                    isochrone_geom = unary_union(edge_geoms)
                    # If the resulting isochrone is a single geometry, it might be a Polygon or MultiPolygon.
                    # unary_union returns the union. We ensure it is a MultiPolygon or Polygon.
                    if isochrone_geom.is_empty:
                        continue
                    
                    isochrone_rows.append({
                        'hospital_name': h_name,
                        'travel_time_min': 15.0,
                        'geometry': isochrone_geom
                    })
        except Exception as e:
            print(f"Error calculating isochrone for {h_name}: {e}")

    isochrones_gdf = gpd.GeoDataFrame(isochrone_rows, crs=target_crs)

    # Save to GPKG
    print("Saving to GPKG...")
    # Using standard to_file with driver='GPKG'
    incidents_out.to_file("emergency_routing.gpkg", layer='incidents', driver="GPKG")
    closest_hospital_gdf.to_file("emergency_routing.gpkg", layer='closest_hospital', driver="GPKG")
    distance_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!")

if __name__ == "__main__":
    solve()

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

def solve():
    # Bounding box: min_lat, min_lon, max_lat, max_lon
    # Prompt: 48.83, 2.30 to 48.88, 2.38
    # osmnx bbox: (north, south, east, west)
    north, south, east, west = 48.88, 48.83, 2.38, 2.30
    
    print("Downloading road network...")
    # Try both ways to call graph_from_bbox depending on version if needed, 
    # but standard is ox.graph_from_bbox(north, south, east, west) 
    # or ox.graph_from_bbox(bbox=(north, south, east, west))
    try:
        G = ox.graph_from_bbox(north, south, east, west, network_type='drive')
    except TypeError:
        G = ox.graph_from_bbox(north=north, south=south, east=east, west=west, network_type='drive')
    
    # Add speeds and travel times
    G = ox.add_edge_speeds(G, fallback=30)
    G = ox.add_edge_travel_times(G)

    print("Downloading hospitals...")
    tags = {'amenity': 'hospital'}
    try:
        hospitals_gdf = ox.features_from_bbox(north, south, east, west, tags=tags)
    except TypeError:
        hospitals_gdf = ox.features_from_bbox(north=north, south=south, east=east, west=west, tags=tags)
    
    # Prepare hospital data
    hospitals_gdf['hospital_name'] = hospitals_gdf.get('name', 'Unnamed Hospital')
    hospitals_gdf['geometry'] = hospitals_gdf.geometry.centroid
    
    print(f"Found {len(hospitals_gdf)} hospitals.")

    # Load incidents
    df_incidents = pd.read_csv('incidents.csv')
    incidents_gdf = gpd.GeoDataFrame(
        df_incidents, 
        geometry=[Point(xy) for xy in zip(df_incidents['longitude'], df_incidents['latitude'])],
        crs="EPSG:4326"
    )

    # Project to Paris metric system (EPSG:2154 - Lambert-93)
    target_crs = "EPSG:2154"
    incidents_gdf = incidents_gdf.to_crs(target_crs)
    hospitals_gdf = hospitals_gdf.to_crs(target_crs)
    G = ox.project_graph(G, to_crs=target_crs)

    # Prepare outputs
    closest_hospital_rows = []
    distance_matrix_rows = []
    isochrone_rows = []

    print("Calculating routes and distance matrix...")
    for _, incident in incidents_gdf.iterrows():
        inc_id = incident['incident_id']
        inc_geom = incident.geometry
        
        inc_node = ox.distance.nearest_nodes(G, inc_geom.x, inc_geom.y)
        
        hosp_data = []
        for _, hosp in hospitals_gdf.iterrows():
            h_name = hosp['hospital_name']
            h_geom = hosp.geometry
            h_node = ox.distance.nearest_nodes(G, h_geom.x, h_geom.y)
            
            try:
                dist = nx.shortest_path_length(G, inc_node, h_node, weight='length')
                hosp_data.append({
                    'hospital_name': h_name,
                    'node': h_node,
                    'network_distance_m': dist
                })
            except nx.NetworkXNoPath:
                continue
        
        hosp_data = sorted(hosp_data, key=lambda x: x['network_distance_m'])
        
        if hosp_data:
            closest = hosp_data[0]
            try:
                path = nx.shortest_path(G, inc_node, closest['node'], weight='length')
                route_edges = list(zip(path[:-1], path[1:]))
                edge_geoms = []
                for u, v in route_edges:
                    edge_data = G.get_edge_data(u, v)
                    best_key = min(edge_data.keys(), key=lambda k: edge_data[k].get('length', float('inf')))
                    edge_geoms.append(edge_data[best_key]['geometry'])
                
                route_line = linemerge(edge_geoms)
                
                closest_hospital_rows.append({
                    'incident_id': inc_id,
                    'hospital_name': closest['hospital_name'],
                    'network_distance_m': float(closest['network_distance_m']),
                    'geometry': route_line
                })
            except Exception as e:
                print(f"Error calculating route for {inc_id}: {e}")

            for rank, h in enumerate(hosp_data[:3], 1):
                distance_matrix_rows.append({
                    'incident_id': inc_id,
                    'hospital_name': h['hospital_name'],
                    'rank': rank,
                    'network_distance_m': float(h['network_distance_m'])
                })

    incidents_out = incidents_gdf[['incident_id', 'geometry']]
    closest_hospital_gdf = gpd.GeoDataFrame(closest_hospital_rows, crs=target_crs)
    distance_matrix_gdf = gpd.GeoDataFrame(distance_matrix_rows, crs=target_crs)
    distance_matrix_gdf = distance_matrix_gdf[['incident_id', 'hospital_name', 'rank', 'network_distance_m']]
    
    print("Calculating isochrones...")
    target_time_seconds = 15 * 60
    for _, hosp in hospitals_gdf.iterrows():
        h_name = hosp['hospital_name']
        h_node = ox.distance.nearest_nodes(G, hosp.geometry.x, hosp.geometry.y)
        
        try:
            subgraph_nodes = nx.single_source_dijkstra_path_length(G, h_node, cutoff=target_time_seconds, weight='travel_time')
            
            if subgraph_nodes:
                edge_geoms = []
                for u, v, k, data in G.edges(subgraph_nodes.keys(), data=True, keys=True):
                    if u in subgraph_nodes and v in subgraph_nodes:
                        edge_geoms.append(data['geometry'].buffer(25))
                
                if edge_geoms:
                    isochrone_geom = unary_union(edge_geoms)
                    if isochrone_geom.is_empty:
                        continue
                    
                    isochrone_rows.append({
                        'hospital_name': h_name,
                        'travel_time_min': 15.0,
                        'geometry': isochrone_geom
                    })
        except Exception as e:
            print(f"Error calculating isochrone for {h_name}: {e}")

    isochrones_gdf = gpd.GeoDataFrame(isochrone_rows, crs=target_crs)

    print("Saving to GPKG...")
    incidents_out.to_file("emergency_routing.gpkg", layer='incidents', driver="GPKG")
    closest_hospital_gdf.to_file("emergency_routing.gpkg", layer='closest_hospital', driver="GPKG")
    distance_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!")

if __name__ == "__main__":
    solve()

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

def solve():
    # Bounding box: min_lat, min_lon, max_lat, max_lon
    # Prompt: 48.83, 2.30 to 48.88, 2.38
    # osmnx bbox: (north, south, east, west)
    north, south, east, west = 48.88, 48.83, 2.38, 2.30
    
    print("Downloading road network...")
    # Correcting the way ox.graph_from_bbox is called.
    # In newer osmnx, it takes a single bbox tuple.
    try:
        G = ox.graph_from_bbox(north, south, east, west, network_type='drive')
    except TypeError:
        G = ox.graph_from_bbox(bbox=(north, south, east, west), network_type='drive')
    
    # Add speeds and travel times
    G = ox.add_edge_speeds(G, fallback=30)
    G = ox.add_edge_travel_times(G)

    print("Downloading hospitals...")
    tags = {'amenity': 'hospital'}
    try:
        hospitals_gdf = ox.features_from_bbox(north, south, east, west, tags=tags)
    except TypeError:
        hospitals_gdf = ox.features_from_bbox(bbox=(north, south, east, west), tags=tags)
    
    # Prepare hospital data
    hospitals_gdf['hospital_name'] = hospitals_gdf.get('name', 'Unnamed Hospital')
    hospitals_gdf['geometry'] = hospitals_gdf.geometry.centroid
    
    print(f"Found {len(hospitals_gdf)} hospitals.")

    # Load incidents
    df_incidents = pd.read_csv('incidents.csv')
    incidents_gdf = gpd.GeoDataFrame(
        df_incidents, 
        geometry=[Point(xy) for xy in zip(df_incidents['longitude'], df_incidents['latitude'])],
        crs="EPSG:4326"
    )

    # Project to Paris metric system (EPSG:2154 - Lambert-93)
    target_crs = "EPSG:2154"
    incidents_gdf = incidents_gdf.to_crs(target_crs)
    hospitals_gdf = hospitals_gdf.to_crs(target_crs)
    G = ox.project_graph(G, to_crs=target_crs)

    # Prepare outputs
    closest_hospital_rows = []
    distance_matrix_rows = []
    isochrone_rows = []

    print("Calculating routes and distance matrix...")
    for _, incident in incidents_gdf.iterrows():
        inc_id = incident['incident_id']
        inc_geom = incident.geometry
        
        inc_node = ox.distance.nearest_nodes(G, inc_geom.x, inc_geom.y)
        
        hosp_data = []
        for _, hosp in hospitals_gdf.iterrows():
            h_name = hosp['hospital_name']
            h_geom = hosp.geometry
            h_node = ox.distance.nearest_nodes(G, h_geom.x, h_geom.y)
            
            try:
                dist = nx.shortest_path_length(G, inc_node, h_node, weight='length')
                hosp_data.append({
                    'hospital_name': h_name,
                    'node': h_node,
                    'network_distance_m': dist
                })
            except nx.NetworkXNoPath:
                continue
        
        hosp_data = sorted(hosp_data, key=lambda x: x['network_distance_m'])
        
        if hosp_data:
            closest = hosp_data[0]
            try:
                path = nx.shortest_path(G, inc_node, closest['node'], weight='length')
                route_edges = list(zip(path[:-1], path[1:]))
                edge_geoms = []
                for u, v in route_edges:
                    edge_data = G.get_edge_data(u, v)
                    best_key = min(edge_data.keys(), key=lambda k: edge_data[k].get('length', float('inf')))
                    edge_geoms.append(edge_data[best_key]['geometry'])
                
                route_line = linemerge(edge_geoms)
                
                closest_hospital_rows.append({
                    'incident_id': inc_id,
                    'hospital_name': closest['hospital_name'],
                    'network_distance_m': float(closest['network_distance_m']),
                    'geometry': route_line
                })
            except Exception as e:
                print(f"Error calculating route for {inc_id}: {e}")

            for rank, h in enumerate(hosp_data[:3], 1):
                distance_matrix_rows.append({
                    'incident_id': inc_id,
                    'hospital_name': h['hospital_name'],
                    'rank': rank,
                    'network_distance_m': float(h['network_distance_m'])
                })

    incidents_out = incidents_gdf[['incident_id', 'geometry']]
    closest_hospital_gdf = gpd.GeoDataFrame(closest_hospital_rows, crs=target_crs)
    distance_matrix_gdf = gpd.GeoDataFrame(distance_matrix_rows, crs=target_crs)
    distance_matrix_gdf = distance_matrix_gdf[['incident_id', 'hospital_name', 'rank', 'network_distance_m']]
    
    print("Calculating isochrones...")
    target_time_seconds = 15 * 60
    for _, hosp in hospitals_gdf.iterrows():
        h_name = hosp['hospital_name']
        h_node = ox.distance.nearest_nodes(G, hosp.geometry.x, hosp.geometry.y)
        
        try:
            subgraph_nodes = nx.single_source_dijkstra_path_length(G, h_node, cutoff=target_time_seconds, weight='travel_time')
            
            if subgraph_nodes:
                edge_geoms = []
                for u, v, k, data in G.edges(subgraph_nodes.keys(), data=True, keys=True):
                    if u in subgraph_nodes and v in subgraph_nodes:
                        edge_geoms.append(data['geometry'].buffer(25))
                
                if edge_geoms:
                    isochrone_geom = unary_union(edge_geoms)
                    if isochrone_geom.is_empty:
                        continue
                    
                    isochrone_rows.append({
                        'hospital_name': h_name,
                        'travel_time_min': 15.0,
                        'geometry': isochrone_geom
                    })
        except Exception as e:
            print(f"Error calculating isochrone for {h_name}: {e}")

    isochrones_gdf = gpd.GeoDataFrame(isochrone_rows, crs=target_crs)

    print("Saving to GPKG...")
    incidents_out.to_file("emergency_routing.gpkg", layer='incidents', driver="GPKG")
    closest_hospital_gdf.to_file("emergency_routing.gpkg", layer='closest_hospital', driver="GPKG")
    distance_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!")

if __name__ == "__main__":
    solve()

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

# Custom Overpass settings for resilience
# We'll try to set the base URL for OSMnx if possible, 
# but OSMnx doesn't make it trivial to change just the base URL.
# Instead, we'll use a retry loop around the osmnx calls.

def robust_osmnx_call(func, *args, **kwargs):
    """Retry wrapper for osmnx calls."""
    endpoints = [
        None, # default
        "https://overpass.kumi.systems/api/interpreter",
        "https://lz4.overpass-api.de/api/interpreter"
    ]
    
    last_err = None
    for i, endpoint in enumerate(endpoints):
        if endpoint:
            # Note: osmnx doesn't easily allow setting a custom endpoint 
            # through a simple parameter in graph_from_bbox.
            # However, it uses requests. We can try to monkeypatch or 
            # rely on the default and hope the retry works.
            # Since we can't easily change the URL in osmnx without more complex hacks,
            # we'll just retry the default first.
            pass
        
        try:
            print(f"Attempt {i+1}...")
            return func(*args, **kwargs)
        except Exception as e:
            print(f"Attempt {i+1} failed: {e}")
            last_err = e
            if i < len(endpoints) - 1:
                time.sleep(5)
            else:
                raise last_err

def solve():
    # Bounding box: north, south, east, west
    # Prompt: 48.83, 2.30 to 48.88, 2.38
    north, south, east, west = 48.88, 48.83, 2.38, 2.30
    bbox_tuple = (north, south, east, west)
    
    print("Downloading road network...")
    # Try both ways to call graph_from_bbox
    try:
        G = robust_osmnx_call(ox.graph_from_bbox, bbox=bbox_tuple, network_type='drive')
    except Exception:
        G = robust_osmnx_call(ox.graph_from_bbox, north=north, south=south, east=east, west=west, network_type='drive')
    
    # Add speeds and travel times
    G = ox.add_edge_speeds(G, fallback=30)
    G = ox.add_edge_travel_times(G)

    print("Downloading hospitals...")
    tags = {'amenity': 'hospital'}
    try:
        hospitals_gdf = robust_osmnx_call(ox.features_from_bbox, bbox=bbox_tuple, tags=tags)
    except Exception:
        hospitals_gdf = robust_osmnx_call(ox.features_from_bbox, north=north, south=south, east=east, west=west, tags=tags)
    
    # Prepare hospital data
    hospitals_gdf['hospital_name'] = hospitals_gdf.get('name', 'Unnamed Hospital')
    hospitals_gdf['geometry'] = hospitals_gdf.geometry.centroid
    
    print(f"Found {len(hospitals_gdf)} hospitals.")

    # Load incidents
    df_incidents = pd.read_csv('incidents.csv')
    incidents_gdf = gpd.GeoDataFrame(
        df_incidents, 
        geometry=[Point(xy) for xy in zip(df_incidents['longitude'], df_incidents['latitude'])],
        crs="EPSG:4326"
    )

    # Project to Paris metric system (EPSG:2154 - Lambert-93)
    target_crs = "EPSG:2154"
    incidents_gdf = incidents_gdf.to_crs(target_crs)
    hospitals_gdf = hospitals_gdf.to_crs(target_crs)
    G = ox.project_graph(G, to_crs=target_crs)

    # Prepare outputs
    closest_hospital_rows = []
    distance_matrix_rows = []
    isochrone_rows = []

    print("Calculating routes and distance matrix...")
    for _, incident in incidents_gdf.iterrows():
        inc_id = incident['incident_id']
        inc_geom = incident.geometry
        
        inc_node = ox.distance.nearest_nodes(G, inc_geom.x, inc_geom.y)
        
        hosp_data = []
        for _, hosp in hospitals_gdf.iterrows():
            h_name = hosp['hospital_name']
            h_geom = hosp.geometry
            h_node = ox.distance.nearest_nodes(G, h_geom.x, h_geom.y)
            
            try:
                dist = nx.shortest_path_length(G, inc_node, h_node, weight='length')
                hosp_data.append({
                    'hospital_name': h_name,
                    'node': h_node,
                    'network_distance_m': dist
                })
            except nx.NetworkXNoPath:
                continue
        
        hosp_data = sorted(hosp_data, key=lambda x: x['network_distance_m'])
        
        if hosp_data:
            closest = hosp_data[0]
            try:
                path = nx.shortest_path(G, inc_node, closest['node'], weight='length')
                route_edges = list(zip(path[:-1], path[1:]))
                edge_geoms = []
                for u, v in route_edges:
                    edge_data = G.get_edge_data(u, v)
                    best_key = min(edge_data.keys(), key=lambda k: edge_data[k].get('length', float('inf')))
                    edge_geoms.append(edge_data[best_key]['geometry'])
                
                route_line = linemerge(edge_geoms)
                
                closest_hospital_rows.append({
                    'incident_id': inc_id,
                    'hospital_name': closest['hospital_name'],
                    'network_distance_m': float(closest['network_distance_m']),
                    'geometry': route_line
                })
            except Exception as e:
                print(f"Error calculating route for {inc_id}: {e}")

            for rank, h in enumerate(hosp_data[:3], 1):
                distance_matrix_rows.append({
                    'incident_id': inc_id,
                    'hospital_name': h['hospital_name'],
                    'rank': rank,
                    'network_distance_m': float(h['network_distance_m'])
                })

    incidents_out = incidents_gdf[['incident_id', 'geometry']]
    closest_hospital_gdf = gpd.GeoDataFrame(closest_hospital_rows, crs=target_crs)
    distance_matrix_gdf = gpd.GeoDataFrame(distance_matrix_rows, crs=target_crs)
    distance_matrix_gdf = distance_matrix_gdf[['incident_id', 'hospital_name', 'rank', 'network_distance_m']]
    
    print("Calculating isochrones...")
    target_time_seconds = 15 * 60
    for _, hosp in hospitals_gdf.iterrows():
        h_name = hosp['hospital_name']
        h_node = ox.distance.nearest_nodes(G, hosp.geometry.x, hosp.geometry.y)
        
        try:
            subgraph_nodes = nx.single_source_dijkstra_path_length(G, h_node, cutoff=target_time_seconds, weight='travel_time')
            
            if subgraph_nodes:
                edge_geoms = []
                for u, v, k, data in G.edges(subgraph_nodes.keys(), data=True, keys=True):
                    if u in subgraph_nodes and v in subgraph_nodes:
                        edge_geoms.append(data['geometry'].buffer(25))
                
                if edge_geoms:
                    isochrone_geom = unary_union(edge_geoms)
                    if isochrone_geom.is_empty:
                        continue
                    
                    isochrone_rows.append({
                        'hospital_name': h_name,
                        'travel_time_min': 15.0,
                        'geometry': isochrone_geom
                    })
        except Exception as e:
            print(f"Error calculating isochrone for {h_name}: {e}")

    isochrones_gdf = gpd.GeoDataFrame(isochrone_rows, crs=target_crs)

    print("Saving to GPKG...")
    incidents_out.to_file("emergency_routing.gpkg", layer='incidents', driver="GPKG")
    closest_hospital_gdf.to_file("emergency_routing.gpkg", layer='closest_hospital', driver="GPKG")
    distance_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!")

if __name__ == "__main__":
    solve()

import osmnx as ox
import networkx as nx
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
from shapely.ops import linemerge, unary_union

def solve():
    # Bounding box: north, south, east, west
    north, south, east, west = 48.88, 48.83, 2.38, 2.30
    
    print("Downloading road network and hospitals via OSMnx...")
    # We'll use the simplified approach first to minimize complexity and connection issues.
    # If these fail, the environment might have connectivity issues with OSM.
    G = ox.graph_from_bbox(north, south, east, west, network_type='drive')
    G = ox.add_edge_speeds(G, fallback=30)
    G = ox.add_edge_travel_times(G)
    
    tags = {'amenity': 'hospital'}
    hospitals_gdf = ox.features_from_bbox(north, south, east, west, tags=tags)
    hospitals_gdf['hospital_name'] = hospitals_gdf.get('name', 'Unnamed Hospital')
    hospitals_gdf['geometry'] = hospitals_gdf.geometry.centroid
    
    # Load incidents
    df_incidents = pd.read_csv('incidents.csv')
    incidents_gdf = gpd.GeoDataFrame(
        df_incidents, 
        geometry=[Point(xy) for xy in zip(df_incidents['longitude'], df_incidents['latitude'])],
        crs="EPSG:4326"
    )

    # Project to Paris metric system (EPSG:2154 - Lambert-93)
    target_crs = "EPSG:2154"
    incidents_gdf = incidents_gdf.to_crs(target_crs)
    hospitals_gdf = hospitals_gdf.to_crs(target_crs)
    G = ox.project_graph(G, to_crs=target_crs)

    # Outputs
    closest_hospital_rows = []
    distance_matrix_rows = []
    isochrone_rows = []

    print("Calculating routes and distance matrix...")
    for _, incident in incidents_gdf.iterrows():
        inc_id = incident['incident_id']
        inc_geom = incident.geometry
        inc_node = ox.distance.nearest_nodes(G, inc_geom.x, inc_geom.y)
        
        hosp_data = []
        for _, hosp in hospitals_gdf.iterrows():
            h_name = hosp['hospital_name']
            h_node = ox.distance.nearest_nodes(G, hosp.geometry.x, hosp.geometry.y)
            try:
                dist = nx.shortest_path_length(G, inc_node, h_node, weight='length')
                hosp_data.append({'hospital_name': h_name, 'node': h_node, 'network_distance_m': dist})
            except nx.NetworkXNoPath:
                continue
        
        hosp_data = sorted(hosp_data, key=lambda x: x['network_distance_m'])
        
        if hosp_data:
            # Closest hospital route
            closest = hosp_data[0]
            try:
                path = nx.shortest_path(G, inc_node, closest['node'], weight='length')
                route_edges = list(zip(path[:-1], path[1:]))
                edge_geoms = []
                for u, v in route_edges:
                    edge_data = G.get_edge_data(u, v)
                    best_key = min(edge_data.keys(), key=lambda k: edge_data[k].get('length', float('inf')))
                    edge_geoms.append(edge_data[best_key]['geometry'])
                route_line = linemerge(edge_geoms)
                closest_hospital_rows.append({
                    'incident_id': inc_id,
                    'hospital_name': closest['hospital_name'],
                    'network_distance_m': float(closest['network_distance_m']),
                    'geometry': route_line
                })
            except:
                pass

            # Distance matrix
            for rank, h in enumerate(hosp_data[:3], 1):
                distance_matrix_rows.append({
                    'incident_id': inc_id,
                    'hospital_name': h['hospital_name'],
                    'rank': rank,
                    'network_distance_m': float(h['network_distance_m'])
                })

    print("Calculating isochrones...")
    target_time_seconds = 15 * 60
    for _, hosp in hospitals_gdf.iterrows():
        h_name = hosp['hospital_name']
        h_node = ox.distance.nearest_nodes(G, hosp.geometry.x, hosp.geometry.y)
        try:
            subgraph_nodes = nx.single_source_dijkstra_path_length(G, h_node, cutoff=target_time_seconds, weight='travel_time')
            if subgraph_nodes:
                edge_geoms = [data['geometry'].buffer(25) for u, v, k, data in G.edges(subgraph_nodes.keys(), data=True, keys=True) 
                              if u in subgraph_nodes and v in subgraph_nodes]
                if edge_geoms:
                    isochrone_geom = unary_union(edge_geoms)
                    if not isochrone_geom.is_empty:
                        isochrone_rows.append({'hospital_name': h_name, 'travel_time_min': 15.0, 'geometry': isochrone_geom})
        except:
            pass

    # Save results
    incidents_out = incidents_gdf[['incident_id', 'geometry']]
    closest_hospital_gdf = gpd.GeoDataFrame(closest_hospital_rows, crs=target_crs)
    distance_matrix_gdf = gpd.GeoDataFrame(distance_matrix_rows, crs=target_crs)[['incident_id', 'hospital_name', 'rank', 'network_distance_m']]
    isochrones_gdf = gpd.GeoDataFrame(isochrone_rows, crs=target_crs)

    print("Writing to emergency_routing.gpkg...")
    incidents_out.to_file("emergency_routing.gpkg", layer='incidents', driver="GPKG")
    closest_hospital_gdf.to_file("emergency_routing.gpkg", layer='closest_hospital', driver="GPKG")
    distance_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!")

if __name__ == "__main__":
    solve()

import osmnx as ox
import networkx as nx
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
from shapely.ops import linemerge, unary_union

def solve():
    # Bounding box: north, south, east, west
    # Prompt: 48.83, 2.30 to 48.88, 2.38
    # osmnx bbox: (north, south, east, west)
    north, south, east, west = 48.88, 48.83, 2.38, 2.30
    
    print("Downloading road network and hospitals via OSMnx...")
    # Using the bbox as a tuple as requested by some osmnx versions
    bbox = (north, south, east, west)
    
    try:
        G = ox.graph_from_bbox(bbox, network_type='drive')
    except TypeError:
        # If that fails, try passing them individually as keyword arguments
        # Note: The error said it takes 1 positional argument (the bbox)
        # or it might take keyword arguments.
        G = ox.graph_from_bbox(north=north, south=south, east=east, west=west, network_type='drive')
    
    G = ox.add_edge_speeds(G, fallback=30)
    G = ox.add_edge_travel_times(G)
    
    tags = {'amenity': 'hospital'}
    try:
        hospitals_gdf = ox.features_from_bbox(bbox, tags=tags)
    except TypeError:
        hospitals_gdf = ox.features_from_bbox(north=north, south=south, east=east, west=west, tags=tags)
    
    hospitals_gdf['hospital_name'] = hospitals_gdf.get('name', 'Unnamed Hospital')
    hospitals_gdf['geometry'] = hospitals_gdf.geometry.centroid
    
    print(f"Found {len(hospitals_gdf)} hospitals.")

    df_incidents = pd.read_csv('incidents.csv')
    incidents_gdf = gpd.GeoDataFrame(
        df_incidents, 
        geometry=[Point(xy) for xy in zip(df_incidents['longitude'], df_incidents['latitude'])],
        crs="EPSG:4326"
    )

    target_crs = "EPSG:2154"
    incidents_gdf = incidents_gdf.to_crs(target_crs)
    hospitals_gdf = hospitals_gdf.to_crs(target_crs)
    G = ox.project_graph(G, to_crs=target_crs)

    closest_hospital_rows = []
    distance_matrix_rows = []
    isochrone_rows = []

    print("Calculating routes and distance matrix...")
    for _, incident in incidents_gdf.iterrows():
        inc_id = incident['incident_id']
        inc_geom = incident.geometry
        inc_node = ox.distance.nearest_nodes(G, inc_geom.x, inc_geom.y)
        
        hosp_data = []
        for _, hosp in hospitals_gdf.iterrows():
            h_name = hosp['hospital_name']
            h_node = ox.distance.nearest_nodes(G, hosp.geometry.x, hosp.geometry.y)
            try:
                dist = nx.shortest_path_length(G, inc_node, h_node, weight='length')
                hosp_data.append({'hospital_name': h_name, 'node': h_node, 'network_distance_m': dist})
            except nx.NetworkXNoPath:
                continue
        
        hosp_data = sorted(hosp_data, key=lambda x: x['network_distance_m'])
        
        if hosp_data:
            closest = hosp_data[0]
            try:
                path = nx.shortest_path(G, inc_node, closest['node'], weight='length')
                route_edges = list(zip(path[:-1], path[1:]))
                edge_geoms = []
                for u, v in route_edges:
                    edge_data = G.get_edge_data(u, v)
                    best_key = min(edge_data.keys(), key=lambda k: edge_data[k].get('length', float('inf')))
                    edge_geoms.append(edge_data[best_key]['geometry'])
                route_line = linemerge(edge_geoms)
                closest_hospital_rows.append({
                    'incident_id': inc_id,
                    'hospital_name': closest['hospital_name'],
                    'network_distance_m': float(closest['network_distance_m']),
                    'geometry': route_line
                })
            except:
                pass

            for rank, h in enumerate(hosp_data[:3], 1):
                distance_matrix_rows.append({
                    'incident_id': inc_id,
                    'hospital_name': h['hospital_name'],
                    'rank': rank,
                    'network_distance_m': float(h['network_distance_m'])
                })

    print("Calculating isochrones...")
    target_time_seconds = 15 * 60
    for _, hosp in hospitals_gdf.iterrows():
        h_name = hosp['hospital_name']
        h_node = ox.distance.nearest_nodes(G, hosp.geometry.x, hosp.geometry.y)
        try:
            subgraph_nodes = nx.single_source_dijkstra_path_length(G, h_node, cutoff=target_time_seconds, weight='travel_time')
            if subgraph_nodes:
                edge_geoms = []
                for u, v, k, data in G.edges(subgraph_nodes.keys(), data=True, keys=True):
                    if u in subgraph_nodes and v in subgraph_nodes:
                        edge_geoms.append(data['geometry'].buffer(25))
                if edge_geoms:
                    isochrone_geom = unary_union(edge_geoms)
                    if not isochrone_geom.is_empty:
                        isochrone_rows.append({'hospital_name': h_name, 'travel_time_min': 15.0, 'geometry': isochrone_geom})
        except:
            pass

    incidents_out = incidents_gdf[['incident_id', 'geometry']]
    closest_hospital_gdf = gpd.GeoDataFrame(closest_hospital_rows, crs=target_crs)
    distance_matrix_gdf = gpd.GeoDataFrame(distance_matrix_rows, crs=target_crs)[['incident_id', 'hospital_name', 'rank', 'network_distance_m']]
    isochrones_gdf = gpd.GeoDataFrame(isochrone_rows, crs=target_crs)

    print("Saving to GPKG...")
    incidents_out.to_file("emergency_routing.gpkg", layer='incidents', driver="GPKG")
    closest_hospital_gdf.to_file("emergency_routing.gpkg", layer='closest_hospital', driver="GPKG")
    distance_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!")

if __name__ == "__main__":
    solve()