spa-l3-paris-emergency-routing

Rerun

Let me break down what's needed:

1. Read incidents from `incidents.csv` - 8 emergency call locations in Paris (bbox 48.83,2.30 to 48.88,2.38)
2. Overpass query for driveable road network and hospitals
3. Compute shortest driving routes from each incident to closest hospital
4. Distance matrix: each incident vs its 3 nearest hospitals
5. 15-minute drive-time isochrones around each hospital
6. Output GPKG with 4 layers: `incidents`, `closest_hospital`, `distance_matrix`, `isochrones_15min`

Assumptions:
- 30 km/h for edges without maxspeed
- Use posted maxspeed otherwise
- CRS: official metric coordinate system for Paris = EPSG:2154 (RGF93 / Lambert-93) is the standard for France. Actually, for Paris specifically, EPSG:2154 is the national metric system. Let me use that.

Let me start by reading the incidents file and understanding the structure.

Let me plan the solution:

1. Read the incidents CSV
2. Query Overpass for road network and hospitals in the bbox
3. Build a routable graph, compute shortest paths
4. Generate the three deliverables:
   - closest_hospital: one route per incident to its closest hospital
   - distance_matrix: each incident vs its 3 nearest hospitals
   - isochrones_15min: 15-min drive-time isochrone around each hospital

Let me use osmnx for Overpass queries and network analysis, as it's pre-installed and wraps OSM network analysis nicely.

Actually, let me think about this carefully. The bbox is 48.83,2.30 to 48.88,2.38 (min_lat, min_lon to max_lat, max_lon). This is in Paris.

I'll use osmnx to:
1. Get the drivable road network
2. Get hospitals (amenity=hospital nodes/ways)
3. Build a graph and compute shortest paths
4. Generate isochrones

Let me write a comprehensive script.

For the isochrones, I can use osmnx's `ox.isochrones_from_graph` or manually compute using networkx. Actually, osmnx has `ox.isochrones_from_nodes` but that's deprecated. Let me use `ox.isochrones_from_graph`.

Wait, let me think about the right approach:

1. Download OSM data for the area
2. Build a graph with speeds (30 km/h default, posted maxspeed otherwise)
3. Compute edge travel times
4. For each incident, find nearest hospitals (by network distance)
5. Compute shortest paths
6. Generate isochrones

Let me use osmnx for graph building and routing.

For hospitals: I'll query Overpass for nodes/ways with amenity=hospital and get their centroids.

Actually, `osmnx` can do all of this. Let me use `ox.graph_from_bbox` for the road network, and `ox.features_from_bbox` for hospitals.

Then for the graph, I need to add speeds (using the maxspeed tag or default 30 km/h) and compute travel times.

For isochrones, I'll use `ox.isochrones_from_graph` which returns polygons for a given travel time.

Let me write the script:

```python
import osmnx as ox
import networkx as nx
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point, LineString, MultiPoint, MultiPolygon
from itertools import combinations
import numpy as np

# Read incidents
incidents_df = pd.read_csv('/work/incidents.csv')

# Bbox
north, south = 48.88, 48.83
east, west = 2.38, 2.30

# Download road network
G = ox.graph_from_bbox(north, south, east, west, network_type='drive', simplify=True)

# Add speeds and travel times
# For edges without maxspeed, default to 30 km/h
for u, v, k, data in G.edges(data=True, keys=True):
    if 'maxspeed' in data:
        # maxspeed could be a list or string
        if isinstance(data['maxspeed'], list):
            maxspeed = float(data['maxspeed'][0])
        else:
            maxspeed = float(data['maxspeed'])
        data['speed_kph'] = maxspeed
    else:
        data['speed_kph'] = 30.0
    # travel time in minutes
    data['travel_time_min'] = (data['length'] / 1000) / data['speed_kph'] * 60

# Download hospitals
hospitals = ox.features_from_bbox(north, south, east, west, tags={'amenity': 'hospital'})

# Process hospitals - get point geometries (centroids)
hospital_points = []
for idx, row in hospitals.iterrows():
    geom = row.geometry
    if geom.geom_type in ('Point', 'MultiPoint'):
        if geom.geom_type == 'MultiPoint':
            centroid = geom.centroid
        else:
            centroid = geom
    else:
        centroid = geom.centroid
    name = row.get('name', f"Hospital_{idx}")
    hospital_points.append({'hospital_name': name, 'geometry': centroid})

hospitals_gdf = gpd.GeoDataFrame(hospital_points, crs='EPSG:4326')

# Get the nearest node on the graph for each point
incident_nodes = []
for idx, row in incidents_df.iterrows():
    node = ox.distance.nearest_nodes(G, row['longitude'], row['latitude'])
    incident_nodes.append(node)

hospital_nodes = []
for idx, row in hospitals_gdf.iterrows():
    node = ox.distance.nearest_nodes(G, row.geometry.x, row.geometry.y)
    hospital_nodes.append(node)

# For each incident, compute shortest path (by travel time) to all hospitals
# Then find closest hospital and top 3

results_closest = []
results_matrix = []

for i, (inc_id, inc_node) in enumerate(zip(incidents_df['incident_id'], incident_nodes)):
    distances = []
    for j, (hosp_name, hosp_node) in enumerate(zip(hospitals_gdf['hospital_name'], hospital_nodes)):
        try:
            # Use travel time as weight
            length = nx.shortest_path_length(G, inc_node, hosp_node, weight='travel_time_min')
            path = nx.shortest_path(G, inc_node, hosp_node, weight='travel_time_min')
            # Get the geometric path
            path_edges = list(zip(path[:-1], path[1:]))
            # Build the LineString
            coords = []
            for u, v in path_edges:
                # Get edge data
                edge_data = G[u][v]
                # Pick the first key
                key = list(edge_data.keys())[0]
                geom = edge_data[key].get('geometry', None)
                if geom:
                    coords.extend(list(geom.coords))
                else:
                    coords.append((G.nodes[u]['x'], G.nodes[u]['y']))
                    coords.append((G.nodes[v]['x'], G.nodes[v]['y']))
            # Simplify by taking unique coords
            unique_coords = [coords[0]]
            for c in coords[1:]:
                if c != unique_coords[-1]:
                    unique_coords.append(c)
            route_line = LineString(unique_coords)
            
            # Network distance in meters
            length_m = nx.shortest_path_length(G, inc_node, hosp_node, weight='length')
            
            distances.append({
                'incident_id': inc_id,
                'hospital_name': hosp_name,
                'network_distance_m': length_m,
                'travel_time_min': length,
                'route_geom': route_line
            })
        except nx.NetworkXNoPath:
            pass
    
    # Sort by travel time
    distances.sort(key=lambda x: x['travel_time_min'])
    
    # Closest hospital
    if distances:
        closest = distances[0]
        results_closest.append(closest)
        
        # Top 3 for matrix
        for rank, d in enumerate(distances[:3], 1):
            results_matrix.append({
                'incident_id': d['incident_id'],
                'hospital_name': d['hospital_name'],
                'rank': rank,
                'network_distance_m': d['network_distance_m']
            })

# Build closest_hospital gdf
closest_gdf = gpd.GeoDataFrame(results_closest, crs='EPSG:4326', geometry='route_geom')
closest_gdf = closest_gdf[['incident_id', 'hospital_name', 'network_distance_m', 'route_geom']]
closest_gdf = closest_gdf.rename(columns={'route_geom': 'geometry'})
closest_gdf.set_geometry('geometry', inplace=True)

# Build distance_matrix gdf
matrix_gdf = gpd.GeoDataFrame(results_matrix, crs='EPSG:4326')
matrix_gdf['geometry'] = None

# Now isochrones - 15 minute drive time around each hospital
# Use osmnx isochrones
isochrone_results = []
for j, (hosp_name, hosp_node) in enumerate(zip(hospitals_gdf['hospital_name'], hospital_nodes)):
    try:
        iso_poly = ox.isochrones_from_graph(G, [hosp_node], edge_buff=25, node_buff=0, 
                                             n_buff=50, trip_times=[15], weight='travel_time_min')
        if len(iso_poly) > 0:
            isochrone_results.append({
                'hospital_name': hosp_name,
                'travel_time_min': 15,
                'geometry': iso_poly[0]
            })
    except Exception as e:
        print(f"Error for {hosp_name}: {e}")

isochrone_gdf = gpd.GeoDataFrame(isochrone_results, crs='EPSG:4326')

# Project to Paris metric CRS (EPSG:2154 - RGF93 / Lambert-93 - official for France)
target_crs = 'EPSG:2154'

incidents_gdf = gpd.GeoDataFrame(
    incidents_df,
    geometry=[Point(xy) for xy in zip(incidents_df['longitude'], incidents_df['latitude'])],
    crs='EPSG:4326'
)

closest_gdf = closest_gdf.to_crs(target_crs)
matrix_gdf = matrix_gdf.to_crs(target_crs)
isochrone_gdf = isochrone_gdf.to_crs(target_crs)
incidents_gdf = incidents_gdf.to_crs(target_crs)

# Write to GPKG
output_path = '/work/emergency_routing.gpkg'

# Layer: incidents
incidents_gdf.to_file(output_path, layer='incidents', driver='GPKG')

# Layer: closest_hospital
closest_gdf.to_file(output_path, layer='closest_hospital', driver='GPKG')

# Layer: distance_matrix
matrix_gdf.to_file(output_path, layer='distance_matrix', driver='GPKG')

# Layer: isochrones_15min
isochrone_gdf.to_file(output_path, layer='isochrones_15min', driver='GPKG')

print("Done!")
```

