spa-l3-paris-emergency-routing
Map (reference ← swipe → agent)
0:00
Need the dispatch coverage model rebuilt before Thursday's briefing. The eight historical emergency-call locations are in `incidents.csv` (columns: incident_id, latitude, longitude, label). Pull the driveable road network and all hospitals from OpenStreetMap via Overpass for the area bbox 48.83,2.30 to 48.88,2.38. I need three deliverables: the shortest driving route from each incident to its closest hospital, a distance matrix of each incident against its three nearest hospitals, and a 15-minute drive-time isochrone around every hospital to visualise coverage gaps. Assume 30 km/h for edges without a maxspeed tag; use posted maxspeed otherwise.
Everything goes into a single GPKG called `emergency_routing.gpkg` with all layers in the official metric coordinate system for Paris and exactly these four layer names. The `incidents` layer has Point (or MultiPoint) geometry and an `incident_id` column. The `closest_hospital` layer has LineString (or MultiLineString) geometry representing the route, with columns `incident_id`, `hospital_name`, and `network_distance_m` (numeric, metres). One row per incident. The `distance_matrix` layer is tabular (geometry may be empty or null) with columns `incident_id`, `hospital_name`, `rank` (integer 1 through 3, ascending by `network_distance_m` within each incident), and `network_distance_m`. The `isochrones_15min` layer has MultiPolygon geometry (Polygon also accepted) with columns `hospital_name` and `travel_time_min`. One row per hospital.
0:00
System info
{
"adapter": "openrouter",
"model": "google/gemma-4-26b-a4b-it",
"base_url": "https://openrouter.ai/api/v1"
}
0:01
Bash
ls -F
ls -F
incidents.csv
0:01
Bash
head -n 5 incidents.csv
head -n 5 incidents.csv
incident_id,latitude,longitude,label INC_001,48.8566,2.3522,near Notre-Dame INC_002,48.8620,2.3360,near Louvre INC_003,48.8462,2.3464,Quartier Latin INC_004,48.8700,2.3431,near Gare du Nord area
1:10
Write
/work/solve.py
content · 279 lines · py
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point, LineString, MultiPolygon, MultiPoint
import osmnx as ox
import networkx as nx
import numpy as np
# Configuration
BBOX = (48.83, 2.30, 48.88, 2.38) # (min_lat, min_lon, max_lat, max_lon)
OUTPUT_FILE = "emergency_routing.gpkg"
CRS = "EPSG:2154" # RGF93 / Lambert-93 (Official metric for France/Paris)
def solve():
# 1. Load incidents
df_incidents = pd.read_csv('incidents.csv')
# BBox order in osmnx/overpass is usually (north, south, east, west) or similar.
# Let's define the bounds properly.
# lat_min, lon_min, lat_max, lon_max = 48.83, 2.30, 48.88, 2.38
# Convert to GeoDataFrame
gdf_incidents = gpd.GeoDataFrame(
df_incidents,
geometry=gpd.points_from_xy(df_incidents.longitude, df_incidents.latitude),
crs="EPSG:4326"
).to_crs(CRS)
# 2. Get Road Network and Hospitals using OSMnx
# BBox: north, south, east, west
north, south, east, west = 48.88, 48.83, 2.38, 2.30
print("Downloading road network...")
# We'll use the bbox to get the graph.
# osmnx.graph_from_bbox(north, south, east, west, network_type='drive')
G = ox.graph_from_bbox(north=north, south=south, east=east, west=west, network_type='drive')
# Apply speed rules
# OSMnx handles maxspeed if present.
# If maxspeed is missing, we set it to 30 km/h (8.333 m/s)
for u, v, k, data in G.edges(data=True, keys=True):
if 'maxspeed' not in data or data['maxspeed'] is None:
data['maxspeed'] = 30
elif isinstance(data['maxspeed'], str):
# Handle cases like "30 km/h"
try:
data['maxspeed'] = float(data['maxspeed'].split()[0])
except:
data['maxspeed'] = 30
# Calculate travel time (seconds)
# speed is in km/h, we need m/s
for u, v, k, data in G.edges(data=True, keys=True):
speed_mps = data['maxspeed'] / 3.6
data['travel_time'] = data['length'] / speed_mps
print("Downloading hospitals...")
# Get hospitals in the bbox
tags = {'amenity': 'hospital'}
# Since graph_from_bbox is used for network, let's use ox.features_from_bbox for hospitals
hospitals = ox.features_from_bbox(north, south, east, west, tags=tags)
# Filter for points or polygons representing hospitals
# We'll take the centroid for routing if they are polygons
hospitals = hospitals.copy()
hospitals['hospital_name'] = hospitals.get('name', 'Unknown Hospital')
# Clean names
hospitals['hospital_name'] = hospitals['hospital_name'].fillna('Unnamed Hospital')
# Convert to metric and extract points for routing
hospitals_gdf = hospitals.to_crs(CRS)
hospitals_gdf['geometry'] = hospitals_gdf.centroid # Use centroid for routing
# 3. Routing
print("Calculating routes and matrix...")
# Project G to CRS
G_proj = ox.project_graph(G, to_crs=CRS)
# We need to find the nearest nodes in the graph for incidents and hospitals
# Pre-calculate nearest nodes for hospitals
nodes, edges = ox.graph_to_gdfs(G_proj)
# Find nearest nodes for incidents and hospitals
# For incidents
inc_nodes = ox.distance.nearest_nodes(G_proj, gdf_incidents.geometry.x, gdf_incidents.geometry.y)
# For hospitals
hosp_nodes = ox.distance.nearest_nodes(G_proj, hospitals_gdf.geometry.x, hospitals_gdf.geometry.y)
# Map hospital index to node
hosp_info = []
for i, (idx, row) in enumerate(hospitals_gdf.iterrows()):
hosp_info.append({
'hospital_name': row['hospital_name'],
'node': hosp_nodes[i],
'original_index': i
})
hosp_df = pd.DataFrame(hosp_info)
# Lists to store results
closest_routes_data = []
distance_matrix_data = []
for i, inc_row in gdf_incidents.iterrows():
inc_id = inc_row['incident_id']
start_node = inc_nodes[i]
# Find distances to all hospitals
dist_to_hospitals = []
for _, h_row in hosp_df.iterrows():
end_node = h_row['node']
try:
# Shortest path distance (using 'length' edge attribute)
path_length = nx.shortest_path_length(G_proj, start_node, end_node, weight='length')
dist_to_hospitals.append({
'hospital_name': h_row['hospital_name'],
'node': end_node,
'network_distance_m': path_length
})
except nx.NetworkXNoPath:
continue
# Sort hospitals by distance
dist_to_hospitals.sort(key=lambda x: x['network_distance_m'])
# 1. Closest hospital route
if dist_to_hospitals:
closest = dist_to_hospitals[0]
path = nx.shortest_path(G_proj, start_node, closest['node'], weight='length')
# Create LineString for the route
route_coords = []
for node in path:
node_data = G_proj.nodes[node]
route_coords.append((node_data['x'], node_data['y']))
closest_routes_data.append({
'incident_id': inc_id,
'hospital_name': closest['hospital_name'],
'network_distance_m': closest['network_distance_m'],
'geometry': LineString(route_coords)
})
# 2. Distance Matrix (top 3)
for rank, h_dist in enumerate(dist_to_hospitals[:3], 1):
distance_matrix_data.append({
'incident_id': inc_id,
'hospital_name': h_dist['hospital_name'],
'rank': rank,
'network_distance_m': h_dist['network_distance_m'],
'geometry': None
})
else:
# Handle case with no path
distance_matrix_data.append({
'incident_id': inc_id,
'hospital_name': None,
'rank': None,
'network_distance_m': None,
'geometry': None
})
# Convert results to GDFs
gdf_closest_hospital = gpd.GeoDataFrame(closest_routes_data, crs=CRS)
gdf_distance_matrix = gpd.GeoDataFrame(distance_matrix_data, crs=CRS)
# Ensure distance_matrix has correct columns even if empty
for col in ['incident_id', 'hospital_name', 'rank', 'network_distance_m']:
if col not in gdf_distance_matrix.columns:
gdf_distance_matrix[col] = None
# 3. Isochrones (15 min = 900 seconds)
print("Calculating isochrones...")
isochrone_data = []
time_limit = 15 * 60 # 15 minutes in seconds
for _, h_row in hosp_df.iterrows():
h_node = h_row['node']
h_name = h_row['hospital_name']
# Get all nodes reachable within time_limit
subgraph = nx.ego_graph(G_proj, h_node, radius=time_limit, weight='travel_time')
if len(subgraph.nodes) > 0:
# Get the boundary (concave hull or union of node areas)
# A simpler way: union of all edges in the subgraph?
# Or just collect node coordinates and use alpha shape/convex hull
# For a "professional" look, let's use the edges to form a polygon.
# Or use the standard method: union of buffered nodes.
# Efficient way: get all nodes in subgraph, buffer them, then union.
# To avoid a massive union of many circles, we can collect nodes and use a concave hull.
node_points = []
for node in subgraph.nodes:
node_data = G_proj.nodes[node]
node_points.append(Point(node_data['x'], node_data['y']))
if node_points:
# Using a simple buffer approach for robustness
# Union of buffers of nodes
# To be faster and less messy, let's buffer the edges.
edge_geoms = []
for u, v, k, data in subgraph.edges(data=True, keys=True):
# Create line from edge
u_data = G_proj.nodes[u]
v_data = G_proj.nodes[v]
edge_geoms.append(LineString([(u_data['x'], u_data['y']), (v_data['x'], v_data['y'])]))
if edge_geoms:
# Buffer the edges slightly to create a "corridor"
# Then union them.
# Or just buffer the nodes and union.
# Let's try buffering nodes with a small radius (e.g. 20m) and unioning.
# Actually, a better way is to use the edges.
combined_edges = unary_union(edge_geoms)
poly = combined_edges.buffer(30) # 30m buffer around the reachable network
# We need to ensure it's a valid polygon/multipolygon
if not poly.is_empty:
isochrone_data.append({
'hospital_name': h_name,
'travel_time_min': 15.0,
'geometry': poly
})
else:
isochrone_data.append({
'hospital_name': h_name,
'travel_time_min': 15.0,
'geometry': None
})
else:
# Only the node itself is reachable
node_data = G_proj.nodes[h_node]
isochrone_data.append({
'hospital_name': h_name,
'travel_time_min': 15.0,
'geometry': Point(node_data['x'], node_data['y']).buffer(30)
})
else:
isochrone_data.append({
'hospital_name': h_name,
'travel_time_min': 15.0,
'geometry': None
})
else:
isochrone_data.append({
'hospital_name': h_name,
'travel_time_min': 15.0,
'geometry': None
})
gdf_isochrones = gpd.GeoDataFrame(isochrone_data, crs=CRS)
# 4. Write to GPKG
print(f"Writing to {OUTPUT_FILE}...")
# Remove rows with empty geometry for certain layers if needed,
# but requirements say "MultiPolygon geometry ... One row per hospital"
# We'll keep them but clean up.
# Save layers
# Requirement: incident_id, hospital_name, network_distance_m
# We have these in closest_routes_data
# Layer 1: incidents
gdf_incidents.to_file(OUTPUT_FILE, layer='incidents', driver="GPKG")
# Layer 2: closest_hospital
gdf_closest_hospital.to_file(OUTPUT_FILE, layer='closest_hospital', driver="GPKG")
# Layer 3: distance_matrix
# columns: incident_id, hospital_name, rank, network_distance_m
gdf_distance_matrix.to_file(OUTPUT_FILE, layer='distance_matrix', driver="GPKG")
# Layer 4: isochrones_15min
# columns: hospital_name, travel_time_min
gdf_isochrones.to_file(OUTPUT_FILE, layer='isochrones_15min', driver="GPKG")
print("Done!")
from shapely.ops import unary_union
if __name__ == "__main__":
solve()
Wrote /work/solve.py (11448 chars).
1:12
Bash
python solve.py
python solve.py
Downloading road network...
Traceback (most recent call last):
File "/work/solve.py", line 278, in <module>
solve()
~~~~~^^
File "/work/solve.py", line 34, in solve
G = ox.graph_from_bbox(north=north, south=south, east=east, west=west, network_type='drive')
TypeError: graph_from_bbox() got an unexpected keyword argument 'north'
1:16
Edit
/work/solve.py
old_string not found in /work/solve.py
1:20
Edit
/work/solve.py
old_string not found in /work/solve.py
1:21
Read
/work/solve.py
1 import pandas as pd
2 import geopandas as gpd
3 from shapely.geometry import Point, LineString, MultiPolygon, MultiPoint
4 import osmnx as ox
5 import networkx as nx
6 import numpy as np
7
8 # Configuration
9 BBOX = (48.83, 2.30, 48.88, 2.38) # (min_lat, min_lon, max_lat, max_lon)
10 OUTPUT_FILE = "emergency_routing.gpkg"
11 CRS = "EPSG:2154" # RGF93 / Lambert-93 (Official metric for France/Paris)
12
13 def solve():
14 # 1. Load incidents
15 df_incidents = pd.read_csv('incidents.csv')
16 # BBox order in osmnx/overpass is usually (north, south, east, west) or similar.
17 # Let's define the bounds properly.
18 # lat_min, lon_min, lat_max, lon_max = 48.83, 2.30, 48.88, 2.38
19
20 # Convert to GeoDataFrame
21 gdf_incidents = gpd.GeoDataFrame(
22 df_incidents,
23 geometry=gpd.points_from_xy(df_incidents.longitude, df_incidents.latitude),
24 crs="EPSG:4326"
25 ).to_crs(CRS)
26
27 # 2. Get Road Network and Hospitals using OSMnx
28 # BBox: north, south, east, west
29 north, south, east, west = 48.88, 48.83, 2.38, 2.30
30
31 print("Downloading road network...")
32 # We'll use the bbox to get the graph.
33 # osmnx.graph_from_bbox(north, south, east, west, network_type='drive')
34 G = ox.graph_from_bbox(north=north, south=south, east=east, west=west, network_type='drive')
35
36 # Apply speed rules
37 # OSMnx handles maxspeed if present.
38 # If maxspeed is missing, we set it to 30 km/h (8.333 m/s)
39 for u, v, k, data in G.edges(data=True, keys=True):
40 if 'maxspeed' not in data or data['maxspeed'] is None:
41 data['maxspeed'] = 30
42 elif isinstance(data['maxspeed'], str):
43 # Handle cases like "30 km/h"
44 try:
45 data['maxspeed'] = float(data['maxspeed'].split()[0])
46 except:
47 data['maxspeed'] = 30
48
49 # Calculate travel time (seconds)
50 # speed is in km/h, we need m/s
51 for u, v, k, data in G.edges(data=True, keys=True):
52 speed_mps = data['maxspeed'] / 3.6
53 data['travel_time'] = data['length'] / speed_mps
54
55 print("Downloading hospitals...")
56 # Get hospitals in the bbox
57 tags = {'amenity': 'hospital'}
58 # Since graph_from_bbox is used for network, let's use ox.features_from_bbox for hospitals
59 hospitals = ox.features_from_bbox(north, south, east, west, tags=tags)
60
61 # Filter for points or polygons representing hospitals
62 # We'll take the centroid for routing if they are polygons
63 hospitals = hospitals.copy()
64 hospitals['hospital_name'] = hospitals.get('name', 'Unknown Hospital')
65 # Clean names
66 hospitals['hospital_name'] = hospitals['hospital_name'].fillna('Unnamed Hospital')
67
68 # Convert to metric and extract points for routing
69 hospitals_gdf = hospitals.to_crs(CRS)
70 hospitals_gdf['geometry'] = hospitals_gdf.centroid # Use centroid for routing
71
72 # 3. Routing
73 print("Calculating routes and matrix...")
74
75 # Project G to CRS
76 G_proj = ox.project_graph(G, to_crs=CRS)
77
78 # We need to find the nearest nodes in the graph for incidents and hospitals
79 # Pre-calculate nearest nodes for hospitals
80 nodes, edges = ox.graph_to_gdfs(G_proj)
81
82 # Find nearest nodes for incidents and hospitals
83 # For incidents
84 inc_nodes = ox.distance.nearest_nodes(G_proj, gdf_incidents.geometry.x, gdf_incidents.geometry.y)
85 # For hospitals
86 hosp_nodes = ox.distance.nearest_nodes(G_proj, hospitals_gdf.geometry.x, hospitals_gdf.geometry.y)
87
88 # Map hospital index to node
89 hosp_info = []
90 for i, (idx, row) in enumerate(hospitals_gdf.iterrows()):
91 hosp_info.append({
92 'hospital_name': row['hospital_name'],
93 'node': hosp_nodes[i],
94 'original_index': i
95 })
96 hosp_df = pd.DataFrame(hosp_info)
97
98 # Lists to store results
99 closest_routes_data = []
100 distance_matrix_data = []
101
102 for i, inc_row in gdf_incidents.iterrows():
103 inc_id = inc_row['incident_id']
104 start_node = inc_nodes[i]
105
106 # Find distances to all hospitals
107 dist_to_hospitals = []
108 for _, h_row in hosp_df.iterrows():
109 end_node = h_row['node']
110 try:
111 # Shortest path distance (using 'length' edge attribute)
112 path_length = nx.shortest_path_length(G_proj, start_node, end_node, weight='length')
113 dist_to_hospitals.append({
114 'hospital_name': h_row['hospital_name'],
115 'node': end_node,
116 'network_distance_m': path_length
117 })
118 except nx.NetworkXNoPath:
119 continue
120
121 # Sort hospitals by distance
122 dist_to_hospitals.sort(key=lambda x: x['network_distance_m'])
123
124 # 1. Closest hospital route
125 if dist_to_hospitals:
126 closest = dist_to_hospitals[0]
127 path = nx.shortest_path(G_proj, start_node, closest['node'], weight='length')
128
129 # Create LineString for the route
130 route_coords = []
131 for node in path:
132 node_data = G_proj.nodes[node]
133 route_coords.append((node_data['x'], node_data['y']))
134
135 closest_routes_data.append({
136 'incident_id': inc_id,
137 'hospital_name': closest['hospital_name'],
138 'network_distance_m': closest['network_distance_m'],
139 'geometry': LineString(route_coords)
140 })
141
142 # 2. Distance Matrix (top 3)
143 for rank, h_dist in enumerate(dist_to_hospitals[:3], 1):
144 distance_matrix_data.append({
145 'incident_id': inc_id,
146 'hospital_name': h_dist['hospital_name'],
147 'rank': rank,
148 'network_distance_m': h_dist['network_distance_m'],
149 'geometry': None
150 })
151 else:
152 # Handle case with no path
153 distance_matrix_data.append({
154 'incident_id': inc_id,
155 'hospital_name': None,
156 'rank': None,
157 'network_distance_m': None,
158 'geometry': None
159 })
160
161 # Convert results to GDFs
162 gdf_closest_hospital = gpd.GeoDataFrame(closest_routes_data, crs=CRS)
163 gdf_distance_matrix = gpd.GeoDataFrame(distance_matrix_data, crs=CRS)
164 # Ensure distance_matrix has correct columns even if empty
165 for col in ['incident_id', 'hospital_name', 'rank', 'network_distance_m']:
166 if col not in gdf_distance_matrix.columns:
167 gdf_distance_matrix[col] = None
168
169 # 3. Isochrones (15 min = 900 seconds)
170 print("Calculating isochrones...")
171 isochrone_data = []
172 time_limit = 15 * 60 # 15 minutes in seconds
173
174 for _, h_row in hosp_df.iterrows():
175 h_node = h_row['node']
176 h_name = h_row['hospital_name']
177
178 # Get all nodes reachable within time_limit
179 subgraph = nx.ego_graph(G_proj, h_node, radius=time_limit, weight='travel_time')
180
181 if len(subgraph.nodes) > 0:
182 # Get the boundary (concave hull or union of node areas)
183 # A simpler way: union of all edges in the subgraph?
184 # Or just collect node coordinates and use alpha shape/convex hull
185 # For a "professional" look, let's use the edges to form a polygon.
186 # Or use the standard method: union of buffered nodes.
187
188 # Efficient way: get all nodes in subgraph, buffer them, then union.
189 # To avoid a massive union of many circles, we can collect nodes and use a concave hull.
190 node_points = []
191 for node in subgraph.nodes:
192 node_data = G_proj.nodes[node]
193 node_points.append(Point(node_data['x'], node_data['y']))
194
195 if node_points:
196 # Using a simple buffer approach for robustness
197 # Union of buffers of nodes
198 # To be faster and less messy, let's buffer the edges.
199 edge_geoms = []
200 for u, v, k, data in subgraph.edges(data=True, keys=True):
201 # Create line from edge
202 u_data = G_proj.nodes[u]
203 v_data = G_proj.nodes[v]
204 edge_geoms.append(LineString([(u_data['x'], u_data['y']), (v_data['x'], v_data['y'])]))
205
206 if edge_geoms:
207 # Buffer the edges slightly to create a "corridor"
208 # Then union them.
209 # Or just buffer the nodes and union.
210 # Let's try buffering nodes with a small radius (e.g. 20m) and unioning.
211 # Actually, a better way is to use the edges.
212 combined_edges = unary_union(edge_geoms)
213 poly = combined_edges.buffer(30) # 30m buffer around the reachable network
214
215 # We need to ensure it's a valid polygon/multipolygon
216 if not poly.is_empty:
217 isochrone_data.append({
218 'hospital_name': h_name,
219 'travel_time_min': 15.0,
220 'geometry': poly
221 })
222 else:
223 isochrone_data.append({
224 'hospital_name': h_name,
225 'travel_time_min': 15.0,
226 'geometry': None
227 })
228 else:
229 # Only the node itself is reachable
230 node_data = G_proj.nodes[h_node]
231 isochrone_data.append({
232 'hospital_name': h_name,
233 'travel_time_min': 15.0,
234 'geometry': Point(node_data['x'], node_data['y']).buffer(30)
235 })
236 else:
237 isochrone_data.append({
238 'hospital_name': h_name,
239 'travel_time_min': 15.0,
240 'geometry': None
241 })
242 else:
243 isochrone_data.append({
244 'hospital_name': h_name,
245 'travel_time_min': 15.0,
246 'geometry': None
247 })
248
249 gdf_isochrones = gpd.GeoDataFrame(isochrone_data, crs=CRS)
250
251 # 4. Write to GPKG
252 print(f"Writing to {OUTPUT_FILE}...")
253 # Remove rows with empty geometry for certain layers if needed,
254 # but requirements say "MultiPolygon geometry ... One row per hospital"
255 # We'll keep them but clean up.
256
257 # Save layers
258 # Requirement: incident_id, hospital_name, network_distance_m
259 # We have these in closest_routes_data
260 # Layer 1: incidents
261 gdf_incidents.to_file(OUTPUT_FILE, layer='incidents', driver="GPKG")
262
263 # Layer 2: closest_hospital
264 gdf_closest_hospital.to_file(OUTPUT_FILE, layer='closest_hospital', driver="GPKG")
265
266 # Layer 3: distance_matrix
267 # columns: incident_id, hospital_name, rank, network_distance_m
268 gdf_distance_matrix.to_file(OUTPUT_FILE, layer='distance_matrix', driver="GPKG")
269
270 # Layer 4: isochrones_15min
271 # columns: hospital_name, travel_time_min
272 gdf_isochrones.to_file(OUTPUT_FILE, layer='isochrones_15min', driver="GPKG")
273
274 print("Done!")
275
276 from shapely.ops import unary_union
277 if __name__ == "__main__":
278 solve()
1:24
Edit
/work/solve.py
Edited /work/solve.py.
1:26
Bash
python solve.py
python solve.py
Downloading road network...
Traceback (most recent call last):
File "/work/solve.py", line 278, in <module>
solve()
~~~~~^^
File "/work/solve.py", line 34, in solve
G = ox.graph_from_bbox(north, south, east, west, network_type='drive')
TypeError: graph_from_bbox() takes 1 positional argument but 4 positional arguments (and 1 keyword-only argument) were given
1:30
Edit
/work/solve.py
Edited /work/solve.py.
5:38
Write
/work/solve.py
content · 205 lines · py
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point, LineString, MultiPolygon, MultiPoint
import osmnx as ox
import networkx as nx
import numpy as np
from shapely.ops import unary_union
# Configuration
OUTPUT_FILE = "emergency_routing.gpkg"
CRS = "EPSG:2154" # RGF93 / Lambert-93 (Official metric for France/Paris)
def solve():
# 1. Load incidents
df_incidents = pd.read_csv('incidents.csv')
# Convert to GeoDataFrame
gdf_incidents = gpd.GeoDataFrame(
df_incidents,
geometry=gpd.points_from_xy(df_incidents.longitude, df_incidents.latitude),
crs="EPSG:4326"
).to_crs(CRS)
# 2. Get Road Network and Hospitals using OSMnx
# BBox: north, south, east, west
north, south, east, west = 48.88, 48.83, 2.38, 2.30
print("Downloading road network...")
G = ox.graph_from_bbox(north=north, south=south, east=east, west=west, network_type='drive')
# Apply speed rules
for u, v, k, data in G.edges(data=True, keys=True):
if 'maxspeed' not in data or data['maxspeed'] is None:
data['maxspeed'] = 30
elif isinstance(data['maxspeed'], str):
try:
# Extracting number from strings like '30 km/h'
val = data['maxspeed'].split()[0]
data['maxspeed'] = float(val)
except:
data['maxspeed'] = 30
else:
# If it's already numeric
pass
# Calculate travel time (seconds)
for u, v, k, data in G.edges(data=True, keys=True):
speed_mps = data['maxspeed'] / 3.6
data['travel_time'] = data['length'] / speed_mps
print("Downloading hospitals...")
tags = {'amenity': 'hospital'}
hospitals = ox.features_from_bbox(north=north, south=south, east=east, west=west, tags=tags)
hospitals = hospitals.copy()
hospitals['hospital_name'] = hospitals.get('name', 'Unknown Hospital')
hospitals['hospital_name'] = hospitals['hospital_name'].fillna('Unnamed Hospital')
# Convert to metric and extract points for routing
hospitals_gdf = hospitals.to_crs(CRS)
# Using centroid for routing
hospitals_gdf['geometry'] = hospitals_gdf.centroid
# 3. Routing
print("Calculating routes and matrix...")
G_proj = ox.project_graph(G, to_crs=CRS)
# Find nearest nodes for incidents and hospitals
inc_nodes = ox.distance.nearest_nodes(G_proj, gdf_incidents.geometry.x, gdf_incidents.geometry.y)
hosp_nodes = ox.distance.nearest_nodes(G_proj, hospitals_gdf.geometry.x, hospitals_gdf.geometry.y)
# Map hospital index to node
hosp_info = []
for i, (idx, row) in enumerate(hospitals_gdf.iterrows()):
hosp_info.append({
'hospital_name': row['hospital_name'],
'node': hosp_nodes.iloc[i],
'original_index': i
})
hosp_df = pd.DataFrame(hosp_info)
closest_routes_data = []
distance_matrix_data = []
for i, inc_row in gdf_incidents.iterrows():
inc_id = inc_row['incident_id']
start_node = inc_nodes.iloc[i]
dist_to_hospitals = []
for _, h_row in hosp_df.iterrows():
end_node = h_row['node']
try:
# Using length weight for distance
path_length = nx.shortest_path_length(G_proj, start_node, end_node, weight='length')
dist_to_hospitals.append({
'hospital_name': h_row['hospital_name'],
'node': end_node,
'network_distance_m': path_length
})
except (nx.NetworkXNoPath, nx.NodeNotFound):
continue
dist_to_hospitals.sort(key=lambda x: x['network_distance_m'])
# 1. Closest hospital route
if dist_to_hospitals:
closest = dist_to_hospitals[0]
try:
path = nx.shortest_path(G_proj, start_node, closest['node'], weight='length')
route_coords = []
for node in path:
node_data = G_proj.nodes[node]
route_coords.append((node_data['x'], node_data['y']))
closest_routes_data.append({
'incident_id': inc_id,
'hospital_name': closest['hospital_name'],
'network_distance_m': float(closest['network_distance_m']),
'geometry': LineString(route_coords)
})
except:
pass
# 2. Distance Matrix (top 3)
for rank, h_dist in enumerate(dist_to_hospitals[:3], 1):
distance_matrix_data.append({
'incident_id': inc_id,
'hospital_name': h_dist['hospital_name'],
'rank': rank,
'network_distance_m': float(h_dist['network_distance_m']),
'geometry': None
})
else:
# If no hospital reachable, still need a placeholder for the matrix logic?
# The prompt implies we'll find routes.
pass
gdf_closest_hospital = gpd.GeoDataFrame(closest_routes_data, crs=CRS)
gdf_distance_matrix = gpd.GeoDataFrame(distance_matrix_data, crs=CRS)
# 3. Isochrones (15 min = 900 seconds)
print("Calculating isochrones...")
isochrone_data = []
time_limit = 15 * 60
for _, h_row in hosp_df.iterrows():
h_node = h_row['node']
h_name = h_row['hospital_name']
try:
# ego_graph with travel_time weight
subgraph = nx.ego_graph(G_proj, h_node, radius=time_limit, weight='travel_time')
if len(subgraph.nodes) > 0:
edge_geoms = []
for u, v, k, data in subgraph.edges(data=True, keys=True):
u_data = G_proj.nodes[u]
v_data = G_proj.nodes[v]
edge_geoms.append(LineString([(u_data['x'], u_data['y']), (v_data['x'], v_data['y'])]))
if edge_geoms:
combined_edges = unary_union(edge_geoms)
poly = combined_edges.buffer(30)
if not poly.is_empty:
isochrone_data.append({
'hospital_name': h_name,
'travel_time_min': 15.0,
'geometry': poly
})
else:
node_data = G_proj.nodes[h_node]
isochrone_data.append({
'hospital_name': h_name,
'travel_time_min': 15.0,
'geometry': Point(node_data['x'], node_data['y']).buffer(30)
})
else:
# fallback if no nodes found
node_data = G_proj.nodes[h_node]
isochrone_data.append({
'hospital_name': h_name,
'travel_time_min': 15.0,
'geometry': Point(node_data['x'], node_data['y']).buffer(30)
})
except:
isochrone_data.append({
'hospital_name': h_name,
'travel_time_min': 15.0,
'geometry': None
})
gdf_isochrones = gpd.GeoDataFrame(isochrone_data, crs=CRS)
# 4. Write to GPKG
print(f"Writing to {OUTPUT_FILE}...")
gdf_incidents.to_file(OUTPUT_FILE, layer='incidents', driver="GPKG")
gdf_closest_hospital.to_file(OUTPUT_FILE, layer='closest_hospital', driver="GPKG")
gdf_distance_matrix.to_file(OUTPUT_FILE, layer='distance_matrix', driver="GPKG")
gdf_isochrones.to_file(OUTPUT_FILE, layer='isochrones_15min', driver="GPKG")
print("Done!")
if __name__ == "__main__":
solve()
Wrote /work/solve.py (7911 chars).
5:39
Bash
python solve.py
python solve.py
Downloading road network...
Traceback (most recent call last):
File "/work/solve.py", line 204, in <module>
solve()
~~~~~^^
File "/work/solve.py", line 29, in solve
G = ox.graph_from_bbox(north=north, south=south, east=east, west=west, network_type='drive')
TypeError: graph_from_bbox() got an unexpected keyword argument 'north'
5:42
Bash
python -c "import osmnx; print(osmnx.__version__)"
python -c "import osmnx; print(osmnx.__version__)"
2.1.0
5:45
Bash
python -c "import osmnx; help(osmnx.graph_from_bbox)"
python -c "import osmnx; help(osmnx.graph_from_bbox)"
Help on function graph_from_bbox in module osmnx.graph:
graph_from_bbox(
bbox: tuple[float, float, float, float],
*,
network_type: str = 'all',
simplify: bool = True,
retain_all: bool = False,
truncate_by_edge: bool = False,
custom_filter: str | list[str] | None = None
) -> nx.MultiDiGraph
Download and create a graph within a lat-lon bounding box.
This function uses filters to query the Overpass API: you can either
specify a pre-defined `network_type` or provide your own `custom_filter`
with Overpass QL.
Use the `settings` module's `useful_tags_node` and `useful_tags_way`
settings to configure which OSM node/way tags are added as graph node/edge
attributes. If you want a fully bidirectional network, ensure your
`network_type` is in `settings.bidirectional_network_types` before
creating your graph. You can also use the `settings` module to retrieve a
snapshot of historical OSM data as of a certain date, or to configure the
Overpass server timeout, memory allocation, and other customizations.
Parameters
----------
bbox
Bounding box as `(left, bottom, right, top)`. Coordinates should be in
unprojected latitude-longitude degrees (EPSG:4326).
network_type
{"all", "all_public", "bike", "drive", "drive_service", "walk"}
What type of street network to retrieve if `custom_filter` is None.
simplify
If True, simplify graph topology via the `simplify_graph` function.
retain_all
If True, return the entire graph even if it is not connected. If
False, retain only the largest weakly connected component.
truncate_by_edge
If True, retain nodes the outside bounding box if at least one of
the node's neighbors lies within the bounding box.
custom_filter
A custom ways filter to be used instead of the `network_type` presets,
e.g. `'["power"~"line"]' or '["highway"~"motorway|trunk"]'`. If `str`,
the intersection of keys/values will be used, e.g., `'[maxspeed=50][lanes=2]'`
will return all ways having both maxspeed of 50 and two lanes. If
`list`, the union of the `list` items will be used, e.g.,
`['[maxspeed=50]', '[lanes=2]']` will return all ways having either
maximum speed of 50 or two lanes. Also pass in a `network_type` that
is in `settings.bidirectional_network_types` if you want the graph to
be fully bidirectional.
Returns
-------
G
The resulting MultiDiGraph.
Notes
-----
Very large query areas use the `utils_geo._consolidate_subdivide_geometry`
function to automatically make multiple requests: see that function's
documentation for caveats.
5:51
Edit
/work/solve.py
Edited /work/solve.py.
5:55
Edit
/work/solve.py
Edited /work/solve.py.
5:56
Bash
python solve.py
python solve.py
Downloading road network...
Traceback (most recent call last):
File "/work/solve.py", line 205, in <module>
solve()
~~~~~^^
File "/work/solve.py", line 49, in solve
speed_mps = data['maxspeed'] / 3.6
~~~~~~~~~~~~~~~~~^~~~~
TypeError: unsupported operand type(s) for /: 'list' and 'float'
6:13
Edit
/work/solve.py
Edited /work/solve.py.
6:14
Bash
python solve.py
python solve.py
Downloading road network... Downloading hospitals... Calculating routes and matrix...
Traceback (most recent call last):
File "/work/solve.py", line 220, in <module>
solve()
~~~~~^^
File "/work/solve.py", line 85, in solve
inc_nodes = ox.distance.nearest_nodes(G_proj, gdf_incidents.geometry.x, gdf_incidents.geometry.y)
File "/usr/local/lib/python3.14/site-packages/osmnx/distance.py", line 368, in nearest_nodes
raise ImportError(msg)
ImportError: scipy must be installed as an optional dependency to search a projected graph.
6:22
Bash
pip install scipy
pip install scipy
Defaulting to user installation because normal site-packages is not writeable Collecting scipy Downloading scipy-1.17.1-cp314-cp314-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl.metadata (62 kB) Requirement already satisfied: numpy<2.7,>=1.26.4 in /usr/local/lib/python3.14/site-packages (from scipy) (2.4.6) Downloading scipy-1.17.1-cp314-cp314-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl (35.2 MB) ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 35.2/35.2 MB 62.3 MB/s 0:00:00 Installing collected packages: scipy Successfully installed scipy-1.17.1
6:31
Bash
python solve.py
python solve.py
Downloading road network... Downloading hospitals... Calculating routes and matrix...
Traceback (most recent call last):
File "/work/solve.py", line 220, in <module>
solve()
~~~~~^^
File "/work/solve.py", line 93, in solve
'node': hosp_nodes.iloc[i],
^^^^^^^^^^^^^^^
AttributeError: 'numpy.ndarray' object has no attribute 'iloc'
6:44
Edit
/work/solve.py
Edited /work/solve.py.
6:46
Bash
python solve.py
python solve.py
Downloading road network... Downloading hospitals... Calculating routes and matrix...
Traceback (most recent call last):
File "/work/solve.py", line 220, in <module>
solve()
~~~~~^^
File "/work/solve.py", line 103, in solve
start_node = inc_nodes.iloc[i]
^^^^^^^^^^^^^^
AttributeError: 'numpy.ndarray' object has no attribute 'iloc'
7:31
Edit
/work/solve.py
old_string not found in /work/solve.py
7:34
Read
/work/solve.py
1 import pandas as pd
2 import geopandas as gpd
3 from shapely.geometry import Point, LineString, MultiPolygon, MultiPoint
4 import osmnx as ox
5 import networkx as nx
6 import numpy as np
7 from shapely.ops import unary_union
8
9 # Configuration
10 OUTPUT_FILE = "emergency_routing.gpkg"
11 CRS = "EPSG:2154" # RGF93 / Lambert-93 (Official metric for France/Paris)
12
13 def solve():
14 # 1. Load incidents
15 df_incidents = pd.read_csv('incidents.csv')
16
17 # Convert to GeoDataFrame
18 gdf_incidents = gpd.GeoDataFrame(
19 df_incidents,
20 geometry=gpd.points_from_xy(df_incidents.longitude, df_incidents.latitude),
21 crs="EPSG:4326"
22 ).to_crs(CRS)
23
24 # 2. Get Road Network and Hospitals using OSMnx
25 # BBox: north, south, east, west
26 north, south, east, west = 48.88, 48.83, 2.38, 2.30
27
28 print("Downloading road network...")
29 # Bbox as (left, bottom, right, top) i.e., (west, south, east, north)
30 G = ox.graph_from_bbox(bbox=(west, south, east, north), network_type='drive')
31
32 # Apply speed rules
33 for u, v, k, data in G.edges(data=True, keys=True):
34 if 'maxspeed' not in data or data['maxspeed'] is None:
35 data['maxspeed'] = 30
36 elif isinstance(data['maxspeed'], str):
37 try:
38 # Extracting number from strings like '30 km/h'
39 val = data['maxspeed'].split()[0]
40 data['maxspeed'] = float(val)
41 except:
42 data['maxspeed'] = 30
43 else:
44 # If it's already numeric
45 pass
46
47 # Calculate travel time (seconds)
48 for u, v, k, data in G.edges(data=True, keys=True):
49 try:
50 speed = data['maxspeed']
51 if isinstance(speed, list):
52 speed = speed[0]
53 if isinstance(speed, str):
54 try:
55 speed = float(speed.split()[0])
56 except:
57 speed = 30
58
59 if not isinstance(speed, (int, float)):
60 speed = 30
61
62 speed_mps = speed / 3.6
63 data['travel_time'] = data['length'] / speed_mps
64 except:
65 data['travel_time'] = data['length'] / (30 / 3.6)
66
67 print("Downloading hospitals...")
68 tags = {'amenity': 'hospital'}
69 hospitals = ox.features_from_bbox(bbox=(west, south, east, north), tags=tags)
70
71 hospitals = hospitals.copy()
72 hospitals['hospital_name'] = hospitals.get('name', 'Unknown Hospital')
73 hospitals['hospital_name'] = hospitals['hospital_name'].fillna('Unnamed Hospital')
74
75 # Convert to metric and extract points for routing
76 hospitals_gdf = hospitals.to_crs(CRS)
77 # Using centroid for routing
78 hospitals_gdf['geometry'] = hospitals_gdf.centroid
79
80 # 3. Routing
81 print("Calculating routes and matrix...")
82 G_proj = ox.project_graph(G, to_crs=CRS)
83
84 # Find nearest nodes for incidents and hospitals
85 inc_nodes = ox.distance.nearest_nodes(G_proj, gdf_incidents.geometry.x, gdf_incidents.geometry.y)
86 hosp_nodes = ox.distance.nearest_nodes(G_proj, hospitals_gdf.geometry.x, hospitals_gdf.geometry.y)
87
88 # Map hospital index to node
89 hosp_info = []
90 for i, (idx, row) in enumerate(hospitals_gdf.iterrows()):
91 hosp_info.append({
92 'hospital_name': row['hospital_name'],
93 'node': hosp_nodes[i],
94 'original_index': i
95 })
96 hosp_df = pd.DataFrame(hosp_info)
97
98 closest_routes_data = []
99 distance_matrix_data = []
100
101 for i, inc_row in gdf_incidents.iterrows():
102 inc_id = inc_row['incident_id']
103 start_node = inc_nodes.iloc[i]
104
105 dist_to_hospitals = []
106 for _, h_row in hosp_df.iterrows():
107 end_node = h_row['node']
108 try:
109 # Using length weight for distance
110 path_length = nx.shortest_path_length(G_proj, start_node, end_node, weight='length')
111 dist_to_hospitals.append({
112 'hospital_name': h_row['hospital_name'],
113 'node': end_node,
114 'network_distance_m': path_length
115 })
116 except (nx.NetworkXNoPath, nx.NodeNotFound):
117 continue
118
119 dist_to_hospitals.sort(key=lambda x: x['network_distance_m'])
120
121 # 1. Closest hospital route
122 if dist_to_hospitals:
123 closest = dist_to_hospitals[0]
124 try:
125 path = nx.shortest_path(G_proj, start_node, closest['node'], weight='length')
126 route_coords = []
127 for node in path:
128 node_data = G_proj.nodes[node]
129 route_coords.append((node_data['x'], node_data['y']))
130
131 closest_routes_data.append({
132 'incident_id': inc_id,
133 'hospital_name': closest['hospital_name'],
134 'network_distance_m': float(closest['network_distance_m']),
135 'geometry': LineString(route_coords)
136 })
137 except:
138 pass
139
140 # 2. Distance Matrix (top 3)
141 for rank, h_dist in enumerate(dist_to_hospitals[:3], 1):
142 distance_matrix_data.append({
143 'incident_id': inc_id,
144 'hospital_name': h_dist['hospital_name'],
145 'rank': rank,
146 'network_distance_m': float(h_dist['network_distance_m']),
147 'geometry': None
148 })
149 else:
150 # If no hospital reachable, still need a placeholder for the matrix logic?
151 # The prompt implies we'll find routes.
152 pass
153
154 gdf_closest_hospital = gpd.GeoDataFrame(closest_routes_data, crs=CRS)
155 gdf_distance_matrix = gpd.GeoDataFrame(distance_matrix_data, crs=CRS)
156
157 # 3. Isochrones (15 min = 900 seconds)
158 print("Calculating isochrones...")
159 isochrone_data = []
160 time_limit = 15 * 60
161
162 for _, h_row in hosp_df.iterrows():
163 h_node = h_row['node']
164 h_name = h_row['hospital_name']
165
166 try:
167 # ego_graph with travel_time weight
168 subgraph = nx.ego_graph(G_proj, h_node, radius=time_limit, weight='travel_time')
169
170 if len(subgraph.nodes) > 0:
171 edge_geoms = []
172 for u, v, k, data in subgraph.edges(data=True, keys=True):
173 u_data = G_proj.nodes[u]
174 v_data = G_proj.nodes[v]
175 edge_geoms.append(LineString([(u_data['x'], u_data['y']), (v_data['x'], v_data['y'])]))
176
177 if edge_geoms:
178 combined_edges = unary_union(edge_geoms)
179 poly = combined_edges.buffer(30)
180 if not poly.is_empty:
181 isochrone_data.append({
182 'hospital_name': h_name,
183 'travel_time_min': 15.0,
184 'geometry': poly
185 })
186 else:
187 node_data = G_proj.nodes[h_node]
188 isochrone_data.append({
189 'hospital_name': h_name,
190 'travel_time_min': 15.0,
191 'geometry': Point(node_data['x'], node_data['y']).buffer(30)
192 })
193 else:
194 # fallback if no nodes found
195 node_data = G_proj.nodes[h_node]
196 isochrone_data.append({
197 'hospital_name': h_name,
198 'travel_time_min': 15.0,
199 'geometry': Point(node_data['x'], node_data['y']).buffer(30)
200 })
201 except:
202 isochrone_data.append({
203 'hospital_name': h_name,
204 'travel_time_min': 15.0,
205 'geometry': None
206 })
207
208 gdf_isochrones = gpd.GeoDataFrame(isochrone_data, crs=CRS)
209
210 # 4. Write to GPKG
211 print(f"Writing to {OUTPUT_FILE}...")
212 gdf_incidents.to_file(OUTPUT_FILE, layer='incidents', driver="GPKG")
213 gdf_closest_hospital.to_file(OUTPUT_FILE, layer='closest_hospital', driver="GPKG")
214 gdf_distance_matrix.to_file(OUTPUT_FILE, layer='distance_matrix', driver="GPKG")
215 gdf_isochrones.to_file(OUTPUT_FILE, layer='isochrones_15min', driver="GPKG")
216
217 print("Done!")
218
219 if __name__ == "__main__":
220 solve()
7:36
Bash
python solve.py
python solve.py
Downloading road network... Downloading hospitals... Calculating routes and matrix...
Traceback (most recent call last):
File "/work/solve.py", line 220, in <module>
solve()
~~~~~^^
File "/work/solve.py", line 103, in solve
start_node = inc_nodes.iloc[i]
^^^^^^^^^^^^^^
AttributeError: 'numpy.ndarray' object has no attribute 'iloc'
8:10
Edit
/work/solve.py
Edited /work/solve.py.
8:11
Bash
python solve.py
python solve.py
Downloading road network... Downloading hospitals... Calculating routes and matrix... Calculating isochrones... Writing to emergency_routing.gpkg... Done!
8:34