Wait, I need to be careful about a few things:

1. `ox.isochrones_from_graph` - let me check the API. In recent osmnx versions, it might be `ox.isochrones()` or something else. Let me think...

In modern osmnx (1.x+), the function is `ox.isochrones_from_graph` which takes a graph, a list of nodes, and trip times. It returns a list of shapely Polygons (one per trip time, or multiple per time if separate is True).

Actually, I think it might return a GeoDataFrame. Let me use a safer approach.

Actually, looking at osmnx docs more carefully, the function signature is:

```python
ox.isochrones_from_graph(G, nodes, trip_times=10, weight='length', ...)
```

It returns a GeoDataFrame with isochrone polygons. Let me use this.

But wait - for isochrones, I need to specify from which nodes to compute. For each hospital, I want a 15-min drive time area. I should use the hospital node.

2. For the route geometry - I need to be more careful about how I extract the linestring from the OSMnx graph. The `nodes` in the graph have x,y attributes, and edges might have geometry (for curved roads). Let me use `ox.utils_graph.route_to_gdf` or just extract coords properly.

Actually, a simpler approach: use `ox.shortest_path` to get the node path, then `ox.utils_graph.get_route_edge_attributes` or just create the geometry from node coordinates.

Let me use `nx.shortest_path` to get node list, then get coordinates from the nodes.

Actually, the most reliable way to get a route LineString is:

```python
route_nodes = nx.shortest_path(G, orig_node, dest_node, weight='travel_time_min')
route_gdf = ox.utils_graph.route_to_gdf(G, route_nodes, weight='travel_time_min')
```

Wait, `route_to_gdf` might have been renamed. Let me check...

In osmnx 1.x, it's `ox.utils_graph.route_to_gdf` which takes a graph and node list and returns a GeoDataFrame with the route geometry.

Actually, I think a simpler approach is:

```python
route_nodes = nx.shortest_path(G, orig, dest, weight='travel_time_min')
route_edges = list(zip(route_nodes[:-1], route_nodes[1:]))
# Build linestring
points = [(G.nodes[n]['x'], G.nodes[n]['y']) for n in route_nodes]
route_line = LineString(points)
```

This gives a simplified route but it's accurate enough for routing purposes. However, for curved roads, edges may have `geometry` attribute with more detailed shape points. Let me use the edge geometries when available.

Let me write the script more carefully.

3. CRS: The official metric coordinate system for Paris (France) is EPSG:2154 (RGF93 / Lambert-93) or EPSG:27561 (NTF Paris / Lambert zone I) or EPSG:27571 (NTF Paris / Lambert zone II). For modern France, EPSG:2154 is the official one. Let me use that.

Actually, for Paris specifically in the context of French official mapping, EPSG:2154 (RGF93 / Lambert-93) is the official national system. Let me use that.

Let me write the full script now.

One concern: the OSMnx query might take some time. Let me set a timeout and handle errors.

Also, I need to be careful about the `ox.settings` for the timeout.

Let me write it all out.

Actually, let me also think about edge cases:
- What if no hospitals are found? (unlikely in Paris)
- What if a shortest path doesn't exist?
- What if isochrone generation fails for a node?

Let me handle these gracefully.

Let me write the script now.

```python
import osmnx as ox
import networkx as nx
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point, LineString, MultiPoint, box
import numpy as np

# Configure osmnx
ox.config(use_cache=True, log_console=True, timeout=180)

# Read incidents
incidents_df = pd.read_csv('/work/incidents.csv')

# Bbox: 48.83,2.30 to 48.88,2.38
north, south = 48.88, 48.83
east, west = 2.38, 2.30

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

print(f"Graph has {len(G.nodes)} nodes and {len(G.edges)} edges")

# Add speeds and travel times
for u, v, k, data in G.edges(data=True, keys=True):
    if 'maxspeed' in data and data['maxspeed']:
        ms = data['maxspeed']
        if isinstance(ms, list):
            ms = ms[0]
        try:
            speed = float(ms.replace(' km/h', '').replace(' mph', '').strip())
        except (ValueError, AttributeError):
            speed = 30.0
        data['speed_kph'] = speed
    else:
        data['speed_kph'] = 30.0
    # travel time in minutes: length(m) / (speed_kph * 1000/3600) / 60 = length / (speed_kph / 3.6) / 60 = length / speed_kph * 3.6 / 60
    # Actually: time_min = distance_m / (speed_kph * 1000 / 3600) / 60
    # = distance_m / (speed_kph / 3.6) / 60
    # = distance_m * 3.6 / speed_kph / 60
    # = distance_m * 0.06 / speed_kph
    data['travel_time_min'] = data['length'] * 0.06 / data['speed_kph']

print("Downloading hospitals...")
# Download hospitals
hospitals_gdf = ox.features_from_bbox(north, south, east, west, tags={'amenity': 'hospital'})
print(f"Found {len(hospitals_gdf)} hospital features")

# Process hospitals - get point geometries
hospital_records = []
for idx, row in hospitals_gdf.iterrows():
    geom = row.geometry
    if geom is None:
        continue
    if geom.geom_type == 'Polygon':
        centroid = geom.centroid
    elif geom.geom_type == 'MultiPolygon':
        centroid = geom.centroid
    elif geom.geom_type == 'Point':
        centroid = geom
    elif geom.geom_type == 'MultiPoint':
        centroid = geom.centroid
    elif geom.geom_type == 'LineString':
        centroid = geom.centroid
    else:
        continue
    
    name = row.get('name', None)
    if name is None or (isinstance(name, float) and np.isnan(name)):
        name = f"Hospital_{idx}"
    
    hospital_records.append({'hospital_name': str(name), 'geometry': centroid})

hospitals_pts = gpd.GeoDataFrame(hospital_records, crs='EPSG:4326')

print(f"Processed {len(hospitals_pts)} hospital points")

# Get nearest nodes
print("Finding nearest graph nodes...")
incident_nodes = []
for _, row in incidents_df.iterrows():
    node = ox.distance.nearest_nodes(G, row['longitude'], row['latitude'])
    incident_nodes.append(node)

hospital_nodes = []
for _, row in hospitals_pts.iterrows():
    node = ox.distance.nearest_nodes(G, row.geometry.x, row.geometry.y)
    hospital_nodes.append(node)

# For each incident, compute shortest paths
results_closest = []
results_matrix = []

for i, (inc_id, inc_node) in enumerate(zip(incidents_df['incident_id'], incident_nodes)):
    print(f"Processing incident {inc_id} ({i+1}/{len(incidents_df)})...")
    distances = []
    for j, (hosp_name, hosp_node) in enumerate(zip(hospitals_pts['hospital_name'], hospital_nodes)):
        try:
            # Travel time
            tt = nx.shortest_path_length(G, inc_node, hosp_node, weight='travel_time_min')
            # Distance in meters
            dist_m = nx.shortest_path_length(G, inc_node, hosp_node, weight='length')
            # Get path nodes
            path_nodes = nx.shortest_path(G, inc_node, hosp_node, weight='travel_time_min')
            
            # Build LineString geometry
            coords = []
            for u, v in zip(path_nodes[:-1], path_nodes[1:]):
                edge_data = G.get_edge_data(u, v)
                if edge_data:
                    # Take first key
                    k = list(edge_data.keys())[0]
                    geom = edge_data[k].get('geometry', None)
                    if geom:
                        # Add coords from edge geometry (skip first point to avoid dupes)
                        edge_coords = list(geom.coords)
                        if not coords:
                            coords.extend(edge_coords)
                        else:
                            # Skip first coord (already the last one in coords)
                            coords.extend(edge_coords[1:])
                    else:
                        if not coords or coords[-1] != (G.nodes[u]['x'], G.nodes[u]['y']):
                            coords.append((G.nodes[u]['x'], G.nodes[u]['y']))
                        coords.append((G.nodes[v]['x'], G.nodes[v]['y']))
            
            route_line = LineString(coords)
            
            distances.append({
                'incident_id': inc_id,
                'hospital_name': hosp_name,
                'network_distance_m': dist_m,
                'travel_time_min': tt,
                'route_geom': route_line
            })
        except nx.NetworkXNoPath:
            print(f"  No path from {inc_id} to {hosp_name}")
        except Exception as e:
            print(f"  Error: {e}")
    
    # Sort by travel time
    distances.sort(key=lambda x: x['travel_time_min'])
    
    if not distances:
        print(f"  No reachable hospitals for {inc_id}!")
        continue
    
    # Closest hospital
    closest = distances[0]
    results_closest.append(closest)
    
    # Top 3 for matrix
    for rank, d in enumerate(distances[:3], 1):
        results_matrix.append({
            'incident_id': d['incident_id'],
            'hospital_name': d['hospital_name'],
            'rank': rank,
            'network_distance_m': d['network_distance_m']
        })

# Build closest_hospital gdf
closest_gdf = gpd.GeoDataFrame(results_closest, crs='EPSG:4326', geometry='route_geom')
closest_gdf = closest_gdf[['incident_id', 'hospital_name', 'network_distance_m', 'route_geom']]
closest_gdf = closest_gdf.rename_geometry('geometry')

# Build distance_matrix gdf
matrix_gdf = gpd.GeoDataFrame(results_matrix, crs='EPSG:4326')
matrix_gdf['geometry'] = None

# Compute isochrones
print("Computing 15-minute isochrones...")
isochrone_results = []
for j, (hosp_name, hosp_node) in enumerate(zip(hospitals_pts['hospital_name'], hospital_nodes)):
    print(f"  Isochrone for {hosp_name}...")
    try:
        iso = ox.isochrones_from_graph(G, [hosp_node], trip_times=[15], weight='travel_time_min')
        if len(iso) > 0:
            poly = iso.geometry.iloc[0]
            # Ensure MultiPolygon
            if poly.geom_type == 'Polygon':
                poly = MultiPolygon([poly])
            isochrone_results.append({
                'hospital_name': hosp_name,
                'travel_time_min': 15,
                'geometry': poly
            })
    except Exception as e:
        print(f"  Error generating isochrone for {hosp_name}: {e}")

isochrone_gdf = gpd.GeoDataFrame(isochrone_results, crs='EPSG:4326')

# Project to French official metric CRS
target_crs = 'EPSG:2154'

incidents_gdf = gpd.GeoDataFrame(
    incidents_df,
    geometry=[Point(x, y) for x, y in zip(incidents_df['longitude'], incidents_df['latitude'])],
    crs='EPSG:4326'
)

incidents_gdf = incidents_gdf.to_crs(target_crs)
closest_gdf = closest_gdf.to_crs(target_crs)
matrix_gdf = matrix_gdf.to_crs(target_crs)
isochrone_gdf = isochrone_gdf.to_crs(target_crs)

# Write to GPKG
output_path = '/work/emergency_routing.gpkg'

print(f"Writing {output_path}...")

incidents_gdf.to_file(output_path, layer='incidents', driver='GPKG')
closest_gdf.to_file(output_path, layer='closest_hospital', driver='GPKG')
matrix_gdf.to_file(output_path, layer='distance_matrix', driver='GPKG')
isochrone_gdf.to_file(output_path, layer='isochrones_15min', driver='GPKG')

print("Done!")
```

Let me also check what version of osmnx is available and what APIs are available. Let me test.

Let me write the comprehensive Python script that:

1. Reads incidents from CSV
2. Downloads road network from OSM via Overpass for the bbox 48.83,2.30 to 48.88,2.38
3. Downloads hospitals from OSM via Overpass
4. Computes shortest routes from each incident to its closest hospital
5. Computes a distance matrix for each incident vs its three nearest hospitals
6. Computes 15-minute drive-time isochrones around each hospital
7. Writes everything to a GPKG with the specified layer names

Key considerations:
- CRS for Paris: EPSG:2154 (Lambert-93, the official metric projection for metropolitan France)
- Speed: 30 km/h default, use posted maxspeed otherwise
- Need networkx for shortest path
- Need osmnx for graph building, isochrones, routing

Let me think about how to do this:

1. Use osmnx to download the street network for the bounding box
2. Use osmnx's features_from_bbox to get hospitals (amenity=hospital)
3. Project to EPSG:2154
4. For each incident, find the nearest hospital using network distance
5. Compute routes using nx.shortest_path
6. For distance matrix, compute the shortest path distances to all hospitals, rank them, take top 3
7. For isochrones, use osmnx's isochrone function with 15 min travel time at the defined speeds

Let me write the script carefully.

Actually, let me think about this more carefully.

For the Overpass query, osmnx can handle the download. Let me use:
- `osmnx.graph_from_bbox(48.83, 48.88, 2.38, 2.30, network_type='drive')` for the road network
- `osmnx.features_from_bbox(48.83, 48.88, 2.38, 2.30, tags={'amenity': 'hospital'})` for hospitals

Wait, but osmnx graph_from_bbox uses (north, south, east, west). The bbox is 48.83,2.30 to 48.88,2.38. So north=48.88, south=48.83, east=2.38, west=2.30.

For speed assignment: osmnx graphs have `maxspeed` attributes on edges. I need to convert these to numeric speed (km/h), defaulting to 30 km/h if missing.

For isochrones, osmnx has `ox.isochrones()` which takes a graph, node, and travel time.

Let me plan the script:

```python
import osmnx as ox
import geopandas as gpd
import pandas as pd
import networkx as nx
import numpy as np
from shapely.geometry import Point, LineString, MultiLineString
import warnings
warnings.filterwarnings('ignore')

# CRS for Paris - official metric system
CRS_METRIC = "EPSG:2154"  # Lambert-93

# 1. Read incidents
incidents_df = pd.read_csv('/work/incidents.csv')
incidents_gdf = gpd.GeoDataFrame(
    incidents_df,
    geometry=gpd.points_from_xy(incidents_df.longitude, incidents_df.latitude),
    crs="EPSG:4326"
).to_crs(CRS_METRIC)

# 2. Download road network
north, south, east, west = 48.88, 48.83, 2.38, 2.30
G = ox.graph_from_bbox(north, south, east, west, network_type='drive', simplify=True)

# Project graph
G_proj = ox.project_graph(G, to_crs=CRS_METRIC)

# Add speed and travel time to edges
for u, v, k, data in G_proj.edges(keys=True, data=True):
    maxspeed = data.get('maxspeed')
    if maxspeed is None:
        speed = 30  # km/h default
    elif isinstance(maxspeed, list):
        # Take the max if list
        speeds = []
        for s in maxspeed:
            try:
                speeds.append(float(s.split()[0]))
            except:
                pass
        speed = max(speeds) if speeds else 30
    else:
        try:
            speed = float(str(maxspeed).split()[0])
        except:
            speed = 30
    data['speed_kph'] = speed
    # travel time in minutes
    data['travel_time_min'] = (data['length'] / 1000) / speed * 60

# 3. Download hospitals
hospitals = ox.features_from_bbox(north, south, east, west, tags={'amenity': 'hospital'})
# Filter to points or polygons
hospitals_gdf = hospitals.copy()
# Get centroid for polygon hospitals
hospitals_gdf['geometry'] = hospitals_gdf.geometry.centroid
hospitals_gdf = hospitals_gdf[['name', 'geometry']].copy()
hospitals_gdf = hospitals_gdf.to_crs(CRS_METRIC)
# Drop those without name
hospitals_gdf = hospitals_gdf.dropna(subset=['name'])
hospitals_gdf = hospitals_gdf.reset_index(drop=True)

# 4. Find nearest nodes for incidents and hospitals
incident_nodes = []
for idx, row in incidents_gdf.iterrows():
    node = ox.nearest_nodes(G_proj, row.geometry.x, row.geometry.y)
    incident_nodes.append(node)

hospital_nodes = []
for idx, row in hospitals_gdf.iterrows():
    node = ox.nearest_nodes(G_proj, row.geometry.x, row.geometry.y)
    hospital_nodes.append(node)

# 5. Compute closest hospital for each incident
# ... use nx.shortest_path and compute distances
```

Wait, I need to be careful. ox.nearest_nodes expects (X, Y) in the graph's CRS (which is already EPSG:2154). That should work.

Let me think about the isochrones. osmnx has `ox.isochrones` which uses `ox.graph_to_gdfs` and networkx. But actually, `ox.isochrones` expects a graph and returns isochrone polygons. Let me check the API.

Actually, in osmnx 2.x, `ox.isochrones` takes the projected graph, node centers, and edge travel times in minutes, then returns polygons.

Let me use `ox.isochrones_from_graph(G_proj, hospital_nodes, travel_times=[15])` or similar.

Actually, looking at the osmnx docs more carefully:

```python
ox.isochrones(G, nodes, edge_buff=25, node_buff=0, clean_intersects=True)
```

No, that's not right. Let me think...

Actually, the modern osmnx API for isochrones is:

```python
# After computing travel times on edges
# Use networkx.single_source_dijkstra_path_length to get reachable nodes within time
# Then use osmnx's `ox.graph_to_gdfs` to get geometries and create isochrone polygons
```

Or I can use `ox.isochrones` which does this:

```python
# From osmnx docs:
# ox.isochrones(G, nodes=[node], travel_times=[15], edge_buff=25)
```

Let me just write the script and handle it carefully.

Actually, in osmnx 2.x, the function is `ox.isochrones_map` or something... Let me check what's available.

Actually, let me think about this differently. I'll use networkx to compute the shortest paths and distances, then manually construct the isochrones.

For isochrones:
1. For each hospital, find the nearest graph node
2. Use nx.single_source_dijkstra_path_length to get all nodes within 15 minutes travel time
3. Get the subgraph of those nodes
4. Use osmnx to create a convex hull or concave hull around those nodes' geometries to form the isochrone polygon

Actually, ox has a function for this. Let me look:

```python
import osmnx as ox
help(ox.isochrones)
```

Let me just write the script and test it. I'll use `ox.isochrones` from the graph module.

Actually, I just recalled - in osmnx 2.x, `ox.isochrones` is available in `ox.isochrones` module. The function signature is:

```python
ox.isochrones(G, nodes, travel_times, edge_buff=25, node_buff=0, clean_intersects=True)
```

Where `nodes` are the graph node IDs, and `travel_times` is a list of times in minutes.

Let me write a comprehensive script.

Actually wait, let me check what's available more carefully.

Let me just write the script and run it; if something fails I'll fix it.

Here's my plan for the script:

```python
#!/usr/bin/env python3
"""
Emergency routing analysis for Paris incidents.
Produces emergency_routing.gpkg with 4 layers.
"""

import osmnx as ox
import geopandas as gpd
import pandas as pd
import networkx as nx
from shapely.geometry import Point, LineString, MultiPolygon
import warnings
warnings.filterwarnings('ignore')

# Configuration
CRS_WGS84 = "EPSG:4326"
CRS_METRIC = "EPSG:2154"  # Lambert-93, official metric for France
OUTPUT_GPKG = "/work/emergency_routing.gpkg"
DEFAULT_SPEED_KPH = 30
ISOCHRONE_TIME = 15  # minutes
BBOX = (48.88, 48.83, 2.38, 2.30)  # north, south, east, west

# 1. Load incidents
incidents_df = pd.read_csv("/work/incidents.csv")
incidents_gdf = gpd.GeoDataFrame(
    incidents_df,
    geometry=gpd.points_from_xy(incidents_df.longitude, incidents_df.latitude),
    crs=CRS_WGS84
).to_crs(CRS_METRIC)

# 2. Download road network from OSM
print("Downloading road network...")
G = ox.graph_from_bbox(*BBOX, network_type='drive', simplify=True)
G_proj = ox.project_graph(G, to_crs=CRS_METRIC)

# Add speed and travel time to edges
for u, v, k, data in G_proj.edges(keys=True, data=True):
    ms = data.get('maxspeed')
    if ms is None:
        speed = DEFAULT_SPEED_KPH
    elif isinstance(ms, list):
        speeds = []
        for s in ms:
            try:
                speeds.append(float(str(s).split()[0]))
            except:
                pass
        speed = max(speeds) if speeds else DEFAULT_SPEED_KPH
    else:
        try:
            speed = float(str(ms).split()[0])
        except:
            speed = DEFAULT_SPEED_KPH
    data['speed_kph'] = speed
    # travel time in minutes: length(m) / 1000 / speed(km/h) * 60
    data['travel_time_min'] = (data['length'] / 1000.0) / speed * 60.0

# 3. Download hospitals
print("Downloading hospitals...")
hospitals = ox.features_from_bbox(*BBOX, tags={'amenity': 'hospital'})
# Keep only needed columns, use centroid
hospitals_gdf = hospitals[['name', 'geometry']].copy()
hospitals_gdf = hospitals_gdf.to_crs(CRS_METRIC)
# Replace polygon/multipolygon geometries with centroids
hospitals_gdf['geometry'] = hospitals_gdf.geometry.centroid
hospitals_gdf = hospitals_gdf.dropna(subset=['name']).reset_index(drop=True)
# Also drop any empty geometry
hospitals_gdf = hospitals_gdf[~hospitals_gdf.geometry.is_empty]

# 4. Map incidents and hospitals to nearest graph nodes
print("Mapping to nearest graph nodes...")
incident_pts = [(row.geometry.x, row.geometry.y) for _, row in incidents_gdf.iterrows()]
hospital_pts = [(row.geometry.x, row.geometry.y) for _, row in hospitals_gdf.iterrows()]

incident_nodes = [ox.nearest_nodes(G_proj, X=x, Y=y) for x, y in incident_pts]
hospital_nodes = [ox.nearest_nodes(G_proj, X=x, Y=y) for x, y in hospital_pts]

# Build a travel-time-weighted graph for routing
# NetworkX shortest paths use the weight attribute
for _, _, _, data in G_proj.edges(keys=True, data=True):
    data['weight'] = data['travel_time_min']

# 5. Compute closest hospital for each incident
print("Computing closest hospitals and routes...")
closest_rows = []
distance_rows = []

for i_idx, inc_node in enumerate(incident_nodes):
    inc_id = incidents_gdf.iloc[i_idx]['incident_id']
    
    # Compute distances to all hospitals
    dists = []
    for h_idx, hos_node in enumerate(hospital_nodes):
        try:
            # Shortest path distance in travel time (minutes)
            dist_time = nx.shortest_path_length(G_proj, inc_node, hos_node, weight='weight')
            # Also get distance in metres (using length)
            path = nx.shortest_path(G_proj, inc_node, hos_node, weight='weight')
            dist_m = sum(
                G_proj.edges[path[i], path[i+1], 0].get('length', 0)
                for i in range(len(path)-1)
            )
            hospital_name = hospitals_gdf.iloc[h_idx]['name']
            dists.append((h_idx, hospital_name, dist_m, dist_time, path))
        except (nx.NetworkXNoPath, Exception):
            continue
    
    if not dists:
        continue
    
    # Sort by distance (metres)
    dists.sort(key=lambda x: x[2])
    
    # Closest hospital - one row
    h_idx, name, dist_m, dist_time, path = dists[0]
    
    # Build route geometry
    route_edges = []
    for i in range(len(path)-1):
        u, v = path[i], path[i+1]
        # Get edge data
        edge_data = G_proj.get_edge_data(u, v)
        if edge_data:
            # Pick the first edge if multiple
            k = list(edge_data.keys())[0]
            geom = edge_data[k].get('geometry')
            if geom is None:
                # Create line from node coordinates
                geom = LineString([
                    (G_proj.nodes[u]['x'], G_proj.nodes[u]['y']),
                    (G_proj.nodes[v]['x'], G_proj.nodes[v]['y'])
                ])
            route_edges.append(geom)
    
    route_geom = LineString()
    if len(route_edges) == 1:
        route_geom = route_edges[0]
    elif len(route_edges) > 1:
        # Merge linestrings
        from shapely.ops import linemerge
        route_geom = linemerge(route_edges)
    
    closest_rows.append({
        'incident_id': inc_id,
        'hospital_name': name,
        'network_distance_m': round(dist_m, 2),
        'geometry': route_geom
    })
    
    # Distance matrix - top 3
    for rank, (_, name, dist_m, _, _) in enumerate(dists[:3], 1):
        distance_rows.append({
            'incident_id': inc_id,
            'hospital_name': name,
            'rank': rank,
            'network_distance_m': round(dist_m, 2)
        })

closest_gdf = gpd.GeoDataFrame(closest_rows, crs=CRS_METRIC)
distance_df = pd.DataFrame(distance_rows)

# 6. Compute 15-min isochrones for each hospital
print("Computing isochrones...")
isochrone_rows = []

for h_idx, hos_node in enumerate(hospital_nodes):
    hospital_name = hospitals_gdf.iloc[h_idx]['name']
    
    try:
        # Use ox.isochrones
        isochrone_polys = ox.isochrones(
            G_proj, 
            nodes=[hos_node], 
            travel_times=[ISOCHRONE_TIME],
            edge_buff=25,
            node_buff=0
        )
        if not isochrone_polys.empty:
            # Take the first (and only) polygon for this time
            poly = isochrone_polys.iloc[0].geometry
            isochrone_rows.append({
                'hospital_name': hospital_name,
                'travel_time_min': ISOCHRONE_TIME,
                'geometry': poly
            })
    except Exception as e:
        print(f"  Isochrone failed for {hospital_name}: {e}")
        continue

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

# 7. Write to GPKG
print("Writing GPKG...")

# Layer 1: incidents
incidents_out = incidents_gdf[['incident_id', 'geometry']].copy()
incidents_out.to_file(OUTPUT_GPKG, layer='incidents', driver='GPKG')

# Layer 2: closest_hospital
closest_gdf.to_file(OUTPUT_GPKG, layer='closest_hospital', driver='GPKG')

# Layer 3: distance_matrix (tabular)
distance_gdf = gpd.GeoDataFrame(distance_df, geometry=[None]*len(distance_df), crs=CRS_METRIC)
distance_gdf.to_file(OUTPUT_GPKG, layer='distance_matrix', driver='GPKG')

# Layer 4: isochrones_15min
isochrones_gdf.to_file(OUTPUT_GPKG, layer='isochrones_15min', driver='GPKG')

print("Done!")
```

Hmm, let me think about a few things:

1. The `ox.nearest_nodes` function signature - it takes (G, X, Y) where X is longitude and Y is latitude. But since we've projected the graph, we need to pass the projected coordinates.

2. For `ox.isochrones` - let me check. In osmnx 2.x, the function might be in a different place. Let me look at the actual API.

Actually, I think in osmnx 2.x, the function is `ox.isochrones` but it takes a graph projected, node IDs, and travel times. Let me just write the script and test.

Let me also think about edge cases with multi-edges. When I do `nx.shortest_path(G_proj, u, v, weight='weight')`, it gives me a list of nodes. Then I need to look up the edge between consecutive nodes, but there might be multiple parallel edges. I need to handle that.

Let me use a more robust approach for getting edge geometries along a path.

Actually, let me simplify: I can use `ox.utils_graph.graph_to_gdfs` to get the edges geodataframe, then merge route edges. But that's complex.

A simpler approach: after getting the path (list of nodes), I can use `ox.plot_graph_route` type logic, or just get the edge geometry by finding the edge with the right geometry.

Let me write a cleaner version:

```python
def get_route_geometry(G, path):
    """Get the merged LineString geometry for a path through the graph."""
    lines = []
    for u, v in zip(path[:-1], path[+1:]):
        # Try each parallel edge
        edge_data = G.get_edge_data(u, v)
        if edge_data is None:
            continue
        # Pick the edge with smallest weight (travel time)
        best_k = min(edge_data.keys(), key=lambda k: edge_data[k].get('weight', float('inf')))
        geom = edge_data[best_k].get('geometry')
        if geom is None:
            geom = LineString([
                (G.nodes[u]['x'], G.nodes[u]['y']),
                (G.nodes[v]['x'], G.nodes[v]['y'])
            ])
        lines.append(geom)
    
    if len(lines) == 1:
        return lines[0]
    elif len(lines) > 1:
        from shapely.ops import linemerge
        return linemerge(lines)
    return None
```

OK let me just write the full script and run it.

One more thing: for the `distance_matrix` layer, it says "geometry may be empty or null". So I can use None geometry.

Let me write and test the script now.

Actually, for isochrones in osmnx 2.x, let me just check if it exists:

I'll use a try-except or just check. Let me write the script and see.

I need to replace the isochrones section with a manual implementation using networkx's ego_graph or shortest path tree to find all nodes reachable within 15 minutes, then construct a convex hull or concave hull around those nodes, possibly with edge buffering.

Let me write a function that does this:

1. For each hospital node, compute the subgraph of nodes reachable within 15 minutes using Dijkstra
2. Collect all nodes and edges within that threshold
3. Buffer the edges to create a polygon via convex hull or union of buffered geometries

Let me also remove the unused imports I had before and also update the function signature for `route_geometry` which uses `x` and `y` on nodes - in osmnx 2.x projected graphs have `x` and `y` attributes on nodes. Let me verify this works.

Also let me reconsider the isochrone approach. The unary_union of buffered edges approach might be slow and produce a blob rather than a clean polygon. A better approach is to:

1. Find all nodes reachable within 15 minutes
2. For each reachable edge (where both endpoints are reachable), take its geometry
3. Buffer and union

But even better: for proper isochrones we typically want to use the approach where we get the whole subgraph within the cutoff and then use the concave hull of the edge nodes. Let me simplify and use a standard approach.

Actually, the revised approach using `nx.single_source_dijkstra_path_length` with all reachable edges buffered and unioned should work fine. Let me run it and see.

But first, I also need to make sure the `route_geometry` function works with projected graph nodes. In osmnx, projected graphs store node coords in `x` and `y` attributes.

Let me also re-check the script - oh wait, I already replaced the isochrone section with one that uses `G_proj.nodes[u]["x"]` - those exist in osmnx projected graphs.

Let me just run it now.

The command timed out after 5 minutes. The Overpass API might be very slow or the requests might be getting stuck. Let me think of an alternative approach.

Instead of using osmnx which relies on Overpass API, I could:

1. Download the OSM data from a Geofabrik extract (Paris/Île-de-France) - but that's a large file
2. Use DuckDB with Overture Maps to get roads and hospitals
3. Use Overture Maps Python library

Let me try using Overture Maps instead. The overturemaps library is pre-installed. Let me check what's available.

Actually, let me first check if we can reach external services via curl or wget.

Let me think about this. Overpass API is not available from this environment. Overture Maps has roads (segment), hospitals (place?), buildings, etc. But Overture's road network (segment) is not the same as a routable graph with connectivity/directionality. 

Let me try another approach - maybe I can download OSM data as a PBF file from Geofabrik and use osmium to extract the data.

Actually, let me check if `curl` or `wget` can reach any of these endpoints:

Let me try a different approach. Let me check if the overpass endpoints work with Python's urllib or requests:

Actually, the issue seems to be network connectivity. Let me try using Overture Maps data since it works via S3 (different servers).

Actually, let me check if there's another Overpass mirror that might work.

Let me build a new approach using Overture Maps for both the road network (segment type) and hospitals (place type). The Overture Maps data has:

- `segment` - for roads/transportation
- `place` - for POIs including hospitals

Let me look at how to use Overture Maps effectively.

Good, so I can use Overture Maps for both the road network (segment subtype='road') and hospitals/places. Let me think about the approach:

1. **Segments (roads)**: Filter `segment` for `subtype == 'road'`. These are line geometries representing road segments. I'll need to build a graph from them. Each segment has `speed_limits` and geometry.

2. **Hospitals**: Use `place` data. Filter for categories that include "hospital" or use the building data with `amenity=hospital`.

3. **Incidents**: Already in CSV.

Let me look at the place data more carefully to find hospitals, and also at the segment data for speed limits.

Let me use `overturemaps.geodataframe` to get the data more easily. The issue is that Overture Maps has its own schema for places - I need to find hospitals specifically. Let me look at the categories structure more closely.

Also for the road network, I need to use the `segment` type with `subtype == 'road'` to get drivable roads. For hospitals, I need to filter places by category.

Let me write a new approach using Overture Maps data:

OK so Overture Maps works. But the issue is that Overture has segments (roads) but they don't form a graph. I'd need to construct a routable graph from the segments. That's quite complex.

Let me think about another approach. Since I can't reach Overpass directly from this environment, maybe I should try a different tactic - use the Overture segment data to build a network, or maybe download data from elsewhere.

Actually, let me try one more time with Overpass via `curl` if it's available... no, curl isn't available. Let me try `wget`.

The Overpass endpoints are unreachable from this environment (network is unreachable). But Overture Maps works because it uses S3. Let me now think about the full approach:

1. **Incidents**: I have them from CSV, easy.
2. **Road network**: Use Overture Maps `segment` data (subtype='road' for driveable roads). But I need to build a routable graph from that.
3. **Hospitals**: Use Overture Maps `place` data with category primary='hospital'.
4. **Routing**: Build a graph from Overture road segments, then run shortest paths.
5. **Isochrones**: Use Dijkstra and buffer reachable edges.

Overture segments have `speed_limits` attribute for speed limits. Let me check the structure.

Actually, let me think more carefully. Building a routable graph from Overture segments is more complex than using osmnx because:
- Overture segments are just linestrings with metadata, not a graph with connectivity.
- I'd need to build a graph from scratch, handling intersections properly.

Alternatively, I could try using DuckDB with the spatial extension to read the Overture data directly from S3, which might be faster.

But actually, the simplest approach given the constraints is:
1. Download Overture road segments
2. Build a networkx graph from them (using endpoints as nodes)
3. Download Overture hospitals
4. Run routing

Let me check the segment structure more carefully - specifically speed_limits and how to build a graph.

Now I have a clearer picture. Overture Maps provides:
1. `segment` (subtype='road') with LineString geometry, connectors, speed_limits
2. `connector` with Point geometry (intersection/junction nodes)
3. `place` with categories including 'hospital' as `primary`

The connectors define the topology - each segment has connector references at positions along the line, and connectors are points that link segments together.

Let me build the solution using Overture Maps data instead of Overpass. I'll need to:
1. Build a routable graph from segments and connectors
2. Get hospital locations from places
3. Compute shortest paths and isochrones

Let me check how connectors link segments together.

The Overture Maps data has segments with connectors. This is a bit complex to build a graph from. Let me think about an alternative approach.

Since the environment can't reach Overpass, I have a few options:

1. Use Overture Maps data to build a road network graph - this is complex but doable
2. Try a different approach - use DuckDB to read Overture Maps parquet data directly
3. Use a local OSM PBF extract if available

Actually, let me think about the Overture approach more carefully. The segments have `connectors` (points at positions along the geometry) and there's a separate `connector` layer. This is the Overture transportation schema - it's a topological network.

Let me build a graph from road segments and connectors. Each connector is a node, each segment edge connects connectors at positions along the geometry.

Actually, this is getting complex. Let me try a different approach - use DuckDB to query Overture Maps data more efficiently, or see if I can use a smaller approach.

Actually, let me reconsider. The instructions say:
- "When fetching from an external API (Overpass, Overture, Geofabrik, etc.)"
- "overturemaps: Overture Maps download (Python + CLI). Wraps the s3://overturemaps-us-west-2/ GeoParquet bucket"

So Overture Maps is a valid data source. Let me build the solution using Overture Maps segments and connectors.

The approach:
1. Load road segments (subtype='road') from Overture within bbox
2. Load connectors from Overture
3. Build a routing graph where connectors are nodes and road segments are edges
4. Load hospitals from Overture places
5. Load incidents from CSV
6. Compute routes, distances, isochrones

The tricky part is building the graph from Overture data. Let me think about it:

For each road segment:
- It has connectors at various positions along its geometry (at fraction 0.0 to 1.0)
- Each connector has an ID that matches to the connector table
- Each connector has a Point geometry (lat/lon)
- Road segments connect pairs of adjacent connectors along the geometry

So I need to:
1. For each road segment, get its connectors sorted by position
2. For each adjacent pair of connectors, create a graph edge with the sub-geometry between them and length/speed info
3. Use connector geometries as node locations

Let me code this up properly.

Let me use a different approach. Since Overpass is unreachable, I'll use Overture Maps data which is working. The Overture segments have connectors (nodes) and geometry. Let me build a routable graph from the road segments.

The segment data has:
- `id`: segment id
- `geometry`: LineString
- `connectors`: array of {connector_id, at} dicts
- `class`: road classification
- `subtype`: road
- `road_flags`: array for oneway, tunnel, etc.
- `speed_limits`: array of {min_speed, max_speed, ...}

And connectors table has:
- `id`: connector id
- `geometry`: Point

Building a proper routable graph from Overture segments is complex but doable. Let me think about this...

The connectors table has point geometries. The segments reference connectors at certain positions along the geometry (the `at` value, 0 to 1). The first and last connectors (at 0.0 and 1.0) are the endpoints of the segment.

For routing, I need:
1. NetworkX graph where nodes = connectors, edges = segments between consecutive connectors
2. Each segment can be split at its connector points

Actually, looking at the connectors more carefully, the `at` values represent the fraction along the segment geometry. So a segment from connector A (at 0.0) to connector B (at 1.0) might have intermediate connectors.

But for routing, the key topology is: segments connect to each other at shared connectors. A road segment goes from its first connector (at=0.0) to its last connector (at=1.0), possibly with intermediate connectors.

Let me simplify: build a graph where each segment creates edges between consecutive connector pairs, with geometry extracted from the sub-line between those connectors.

Actually, let me take an even simpler approach: each road segment connects its first connector (at=0.0) to its last connector (at=1.0). The intermediate connectors are just shape points. The distance is the full segment length. This works for routing because the intermediate connectors don't connect to other segments (the topology is at the endpoints).

Wait, but intermediate connectors CAN connect to other segments. Let me check...

Actually let me look at this more carefully. In Overture, connectors are at intersections. A segment's connectors list includes all connectors that lie on that segment. If the same connector_id appears in multiple segments' connector lists, that's an intersection.

So I need to treat each connector as a node, and each consecutive pair within a segment's sorted connectors as an edge.

Let me also check: does the segment geometry have vertices that correspond to the connector positions?

Let me build this step by step.

OK, so Overture Maps for this region doesn't have oneway road flags - interesting. Let me check the `access_restrictions` field and also check if connectors have the info we need to build a graph.

Since Overpass is unreachable, I'll use Overture Maps for:
1. Roads (segment subtype 'road') - build a routable graph
2. Hospitals (place with category 'hospital') 
3. Use networkx for shortest paths

Let me also check if overturemaps has a `segment` connector that links them.

For the routing graph, I'll need to:
1. Get all road segments (LineStrings)
2. Get all connectors (Points)
3. Build a graph where connectors are nodes and road segments are edges
4. For each edge, compute length in meters and travel time

Let me check how connectors link to segments.

Now I have a good understanding of Overture Maps data. Let me build the solution:

1. **Incidents** - from CSV, simple
2. **Road network** - from Overture Maps segments with subtype='road', filter to drivable classes
3. **Hospitals** - from Overture Maps places with category 'hospital'
4. **Routing** - Build a NetworkX graph from the road network using connectors as nodes
5. **Isochrones** - Compute travel time from each hospital

Let me design the approach:

For the road network, Overture Maps segments have:
- `id`: segment ID
- `geometry`: LineString in WGS84
- `connectors`: list of {connector_id, at} at fractional positions along the geometry
- `class`: road class (residential, primary, etc.)
- `speed_limits`: optional list of speed limits (with value and unit)
- `road_flags`: optional flags

For routing, I need to:
1. Build a graph where nodes = connectors (from the connectors table joined with segments)
2. Edges = segment portions between consecutive connectors
3. Use speed limits to compute travel times

But wait - the connectors table has point geometries. Each connector_id is a unique node. Each segment has a list of connectors along it with fractional positions. So I need to:
- For each connector in a segment, get its geometry from the connectors table
- Create edges between consecutive connectors on the same segment
- Use the sub-geometry and sub-length of the segment for each edge
- Apply the speed limit

Let me think about the driveable road classes. From the class distribution:
- residential, service, tertiary, secondary, primary, living_street, unclassified - these are driveable
- footway, steps, cycleway, pedestrian, path - not driveable
- unknown - skip

Let me also check what the 'road' filter already gives vs the class distribution. The subtype='road' already filters some things. Let me see what classes are considered road vs non-road in Overture.

Actually, I see that the data already has `subtype='road'` which gives 36880 segments. But many are footway, steps, etc. Let me filter to only drivable classes.

For the graph building, I need to be more careful. Let me look at how the connectors work.

Actually, I think a simpler approach would work: use the full road segment geometries directly, connecting segments that share connectors. But that requires building a proper graph.

Let me think about this differently. The Overture Maps connectors are the junction points. Each connector_id is a unique intersection point. Each segment references connectors at fractional positions along its geometry.

So the graph structure is:
- Nodes: connector points (from connectors table)
- Edges: portions of road segments between consecutive connectors

For each road segment:
1. Get its connectors, sorted by `at` position
2. For each consecutive pair of connectors, create an edge (sub-line from position at_i to at_{i+1})
3. Compute the sub-length, apply speed limit, compute travel time

This is the correct approach but complex. Let me implement it carefully.

Actually, for routing purposes, I don't need to split segments at every connector - I can just use the full segment between its first and last connector as one edge. The connectors at intermediate points would be other roads connecting at those points. But this won't create proper junctions unless we actually split.

Let me take the simpler approach: treat each road segment as an edge between its first and last connector (the endpoints), and use all connectors as potential junction nodes. This is how Overture is designed - segments connect at connectors, and multiple segments can reference the same connector.

Wait, but some segments have intermediate connectors (at positions not 0.0 or 1.0). These represent places where other roads connect along the segment. For routing purposes I should split the segment at those points.

Let me implement a cleaner approach:

1. For each road segment, get its sorted connectors
2. For each consecutive pair (c_i, c_{i+1}), create a graph edge
3. Get the sub-geometry by interpolating along the segment LineString
4. Compute sub-length from the geometry

This is doable. Let me code it up.

Actually, for the purposes of this problem, a simpler approach would work well enough. Let me:

1. Get all road segments with drivable classes
2. Extract the start and end points of each segment (from the geometry)
3. Build a graph where edges are segments and nodes are the start/end points
4. For junctions, use the connectors to identify which segments share a node

Wait, the connectors already handle this. Let me use the full connector table.

Let me simplify even further: I'll treat each road segment as a graph edge connecting its endpoint connectors (the connectors at at=0.0 and at=1.0). For intermediate connectors, I'll also split the segment. This is complex but correct.

Let me just write the script.

Wait, I left references to `_interpolate_subline` inside the main loop before it was defined. Let me check the flow of the script more carefully.

Actually I defined `_interpolate_subline` after the loop and also had a helper `p_interpolate`. The function definitions at module level are hoisted, so they will be available. But wait, I defined `_interpolate_subline` and `p_interpolate` inside the script flow (not as top-level functions, they're defined after the print statements). In Python, function definitions are executed at runtime, so they need to be defined before they're called. Currently `_interpolate_subline` is called inside the big loop in section 5, but it's defined later in the file. That won't work!

Let me restructure the script to put all helper functions at the top.

The issue is that nodes might not have x/y if they were auto-created via add_edge. Let me also make sure that when we add edges, we ensure the endpoints have x/y set. The simplest fix: after building the graph, go through each node and if it doesn't have x/y, look it up from the connector geometries.

Actually, the issue is clearer now. `G.add_edge(cid_a, cid_b, ...)` will automatically add nodes cid_a and cid_b if they don't exist, but they won't have x/y attributes. The loop at line 244 only adds nodes for sorted_conns, but many nodes are auto-created via edge addition.

Fix: Ensure nodes are created before edges, or after graph building, ensure all nodes have x/y.

The issue is that Overture categorizes many things as 'hospital' that aren't actually hospitals (doctors offices, clinics, etc). But the key point is: 186 hospitals is too many for computing isochrones. Let me take a different approach:

1. The task says "all hospitals" in the bbox from OSM. In Overture, let me filter those with `basic_category == 'hospital'` and `taxonomy.primary == 'hospital'` - that's still 186.

Let me check if there's a way to identify only major hospitals. Actually, let me just limit the isochrones to the hospitals that are actually relevant - i.e., those that show up in the closest or top-3 for any incident. That would reduce the isochrones from 186 to at most 8*3=24. But the task says "every hospital" in the area.

Let me try a more efficient isochrone approach. The current one buffers each edge geometry individually - that's slow for 186 hospitals with 43000 edges each time.

Better approach: 
1. Precompute the network distances from each node to each hospital using multi-source Dijkstra
2. For isochrones, we can do a faster approach: use the graph to find reachable nodes, then compute the convex hull of those node positions, buffered.

Actually, the issue is the buffering of 43000 edges × 186 hospitals = 8 million operations. Let me instead:
- For each hospital, get reachable nodes
- Take the convex hull of those nodes + buffer
- Much faster

Or even simpler: compute reachable nodes, then use a concave hull approach. Let me just buffer the reachable node coordinates.

Also, let me filter to only proper hospitals - those that have meaningful names with "Hôpital", "Hospital", "Clinique" etc.

The problem is 186 hospitals, many of which are not real hospitals. The isochrone computation loops through all edges for each hospital, which is slow. Let me:

1. Filter hospitals to real ones only
2. Use a convex hull of reachable nodes for isochrones instead of buffering every edge

Actually, the main bottleneck is the isochrone for 186 hospitals. Each one requires:
- single_source_dijkstra_path_length (fast)
- Loop through ALL edges (43000 edges × 186 = 8M iterations, each with a set lookup)

And then buffering all those edge geometries.

Let me rewrite the isochrone section to be much more efficient:
1. Use convex hull of reachable node coordinates
2. Filter hospitals to real ones

Actually, 186 isochrones could run in reasonable time if I use a more efficient approach. Let me instead use node convex hulls rather than edge buffering, which will be much faster.

Also, many of these "hospitals" are not real hospitals - they're individual practitioners or clinics mis-categorized. Let me try to filter to names that are more clearly hospitals.

Wait, let me first check if the script actually completed or just printed those names and then timed out during the isochrone computation itself. Looking at the output, it printed ~40 hospital names with checkmarks and then timed out. So it was in the middle of the isochrone loop when it hit the 8-minute timeout.

Let me rewrite with a faster isochrone approach - using convex hulls of reachable nodes instead of edge buffering.