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:00
Bash
ls -F
incidents.csv
0:01
Read
incidents.csv
1 incident_id,latitude,longitude,label
2 INC_001,48.8566,2.3522,near Notre-Dame
3 INC_002,48.8620,2.3360,near Louvre
4 INC_003,48.8462,2.3464,Quartier Latin
5 INC_004,48.8700,2.3431,near Gare du Nord area
6 INC_005,48.8530,2.3325,Luxembourg Garden area
7 INC_006,48.8600,2.3700,Bastille area
8 INC_007,48.8450,2.3100,near Montparnasse
9 INC_008,48.8750,2.3600,near Republique
0:38
Write
solve.py
content · 197 lines · py
import osmnx as ox
import networkx as nx
import geopandas as gpd
import pandas as pd
from shapely.geometry import Point, LineString, MultiLineString, Polygon, MultiPolygon
import numpy as np
# Constants
BBOX = (48.83, 2.30, 48.88, 2.38) # min_lat, min_lon, max_lat, max_lon
CRS = "EPSG:2154" # RGF93 / Lambert-93 (Official metric CRS for France/Paris)
SPEED_DEFAULT = 30 / 3.6 # 30 km/h to m/s
def solve():
# 1. Load Incidents
incidents_df = pd.read_csv('incidents.csv')
# Create GeoDataFrame
incidents_gdf = gpd.GeoDataFrame(
incidents_df,
geometry=gpd.points_from_xy(incidents_df.longitude, incidents_df.latitude),
crs="EPSG:4326"
).to_crs(CRS)
# 2. Fetch OSM Data (Network and Hospitals)
print("Fetching OSM network and hospitals...")
# Use osmnx to get the graph
# bbox: north, south, east, west
north, south, east, west = BBOX[2], BBOX[0], BBOX[3], BBOX[1]
G = ox.graph_from_bbox(north, south, east, west, network_type='drive')
# Apply speeds
G = ox.add_edge_speeds(G, fallback=30)
G = ox.add_edge_travel_times(G)
# Fetch Hospitals using Overpass (via osmnx or direct query if needed)
# osmnx doesn't have a direct 'hospitals' function but we can use query
tags = {"amenity": "hospital"}
hospitals_gdf = ox.features_from_bbox(north, south, east, west, tags=tags)
# Clean hospitals: ensure they are points/polygons and have names
# Some might be Multipolygons or Points
hospitals_gdf = hospitals_gdf.copy()
if 'name' not in hospitals_gdf.columns:
hospitals_gdf['name'] = 'Unknown Hospital'
hospitals_gdf['hospital_name'] = hospitals_gdf['name'].fillna('Unknown Hospital')
# We want the centroid for routing if they are polygons
hospitals_gdf['geometry'] = hospitals_gdf.geometry.centroid
hospitals_gdf = hospitals_gdf.to_crs(CRS)
# 3. Routing Analysis
print("Calculating routes...")
# Pre-calculate nearest nodes for incidents and hospitals
inc_nodes = ox.distance.nearest_nodes(G, incidents_gdf.geometry.x, incidents_gdf.geometry.y)
hosp_nodes = ox.distance.nearest_nodes(G, hospitals_gdf.geometry.x, hospitals_gdf.geometry.y)
closest_routes_data = []
distance_matrix_data = []
for idx, inc in incidents_gdf.iterrows():
inc_id = inc['incident_id']
u = inc_nodes[idx]
# Calculate distances to all hospitals for this incident
hosp_distances = []
for h_idx, hosp in hospitals_gdf.iterrows():
v = hosp_nodes[h_idx]
try:
# weight='length' to get network distance
d = nx.shortest_path_length(G, u, v, weight='length')
hosp_distances.append({
'hospital_name': hosp['hospital_name'],
'network_distance_m': d,
'hosp_idx': h_idx,
'v_node': v
})
except nx.NetworkXNoPath:
continue
# Sort by distance
hosp_distances.sort(key=lambda x: x['network_distance_m'])
# Closest hospital for 'closest_hospital' layer
if hosp_distances:
closest = hosp_distances[0]
# Get path for the line geometry
path = nx.shortest_path(G, u, closest['v_node'], weight='length')
# Convert path to LineString
nodes_coords = []
for node in path:
node_data = G.nodes[node]
# Get lat/lon from graph node to ensure accurate geometry
# Actually G.nodes[node] contains 'x' and 'y' (projected)
nodes_coords.append((G.nodes[node]['x'], G.nodes[node]['y']))
route_line = LineString(nodes_coords)
closest_routes_data.append({
'incident_id': inc_id,
'hospital_name': closest['hospital_name'],
'network_distance_m': closest['network_distance_m'],
'geometry': route_line
})
# Distance matrix (top 3)
for rank, h_info in enumerate(hosp_distances[:3], 1):
distance_matrix_data.append({
'incident_id': inc_id,
'hospital_name': h_info['hospital_name'],
'rank': rank,
'network_distance_m': h_info['network_distance_m']
})
else:
# Fallback if no hospital reachable
distance_matrix_data.append({
'incident_id': inc_id,
'hospital_name': 'None',
'rank': 1,
'network_distance_m': np.nan
})
# Convert to GeoDataFrames
closest_hospital_gdf = gpd.GeoDataFrame(closest_routes_data, crs=CRS)
distance_matrix_gdf = gpd.GeoDataFrame(distance_matrix_data, crs=CRS)
# distance_matrix_gdf.geometry = None # explicitly empty
# 4. Isochrones (15 min)
print("Calculating isochrones...")
isochrone_data = []
# 15 minutes = 900 seconds
time_limit = 15 * 60
for idx, hosp in hospitals_gdf.iterrows():
h_name = hosp['hospital_name']
v = hosp_nodes[idx]
# Subgraph of reachable nodes within time limit
# Use Dijkstra to find all nodes within time_limit
subgraph_nodes = nx.single_source_dijkstra_path_length(G, v, cutoff=time_limit, weight='travel_time')
if subgraph_nodes:
# Get geometries of these nodes
node_points = [Point(G.nodes[n]['x'], G.nodes[n]['y']) for n in subgraph_nodes.keys()]
# Create a concave hull or buffer?
# Usually, isochrones are polygons around the nodes.
# A simple approach is to buffer the nodes and take the union,
# or better, use the edges.
# Let's use the nodes and create a buffer/concave hull.
# For speed and robustness, we'll buffer the points and dissolve.
points_gdf = gpd.GeoDataFrame(geometry=node_points, crs=CRS)
# Buffer points by a small amount to ensure connectivity in union
isochrone_poly = points_gdf.buffer(20).union_all()
# We use union_all() for newer geopandas, or unary_union for older.
# If union_all is not available:
# try: isochrone_poly = points_gdf.buffer(20).unary_union
# except: isochrone_poly = points_gdf.buffer(20).union_all()
isochrone_data.append({
'hospital_name': h_name,
'travel_time_min': 15.0,
'geometry': isochrone_poly
})
else:
# If no nodes reachable, we still need a row
isochrone_data.append({
'hospital_name': h_name,
'travel_time_min': 0.0,
'geometry': Point(hosp.geometry.x, hosp.geometry.y).buffer(0)
})
isochrones_gdf = gpd.GeoDataFrame(isochrone_data, crs=CRS)
# 5. Save to GPKG
print("Saving to GPKG...")
# Prepare incidents layer
incidents_output = incidents_gdf[['incident_id', 'geometry']]
with gpd.GeoDataFrame(columns=incidents_output.columns, crs=CRS) as empty_gdf:
# We can't use 'with' like that for writing. Use the engine.
pass
# Write layers
# Layer 1: incidents
incidents_output.to_file('emergency_routing.gpkg', layer='incidents', driver='GPKG')
# Layer 2: closest_hospital
closest_hospital_gdf.to_file('emergency_routing.gpkg', layer='closest_hospital', driver='GPKG')
# Layer 3: distance_matrix
distance_matrix_gdf.to_file('emergency_routing.gpkg', layer='distance_matrix', driver='GPKG')
# Layer 4: isochrones_15min
isochrones_gdf.to_file('emergency_routing.gpkg', layer='isochrones_15min', driver='GPKG')
print("Done!")
if __name__ == "__main__":
solve()
Wrote solve.py (7969 chars).
0:41
Edit
solve.py
old_string not found in solve.py
0:45
Edit
solve.py
old_string not found in solve.py
0:46
Read
solve.py
1 import osmnx as ox
2 import networkx as nx
3 import geopandas as gpd
4 import pandas as pd
5 from shapely.geometry import Point, LineString, MultiLineString, Polygon, MultiPolygon
6 import numpy as np
7
8 # Constants
9 BBOX = (48.83, 2.30, 48.88, 2.38) # min_lat, min_lon, max_lat, max_lon
10 CRS = "EPSG:2154" # RGF93 / Lambert-93 (Official metric CRS for France/Paris)
11 SPEED_DEFAULT = 30 / 3.6 # 30 km/h to m/s
12
13 def solve():
14 # 1. Load Incidents
15 incidents_df = pd.read_csv('incidents.csv')
16 # Create GeoDataFrame
17 incidents_gdf = gpd.GeoDataFrame(
18 incidents_df,
19 geometry=gpd.points_from_xy(incidents_df.longitude, incidents_df.latitude),
20 crs="EPSG:4326"
21 ).to_crs(CRS)
22
23 # 2. Fetch OSM Data (Network and Hospitals)
24 print("Fetching OSM network and hospitals...")
25 # Use osmnx to get the graph
26 # bbox: north, south, east, west
27 north, south, east, west = BBOX[2], BBOX[0], BBOX[3], BBOX[1]
28 G = ox.graph_from_bbox(north, south, east, west, network_type='drive')
29
30 # Apply speeds
31 G = ox.add_edge_speeds(G, fallback=30)
32 G = ox.add_edge_travel_times(G)
33
34 # Fetch Hospitals using Overpass (via osmnx or direct query if needed)
35 # osmnx doesn't have a direct 'hospitals' function but we can use query
36 tags = {"amenity": "hospital"}
37 hospitals_gdf = ox.features_from_bbox(north, south, east, west, tags=tags)
38
39 # Clean hospitals: ensure they are points/polygons and have names
40 # Some might be Multipolygons or Points
41 hospitals_gdf = hospitals_gdf.copy()
42 if 'name' not in hospitals_gdf.columns:
43 hospitals_gdf['name'] = 'Unknown Hospital'
44 hospitals_gdf['hospital_name'] = hospitals_gdf['name'].fillna('Unknown Hospital')
45
46 # We want the centroid for routing if they are polygons
47 hospitals_gdf['geometry'] = hospitals_gdf.geometry.centroid
48 hospitals_gdf = hospitals_gdf.to_crs(CRS)
49
50 # 3. Routing Analysis
51 print("Calculating routes...")
52
53 # Pre-calculate nearest nodes for incidents and hospitals
54 inc_nodes = ox.distance.nearest_nodes(G, incidents_gdf.geometry.x, incidents_gdf.geometry.y)
55 hosp_nodes = ox.distance.nearest_nodes(G, hospitals_gdf.geometry.x, hospitals_gdf.geometry.y)
56
57 closest_routes_data = []
58 distance_matrix_data = []
59
60 for idx, inc in incidents_gdf.iterrows():
61 inc_id = inc['incident_id']
62 u = inc_nodes[idx]
63
64 # Calculate distances to all hospitals for this incident
65 hosp_distances = []
66 for h_idx, hosp in hospitals_gdf.iterrows():
67 v = hosp_nodes[h_idx]
68 try:
69 # weight='length' to get network distance
70 d = nx.shortest_path_length(G, u, v, weight='length')
71 hosp_distances.append({
72 'hospital_name': hosp['hospital_name'],
73 'network_distance_m': d,
74 'hosp_idx': h_idx,
75 'v_node': v
76 })
77 except nx.NetworkXNoPath:
78 continue
79
80 # Sort by distance
81 hosp_distances.sort(key=lambda x: x['network_distance_m'])
82
83 # Closest hospital for 'closest_hospital' layer
84 if hosp_distances:
85 closest = hosp_distances[0]
86 # Get path for the line geometry
87 path = nx.shortest_path(G, u, closest['v_node'], weight='length')
88 # Convert path to LineString
89 nodes_coords = []
90 for node in path:
91 node_data = G.nodes[node]
92 # Get lat/lon from graph node to ensure accurate geometry
93 # Actually G.nodes[node] contains 'x' and 'y' (projected)
94 nodes_coords.append((G.nodes[node]['x'], G.nodes[node]['y']))
95
96 route_line = LineString(nodes_coords)
97
98 closest_routes_data.append({
99 'incident_id': inc_id,
100 'hospital_name': closest['hospital_name'],
101 'network_distance_m': closest['network_distance_m'],
102 'geometry': route_line
103 })
104
105 # Distance matrix (top 3)
106 for rank, h_info in enumerate(hosp_distances[:3], 1):
107 distance_matrix_data.append({
108 'incident_id': inc_id,
109 'hospital_name': h_info['hospital_name'],
110 'rank': rank,
111 'network_distance_m': h_info['network_distance_m']
112 })
113 else:
114 # Fallback if no hospital reachable
115 distance_matrix_data.append({
116 'incident_id': inc_id,
117 'hospital_name': 'None',
118 'rank': 1,
119 'network_distance_m': np.nan
120 })
121
122 # Convert to GeoDataFrames
123 closest_hospital_gdf = gpd.GeoDataFrame(closest_routes_data, crs=CRS)
124 distance_matrix_gdf = gpd.GeoDataFrame(distance_matrix_data, crs=CRS)
125 # distance_matrix_gdf.geometry = None # explicitly empty
126
127 # 4. Isochrones (15 min)
128 print("Calculating isochrones...")
129 isochrone_data = []
130 # 15 minutes = 900 seconds
131 time_limit = 15 * 60
132
133 for idx, hosp in hospitals_gdf.iterrows():
134 h_name = hosp['hospital_name']
135 v = hosp_nodes[idx]
136
137 # Subgraph of reachable nodes within time limit
138 # Use Dijkstra to find all nodes within time_limit
139 subgraph_nodes = nx.single_source_dijkstra_path_length(G, v, cutoff=time_limit, weight='travel_time')
140
141 if subgraph_nodes:
142 # Get geometries of these nodes
143 node_points = [Point(G.nodes[n]['x'], G.nodes[n]['y']) for n in subgraph_nodes.keys()]
144 # Create a concave hull or buffer?
145 # Usually, isochrones are polygons around the nodes.
146 # A simple approach is to buffer the nodes and take the union,
147 # or better, use the edges.
148 # Let's use the nodes and create a buffer/concave hull.
149 # For speed and robustness, we'll buffer the points and dissolve.
150 points_gdf = gpd.GeoDataFrame(geometry=node_points, crs=CRS)
151 # Buffer points by a small amount to ensure connectivity in union
152 isochrone_poly = points_gdf.buffer(20).union_all()
153
154 # We use union_all() for newer geopandas, or unary_union for older.
155 # If union_all is not available:
156 # try: isochrone_poly = points_gdf.buffer(20).unary_union
157 # except: isochrone_poly = points_gdf.buffer(20).union_all()
158
159 isochrone_data.append({
160 'hospital_name': h_name,
161 'travel_time_min': 15.0,
162 'geometry': isochrone_poly
163 })
164 else:
165 # If no nodes reachable, we still need a row
166 isochrone_data.append({
167 'hospital_name': h_name,
168 'travel_time_min': 0.0,
169 'geometry': Point(hosp.geometry.x, hosp.geometry.y).buffer(0)
170 })
171
172 isochrones_gdf = gpd.GeoDataFrame(isochrone_data, crs=CRS)
173
174 # 5. Save to GPKG
175 print("Saving to GPKG...")
176 # Prepare incidents layer
177 incidents_output = incidents_gdf[['incident_id', 'geometry']]
178
179 with gpd.GeoDataFrame(columns=incidents_output.columns, crs=CRS) as empty_gdf:
180 # We can't use 'with' like that for writing. Use the engine.
181 pass
182
183 # Write layers
184 # Layer 1: incidents
185 incidents_output.to_file('emergency_routing.gpkg', layer='incidents', driver='GPKG')
186 # Layer 2: closest_hospital
187 closest_hospital_gdf.to_file('emergency_routing.gpkg', layer='closest_hospital', driver='GPKG')
188 # Layer 3: distance_matrix
189 distance_matrix_gdf.to_file('emergency_routing.gpkg', layer='distance_matrix', driver='GPKG')
190 # Layer 4: isochrones_15min
191 isochrones_gdf.to_file('emergency_routing.gpkg', layer='isochrones_15min', driver='GPKG')
192
193 print("Done!")
194
195 if __name__ == "__main__":
196 solve()
0:52
://thought thought
0:52
Edit
solve.py
old_string not found in solve.py
1:25
Write
solve.py
content · 190 lines · py
import osmnx as ox
import networkx as nx
import geopandas as gpd
import pandas as pd
from shapely.geometry import Point, LineString, MultiLineString, Polygon, MultiPolygon
import numpy as np
# Constants
BBOX = (48.83, 2.30, 48.88, 2.38) # min_lat, min_lon, max_lat, max_lon
CRS = "EPSG:2154" # RGF93 / Lambert-93 (Official metric CRS for France/Paris)
SPEED_DEFAULT = 30 / 3.6 # 30 km/h to m/s
def solve():
# 1. Load Incidents
incidents_df = pd.read_csv('incidents.csv')
# Create GeoDataFrame
incidents_gdf = gpd.GeoDataFrame(
incidents_df,
geometry=gpd.points_from_xy(incidents_df.longitude, incidents_df.latitude),
crs="EPSG:4326"
).to_crs(CRS)
# 2. Fetch OSM Data (Network and Hospitals)
print("Fetching OSM network and hospitals...")
# Use osmnx to get the graph
# bbox: north, south, east, west
north, south, east, west = BBOX[2], BBOX[0], BBOX[3], BBOX[1]
G = ox.graph_from_bbox(north, south, east, west, network_type='drive')
# Apply speeds
G = ox.add_edge_speeds(G, fallback=30)
G = ox.add_edge_travel_times(G)
# Fetch Hospitals using Overpass (via osmnx or direct query if needed)
# osmnx doesn't have a direct 'hospitals' function but we can use query
tags = {"amenity": "hospital"}
hospitals_gdf = ox.features_from_bbox(north, south, east, west, tags=tags)
# Clean hospitals: ensure they are points/polygons and have names
# Some might be Multipolygons or Points
hospitals_gdf = hospitals_gdf.copy()
if 'name' not in hospitals_gdf.columns:
hospitals_gdf['name'] = 'Unknown Hospital'
hospitals_gdf['hospital_name'] = hospitals_gdf['name'].fillna('Unknown Hospital')
# We want the centroid for routing if they are polygons
hospitals_gdf['geometry'] = hospitals_gdf.geometry.centroid
hospitals_gdf = hospitals_gdf.to_crs(CRS)
# 3. Routing Analysis
print("Calculating routes...")
# Pre-calculate nearest nodes for incidents and hospitals
inc_nodes = ox.distance.nearest_nodes(G, incidents_gdf.geometry.x, incidents_gdf.geometry.y)
hosp_nodes = ox.distance.nearest_nodes(G, hospitals_gdf.geometry.x, hospitals_gdf.geometry.y)
closest_routes_data = []
distance_matrix_data = []
for idx, inc in incidents_gdf.iterrows():
inc_id = inc['incident_id']
u = inc_nodes[idx]
# Calculate distances to all hospitals for this incident
hosp_distances = []
for h_idx, hosp in hospitals_gdf.iterrows():
v = hosp_nodes[h_idx]
try:
# weight='length' to get network distance
d = nx.shortest_path_length(G, u, v, weight='length')
hosp_distances.append({
'hospital_name': hosp['hospital_name'],
'network_distance_m': d,
'hosp_idx': h_idx,
'v_node': v
})
except nx.NetworkXNoPath:
continue
# Sort by distance
hosp_distances.sort(key=lambda x: x['network_distance_m'])
# Closest hospital for 'closest_hospital' layer
if hosp_distances:
closest = hosp_distances[0]
# Get path for the line geometry
path = nx.shortest_path(G, u, closest['v_node'], weight='length')
# Convert path to LineString
nodes_coords = []
for node in path:
# Get lat/lon from graph node to ensure accurate geometry
nodes_coords.append((G.nodes[node]['x'], G.nodes[node]['y']))
route_line = LineString(nodes_coords)
closest_routes_data.append({
'incident_id': inc_id,
'hospital_name': closest['hospital_name'],
'network_distance_m': float(closest['network_distance_m']),
'geometry': route_line
})
# Distance matrix (top 3)
for rank, h_info in enumerate(hosp_distances[:3], 1):
distance_matrix_data.append({
'incident_id': inc_id,
'hospital_name': h_info['hospital_name'],
'rank': rank,
'network_distance_m': float(h_info['network_distance_m'])
})
else:
# Fallback if no hospital reachable
distance_matrix_data.append({
'incident_id': inc_id,
'hospital_name': 'None',
'rank': 1,
'network_distance_m': np.nan
})
# Convert to GeoDataFrames
closest_hospital_gdf = gpd.GeoDataFrame(closest_routes_data, crs=CRS)
distance_matrix_gdf = gpd.GeoDataFrame(distance_matrix_data, crs=CRS)
# distance_matrix_gdf.geometry = None # explicitly empty
# 4. Isochrones (15 min)
print("Calculating isochrones...")
isochrone_data = []
# 15 minutes = 900 seconds
time_limit = 15 * 60
for idx, hosp in hospitals_gdf.iterrows():
h_name = hosp['hospital_name']
v = hosp_nodes[idx]
# Subgraph of reachable nodes within time limit
# Use Dijkstra to find all nodes within time_limit
subgraph_nodes = nx.single_source_dijkstra_path_length(G, v, cutoff=time_limit, weight='travel_time')
if subgraph_nodes:
# Get geometries of these nodes
node_points = [Point(G.nodes[n]['x'], G.nodes[n]['y']) for n in subgraph_nodes.keys()]
# Create a concave hull or buffer?
# Usually, isochrones are polygons around the nodes.
# A simple approach is to buffer the nodes and take the union,
# or better, use the edges.
# Let's use the nodes and create a buffer/concave hull.
# For speed and robustness, we'll buffer the points and dissolve.
points_gdf = gpd.GeoDataFrame(geometry=node_points, crs=CRS)
# Buffer points by a small amount to ensure connectivity in union
try:
isochrone_poly = points_gdf.buffer(20).union_all()
except AttributeError:
isochrone_poly = points_gdf.buffer(20).unary_union
isochrone_data.append({
'hospital_name': h_name,
'travel_time_min': 15.0,
'geometry': isochrone_poly
})
else:
# If no nodes reachable, we still need a row
isochrone_data.append({
'hospital_name': h_name,
'travel_time_min': 0.0,
'geometry': Point(hosp.geometry.x, hosp.geometry.y).buffer(0)
})
isochrones_gdf = gpd.GeoDataFrame(isochrone_data, crs=CRS)
# 5. Save to GPKG
print("Saving to GPKG...")
# Prepare incidents layer
incidents_output = incidents_gdf[['incident_id', 'geometry']]
# Write layers
# Layer 1: incidents
incidents_output.to_file('emergency_routing.gpkg', layer='incidents', driver='GPKG')
# Layer 2: closest_hospital
closest_hospital_gdf.to_file('emergency_routing.gpkg', layer='closest_hospital', driver='GPKG')
# $\_fix_
# Layer 3: distance_matrix
distance_matrix_gdf.to_file('emergency_routing.gpkg', layer='distance_matrix', driver='GPKG')
# Layer 4: isochrones_15min
isochrones_gdf.to_file('emergency_routing.gpkg', layer='isochrones_15min', driver='GPKG')
print("Done!")
if __name__ == "__main__":
solve()
Wrote solve.py (7557 chars).
1:27
Bash
python solve.py
Fetching OSM network and hospitals...
Traceback (most recent call last):
File "/work/solve.py", line 189, in <module>
solve()
~~~~~^^
File "/work/solve.py", line 28, 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:36
Edit
solve.py
old_string not found in solve.py
1:37
Read
solve.py
1 import osmnx as ox
2 import networkx as nx
3 import geopandas as gpd
4 import pandas as pd
5 from shapely.geometry import Point, LineString, MultiLineString, Polygon, MultiPolygon
6 import numpy as np
7
8 # Constants
9 BBOX = (48.83, 2.30, 48.88, 2.38) # min_lat, min_lon, max_lat, max_lon
10 CRS = "EPSG:2154" # RGF93 / Lambert-93 (Official metric CRS for France/Paris)
11 SPEED_DEFAULT = 30 / 3.6 # 30 km/h to m/s
12
13 def solve():
14 # 1. Load Incidents
15 incidents_df = pd.read_csv('incidents.csv')
16 # Create GeoDataFrame
17 incidents_gdf = gpd.GeoDataFrame(
18 incidents_df,
19 geometry=gpd.points_from_xy(incidents_df.longitude, incidents_df.latitude),
20 crs="EPSG:4326"
21 ).to_crs(CRS)
22
23 # 2. Fetch OSM Data (Network and Hospitals)
24 print("Fetching OSM network and hospitals...")
25 # Use osmnx to get the graph
26 # bbox: north, south, east, west
27 north, south, east, west = BBOX[2], BBOX[0], BBOX[3], BBOX[1]
28 G = ox.graph_from_bbox(north, south, east, west, network_type='drive')
29
30 # Apply speeds
31 G = ox.add_edge_speeds(G, fallback=30)
32 G = ox.add_edge_travel_times(G)
33
34 # Fetch Hospitals using Overpass (via osmnx or direct query if needed)
35 # osmnx doesn't have a direct 'hospitals' function but we can use query
36 tags = {"amenity": "hospital"}
37 hospitals_gdf = ox.features_from_bbox(north, south, east, west, tags=tags)
38
39 # Clean hospitals: ensure they are points/polygons and have names
40 # Some might be Multipolygons or Points
41 hospitals_gdf = hospitals_gdf.copy()
42 if 'name' not in hospitals_gdf.columns:
43 hospitals_gdf['name'] = 'Unknown Hospital'
44 hospitals_gdf['hospital_name'] = hospitals_gdf['name'].fillna('Unknown Hospital')
45
46 # We want the centroid for routing if they are polygons
47 hospitals_gdf['geometry'] = hospitals_gdf.geometry.centroid
48 hospitals_gdf = hospitals_gdf.to_crs(CRS)
49
50 # 3. Routing Analysis
51 print("Calculating routes...")
52
53 # Pre-calculate nearest nodes for incidents and hospitals
54 inc_nodes = ox.distance.nearest_nodes(G, incidents_gdf.geometry.x, incidents_gdf.geometry.y)
55 hosp_nodes = ox.distance.nearest_nodes(G, hospitals_gdf.geometry.x, hospitals_gdf.geometry.y)
56
57 closest_routes_data = []
58 distance_matrix_data = []
59
60 for idx, inc in incidents_gdf.iterrows():
61 inc_id = inc['incident_id']
62 u = inc_nodes[idx]
63
64 # Calculate distances to all hospitals for this incident
65 hosp_distances = []
66 for h_idx, hosp in hospitals_gdf.iterrows():
67 v = hosp_nodes[h_idx]
68 try:
69 # weight='length' to get network distance
70 d = nx.shortest_path_length(G, u, v, weight='length')
71 hosp_distances.append({
72 'hospital_name': hosp['hospital_name'],
73 'network_distance_m': d,
74 'hosp_idx': h_idx,
75 'v_node': v
76 })
77 except nx.NetworkXNoPath:
78 continue
79
80 # Sort by distance
81 hosp_distances.sort(key=lambda x: x['network_distance_m'])
82
83 # Closest hospital for 'closest_hospital' layer
84 if hosp_distances:
85 closest = hosp_distances[0]
86 # Get path for the line geometry
87 path = nx.shortest_path(G, u, closest['v_node'], weight='length')
88 # Convert path to LineString
89 nodes_coords = []
90 for node in path:
91 # Get lat/lon from graph node to ensure accurate geometry
92 nodes_coords.append((G.nodes[node]['x'], G.nodes[node]['y']))
93
94 route_line = LineString(nodes_coords)
95
96 closest_routes_data.append({
97 'incident_id': inc_id,
98 'hospital_name': closest['hospital_name'],
99 'network_distance_m': float(closest['network_distance_m']),
100 'geometry': route_line
101 })
102
103 # Distance matrix (top 3)
104 for rank, h_info in enumerate(hosp_distances[:3], 1):
105 distance_matrix_data.append({
106 'incident_id': inc_id,
107 'hospital_name': h_info['hospital_name'],
108 'rank': rank,
109 'network_distance_m': float(h_info['network_distance_m'])
110 })
111 else:
112 # Fallback if no hospital reachable
113 distance_matrix_data.append({
114 'incident_id': inc_id,
115 'hospital_name': 'None',
116 'rank': 1,
117 'network_distance_m': np.nan
118 })
119
120 # Convert to GeoDataFrames
121 closest_hospital_gdf = gpd.GeoDataFrame(closest_routes_data, crs=CRS)
122 distance_matrix_gdf = gpd.GeoDataFrame(distance_matrix_data, crs=CRS)
123 # distance_matrix_gdf.geometry = None # explicitly empty
124
125 # 4. Isochrones (15 min)
126 print("Calculating isochrones...")
127 isochrone_data = []
128 # 15 minutes = 900 seconds
129 time_limit = 15 * 60
130
131 for idx, hosp in hospitals_gdf.iterrows():
132 h_name = hosp['hospital_name']
133 v = hosp_nodes[idx]
134
135 # Subgraph of reachable nodes within time limit
136 # Use Dijkstra to find all nodes within time_limit
137 subgraph_nodes = nx.single_source_dijkstra_path_length(G, v, cutoff=time_limit, weight='travel_time')
138
139 if subgraph_nodes:
140 # Get geometries of these nodes
141 node_points = [Point(G.nodes[n]['x'], G.nodes[n]['y']) for n in subgraph_nodes.keys()]
142 # Create a concave hull or buffer?
143 # Usually, isochrones are polygons around the nodes.
144 # A simple approach is to buffer the nodes and take the union,
145 # or better, use the edges.
146 # Let's use the nodes and create a buffer/concave hull.
147 # For speed and robustness, we'll buffer the points and dissolve.
148 points_gdf = gpd.GeoDataFrame(geometry=node_points, crs=CRS)
149 # Buffer points by a small amount to ensure connectivity in union
150 try:
151 isochrone_poly = points_gdf.buffer(20).union_all()
152 except AttributeError:
153 isochrone_poly = points_gdf.buffer(20).unary_union
154
155 isochrone_data.append({
156 'hospital_name': h_name,
157 'travel_time_min': 15.0,
158 'geometry': isochrone_poly
159 })
160 else:
161 # If no nodes reachable, we still need a row
162 isochrone_data.append({
163 'hospital_name': h_name,
164 'travel_time_min': 0.0,
165 'geometry': Point(hosp.geometry.x, hosp.geometry.y).buffer(0)
166 })
167
168 isochrones_gdf = gpd.GeoDataFrame(isochrone_data, crs=CRS)
169
170 # 5. Save to GPKG
171 print("Saving to GPKG...")
172 # Prepare incidents layer
173 incidents_output = incidents_gdf[['incident_id', 'geometry']]
174
175 # Write layers
176 # Layer 1: incidents
177 incidents_output.to_file('emergency_routing.gpkg', layer='incidents', driver='GPKG')
178 # Layer 2: closest_hospital
179 closest_hospital_gdf.to_file('emergency_routing.gpkg', layer='closest_hospital', driver='GPKG')
180 # $\_fix_
181 # Layer 3: distance_matrix
182 distance_matrix_gdf.to_file('emergency_routing.gpkg', layer='distance_matrix', driver='GPKG')
183 # Layer 4: isochrones_15min
184 isochrones_gdf.to_file('emergency_routing.gpkg', layer='isochrones_15min', driver='GPKG')
185
186 print("Done!")
187
188 if __name__ == "__main__":
189 solve()
1:41
Edit
solve.py
old_string not found in solve.py
2:13
Write
solve.py
content · 189 lines · py
import osmnx as ox
import networkx as nx
import geopandas as gpd
import pandas as pd
from shapely.geometry import Point, LineString, MultiLineString, Polygon, MultiPolygon
import numpy as np
# Constants
BBOX = (48.83, 2.30, 48.88, 2.38) # min_lat, min_lon, max_lat, max_lon
CRS = "EPSG:2154" # RGF93 / Lambert-93 (Official metric CRS for France/Paris)
SPEED_DEFAULT = 30 / 3.6 # 30 km/h to m/s
def solve():
# 1. Load Incidents
incidents_df = pd.read_csv('incidents.csv')
# Create GeoDataFrame
incidents_gdf = gpd.GeoDataFrame(
incidents_df,
geometry=gpd.points_from_xy(incidents_df.longitude, incidents_df.latitude),
crs="EPSG:4326"
).to_crs(CRS)
# 2. Fetch OSM Data (Network and Hospitals)
print("Fetching OSM network and hospitals...")
# Use osmnx to get the graph
# bbox: north, south, east, west
north, south, east, west = BBOX[2], BBOX[0], BBOX[3], BBOX[1]
G = ox.graph_from_bbox(north=north, south=south, east=east, west=west, network_type='drive')
# Apply speeds
G = ox.add_edge_speeds(G, fallback=30)
G = ox.add_edge_travel_times(G)
# Fetch Hospitals using Overpass (via osmnx or direct query if needed)
# osmnx doesn't have a direct 'hospitals' function but we can use query
tags = {"amenity": "hospital"}
hospitals_gdf = ox.features_from_bbox(north=north, south=south, east=east, west=west, tags=tags)
# Clean hospitals: ensure they are points/polygons and have names
# Some might be Multipolygons or Points
hospitals_gdf = hospitals_gdf.copy()
if 'name' not in hospitals_gdf.columns:
hospitals_gdf['name'] = 'Unknown Hospital'
hospitals_gdf['hospital_name'] = hospitals_gdf['name'].fillna('Unknown Hospital')
# We want the centroid for routing if they are polygons
hospitals_gdf['geometry'] = hospitals_gdf.geometry.centroid
hospitals_gdf = hospitals_gdf.to_crs(CRS)
# 3. Routing Analysis
print("Calculating routes...")
# Pre-calculate nearest nodes for incidents and hospitals
inc_nodes = ox.distance.nearest_nodes(G, incidents_gdf.geometry.x, incidents_gdf.geometry.y)
hosp_nodes = ox.distance.nearest_nodes(G, hospitals_gdf.geometry.x, hospitals_gdf.geometry.y)
closest_routes_data = []
distance_matrix_data = []
for idx, inc in incidents_gdf.iterrows():
inc_id = inc['incident_id']
u = inc_nodes[idx]
# Calculate distances to all hospitals for this incident
hosp_distances = []
for h_idx, hosp in hospitals_gdf.iterrows():
v = hosp_nodes[h_idx]
try:
# weight='length' to get network distance
d = nx.shortest_path_length(G, u, v, weight='length')
hosp_distances.append({
'hospital_name': hosp['hospital_name'],
'network_distance_m': d,
'hosp_idx': h_idx,
'v_node': v
})
except nx.NetworkXNoPath:
continue
# Sort by distance
hosp_distances.sort(key=lambda x: x['network_distance_m'])
# Closest hospital for 'closest_hospital' layer
if hosp_distances:
closest = hosp_distances[0]
# Get path for the line geometry
path = nx.shortest_path(G, u, closest['v_node'], weight='length')
# Convert path to LineString
nodes_coords = []
for node in path:
# Get lat/lon from graph node to ensure accurate geometry
nodes_coords.append((G.nodes[node]['x'], G.nodes[node]['y']))
route_line = LineString(nodes_coords)
closest_routes_data.append({
'incident_id': inc_id,
'hospital_name': closest['hospital_name'],
'network_distance_m': float(closest['network_distance_m']),
'geometry': route_line
})
# Distance matrix (top 3)
for rank, h_info in enumerate(hosp_distances[:3], 1):
distance_matrix_data.append({
'incident_id': inc_id,
'hospital_name': h_info['hospital_name'],
'rank': rank,
'network_distance_m': float(h_info['network_distance_m'])
})
else:
# Fallback if no hospital reachable
distance_matrix_data.append({
'incident_id': inc_id,
'hospital_name': 'None',
'rank': 1,
'network_distance_m': np.nan
})
# Convert to GeoDataFrames
closest_hospital_gdf = gpd.GeoDataFrame(closest_routes_data, crs=CRS)
distance_matrix_gdf = gpd.GeoDataFrame(distance_matrix_data, crs=CRS)
# distance_matrix_gdf.geometry = None # explicitly empty
# 4. Isochrones (15 min)
print("Calculating isochrones...")
isochrone_data = []
# 15 minutes = 900 seconds
time_limit = 15 * 60
for idx, hosp in hospitals_gdf.iterrows():
h_name = hosp['hospital_name']
v = hosp_nodes[idx]
# Subgraph of reachable nodes within time limit
# Use Dijkstra to find all nodes within time_limit
subgraph_nodes = nx.single_source_dijkstra_path_length(G, v, cutoff=time_limit, weight='travel_time')
if subgraph_nodes:
# Get geometries of these nodes
node_points = [Point(G.nodes[n]['x'], G.nodes[n]['y']) for n in subgraph_nodes.keys()]
# Create a concave hull or buffer?
# Usually, isochrones are polygons around the nodes.
# A simple approach is to buffer the nodes and take the union,
# or better, use the edges.
# Let's use the nodes and create a buffer/concave hull.
# For speed and robustness, we'll buffer the points and dissolve.
points_gdf = gpd.GeoDataFrame(geometry=node_points, crs=CRS)
# Buffer points by a small amount to ensure connectivity in union
try:
isochrone_poly = points_gdf.buffer(20).union_all()
except AttributeError:
isochrone_poly = points_gdf.buffer(20).unary_union
isochrone_data.append({
'hospital_name': h_name,
'travel_time_min': 15.0,
'geometry': isochrone_poly
})
else:
# If no nodes reachable, we still need a row
isochrone_data.append({
'hospital_name': h_name,
'travel_time_min': 0.0,
'geometry': Point(hosp.geometry.x, hosp.geometry.y).buffer(0)
})
isochrones_gdf = gpd.GeoDataFrame(isochrone_data, crs=CRS)
# 5. Save to GPKG
print("Saving to GPKG...")
# Prepare incidents layer
incidents_output = incidents_gdf[['incident_id', 'geometry']]
# Write layers
# Layer 1: incidents
incidents_output.to_file('emergency_routing.gpkg', layer='incidents', driver='GPKG')
# Layer 2: closest_hospital
closest_hospital_gdf.to_file('emergency_routing.gpkg', layer='closest_hospital', driver='GPKG')
# Layer 3: distance_matrix
distance_matrix_gdf.to_file('emergency_routing.gpkg', layer='distance_matrix', driver='GPKG')
# Layer 4: isochrones_15min
isochrones_gdf.to_file('emergency_routing.gpkg', layer='isochrones_15min', driver='GPKG')
print("Done!")
if __name__ == "__main__":
solve()
Wrote solve.py (7587 chars).
2:15
Bash
python solve.py
Fetching OSM network and hospitals...
Traceback (most recent call last):
File "/work/solve.py", line 188, in <module>
solve()
~~~~~^^
File "/work/solve.py", line 28, 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'
2:20
Bash
python -c "import osmnx; print(osmnx.__version__)"
2.1.0
2:24
Bash
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.
2:34
Edit
solve.py
Edited solve.py.
2:36
Bash
python solve.py
<exit_code>-1</exit_code> <error>Command timed out after 120000ms</error>
4:41
(no output)
Bash
# No way to increase timeout in Bash tool directly,
# No way to increase timeout in Bash tool directly, # but I can try to run it and see if it's just slow. # I will try to optimize the code if I can, but first let's see if it works with more time. # Since I can't, I'll try to make the code more efficient. # One optimization: for isochrones, we don't need to buffer all points if they are very close. # Also, checking if the hospital is already in the graph nodes. # The main issue might be the Overpass query or the Dijkstra. # Let's try to reduce the bbox slightly or check if it's a large area. # The area is about 0.05 deg, which is ~5km. That should be fast. # Maybe the hospital query is slow.
5:02
Edit
solve.py
old_string not found in solve.py
5:04
Read
solve.py
1 import osmnx as ox
2 import networkx as nx
3 import geopandas as gpd
4 import pandas as pd
5 from shapely.geometry import Point, LineString, MultiLineString, Polygon, MultiPolygon
6 import numpy as np
7
8 # Constants
9 BBOX = (48.83, 2.30, 48.88, 2.38) # min_lat, min_lon, max_lat, max_lon
10 CRS = "EPSG:2154" # RGF93 / Lambert-93 (Official metric CRS for France/Paris)
11 SPEED_DEFAULT = 30 / 3.6 # 30 km/h to m/s
12
13 def solve():
14 # 1. Load Incidents
15 incidents_df = pd.read_csv('incidents.csv')
16 # Create GeoDataFrame
17 incidents_gdf = gpd.GeoDataFrame(
18 incidents_df,
19 geometry=gpd.points_from_xy(incidents_df.longitude, incidents_df.latitude),
20 crs="EPSG:4326"
21 ).to_crs(CRS)
22
23 # 2. Fetch OSM Data (Network and Hospitals)
24 print("Fetching OSM network and hospitals...")
25 # Use osmnx to get the graph
26 # bbox: (left, bottom, right, top) -> (west, south, east, north)
27 west, south, east, north = BBOX[1], BBOX[0], BBOX[3], BBOX[2]
28 G = ox.graph_from_bbox(bbox=(west, south, east, north), network_type='drive')
29
30 # Apply speeds
31 G = ox.add_edge_speeds(G, fallback=30)
32 G = ox.add_edge_travel_times(G)
33
34 # Fetch Hospitals using Overpass (via osmnx or direct query if needed)
35 # osmnx doesn't have a direct 'hospitals' function but we can use query
36 tags = {"amenity": "hospital"}
37 hospitals_gdf = ox.features_from_bbox(bbox=(west, south, east, north), tags=tags)
38
39 # Clean hospitals: ensure they are points/polygons and have names
40 # Some might be Multipolygons or Points
41 hospitals_gdf = hospitals_gdf.copy()
42 if 'name' not in hospitals_gdf.columns:
43 hospitals_gdf['name'] = 'Unknown Hospital'
44 hospitals_gdf['hospital_name'] = hospitals_gdf['name'].fillna('Unknown Hospital')
45
46 # We want the centroid for routing if they are polygons
47 hospitals_gdf['geometry'] = hospitals_gdf.geometry.centroid
48 hospitals_gdf = hospitals_gdf.to_crs(CRS)
49
50 # 3. Routing Analysis
51 print("Calculating routes...")
52
53 # Pre-calculate nearest nodes for incidents and hospitals
54 inc_nodes = ox.distance.nearest_nodes(G, incidents_gdf.geometry.x, incidents_gdf.geometry.y)
55 hosp_nodes = ox.distance.nearest_nodes(G, hospitals_gdf.geometry.x, hospitals_gdf.geometry.y)
56
57 closest_routes_data = []
58 distance_matrix_data = []
59
60 for idx, inc in incidents_gdf.iterrows():
61 inc_id = inc['incident_id']
62 u = inc_nodes[idx]
63
64 # Calculate distances to all hospitals for this incident
65 hosp_distances = []
66 for h_idx, hosp in hospitals_gdf.iterrows():
67 v = hosp_nodes[h_idx]
68 try:
69 # weight='length' to get network distance
70 d = nx.shortest_path_length(G, u, v, weight='length')
71 hosp_distances.append({
72 'hospital_name': hosp['hospital_name'],
73 'network_distance_m': d,
74 'hosp_idx': h_idx,
75 'v_node': v
76 })
77 except nx.NetworkXNoPath:
78 continue
79
80 # Sort by distance
81 hosp_distances.sort(key=lambda x: x['network_distance_m'])
82
83 # Closest hospital for 'closest_hospital' layer
84 if hosp_distances:
85 closest = hosp_distances[0]
86 # Get path for the line geometry
87 path = nx.shortest_path(G, u, closest['v_node'], weight='length')
88 # Convert path to LineString
89 nodes_coords = []
90 for node in path:
91 # Get lat/lon from graph node to ensure accurate geometry
92 nodes_coords.append((G.nodes[node]['x'], G.nodes[node]['y']))
93
94 route_line = LineString(nodes_coords)
95
96 closest_routes_data.append({
97 'incident_id': inc_id,
98 'hospital_name': closest['hospital_name'],
99 'network_distance_m': float(closest['network_distance_m']),
100 'geometry': route_line
101 })
102
103 # Distance matrix (top 3)
104 for rank, h_info in enumerate(hosp_distances[:3], 1):
105 distance_matrix_data.append({
106 'incident_id': inc_id,
107 'hospital_name': h_info['hospital_name'],
108 'rank': rank,
109 'network_distance_m': float(h_info['network_distance_m'])
110 })
111 else:
112 # Fallback if no hospital reachable
113 distance_matrix_data.append({
114 'incident_id': inc_id,
115 'hospital_name': 'None',
116 'rank': 1,
117 'network_distance_m': np.nan
118 })
119
120 # Convert to GeoDataFrames
121 closest_hospital_gdf = gpd.GeoDataFrame(closest_routes_data, crs=CRS)
122 distance_matrix_gdf = gpd.GeoDataFrame(distance_matrix_data, crs=CRS)
123 # distance_matrix_gdf.geometry = None # explicitly empty
124
125 # 4. Isochrones (15 min)
126 print("Calculating isochrones...")
127 isochrone_data = []
128 # 15 minutes = 900 seconds
129 time_limit = 15 * 60
130
131 for idx, hosp in hospitals_gdf.iterrows():
132 h_name = hosp['hospital_name']
133 v = hosp_nodes[idx]
134
135 # Subgraph of reachable nodes within time limit
136 # Use Dijkstra to find all nodes within time_limit
137 subgraph_nodes = nx.single_source_dijkstra_path_length(G, v, cutoff=time_limit, weight='travel_time')
138
139 if subgraph_nodes:
140 # Get geometries of these nodes
141 node_points = [Point(G.nodes[n]['x'], G.nodes[n]['y']) for n in subgraph_nodes.keys()]
142 # Create a concave hull or buffer?
143 # Usually, isochrones are polygons around the nodes.
144 # A simple approach is to buffer the nodes and take the union,
145 # or better, use the edges.
146 # Let's use the nodes and create a buffer/concave hull.
147 # For speed and robustness, we'll buffer the points and dissolve.
148 points_gdf = gpd.GeoDataFrame(geometry=node_points, crs=CRS)
149 # Buffer points by a small amount to ensure connectivity in union
150 try:
151 isochrone_poly = points_gdf.buffer(20).union_all()
152 except AttributeError:
153 isochrone_poly = points_gdf.buffer(20).unary_union
154
155 isochrone_data.append({
156 'hospital_name': h_name,
157 'travel_time_min': 15.0,
158 'geometry': isochrone_poly
159 })
160 else:
161 # If no nodes reachable, we still need a row
162 isochrone_data.append({
163 'hospital_name': h_name,
164 'travel_time_min': 0.0,
165 'geometry': Point(hosp.geometry.x, hosp.geometry.y).buffer(0)
166 })
167
168 isochrones_gdf = gpd.GeoDataFrame(isochrone_data, crs=CRS)
169
170 # 5. Save to GPKG
171 print("Saving to GPKG...")
172 # Prepare incidents layer
173 incidents_output = incidents_gdf[['incident_id', 'geometry']]
174
175 # Write layers
176 # Layer 1: incidents
177 incidents_output.to_file('emergency_routing.gpkg', layer='incidents', driver='GPKG')
178 # Layer 2: closest_hospital
179 closest_hospital_gdf.to_file('emergency_routing.gpkg', layer='closest_hospital', driver='GPKG')
180 # Layer 3: distance_matrix
181 distance_matrix_gdf.to_file('emergency_routing.gpkg', layer='distance_matrix', driver='GPKG')
182 # Layer 4: isochrones_15min
183 isochrones_gdf.to_file('emergency_routing.gpkg', layer='isochrones_15min', driver='GPKG')
184
185 print("Done!")
186
187 if __name__ == "__main__":
188 solve()
5:36
Write
solve.py
content · 175 lines · py
import osmnx as ox
import networkx as nx
import geopandas as gpd
import pandas as pd
from shapely.geometry import Point, LineString, MultiLineString, Polygon, MultiPolygon
import numpy as np
# Constants
BBOX = (48.83, 2.30, 48.88, 2.38) # min_lat, min_lon, max_lat, max_lon
CRS = "EPSG:2154" # RGF93 / Lambert-93 (Official metric CRS for France/Paris)
SPEED_DEFAULT = 30 / 3.6 # 30 km/h to m/s
def solve():
# 1. Load Incidents
incidents_df = pd.read_csv('incidents.csv')
# Create GeoDataFrame
incidents_gdf = gpd.GeoDataFrame(
incidents_df,
geometry=gpd.points_from_xy(incidents_df.longitude, incidents_df.latitude),
crs="EPSG:4326"
).to_crs(CRS)
# 2. Fetch OSM Data (Network and Hospitals)
print("Fetching OSM network and hospitals...")
# Use osmnx to get the graph
# bbox: (left, bottom, right, top) -> (west, south, east, north)
west, south, east, north = BBOX[1], BBOX[0], BBOX[3], BBOX[2]
G = ox.graph_from_bbox(bbox=(west, south, east, north), network_type='drive')
# Apply speeds
G = ox.add_edge_speeds(G, fallback=30)
G = ox.add_edge_travel_times(G)
# Fetch Hospitals using Overpass (via osmnx or direct query if needed)
tags = {"amenity": "hospital"}
hospitals_gdf = ox.features_from_bbox(bbox=(west, south, east, north), tags=tags)
# Clean hospitals: ensure they are points/polygons and have names
hospitals_gdf = hospitals_gdf.copy()
if 'name' not in hospitals_gdf.columns:
hospitals_gdf['name'] = 'Unknown Hospital'
hospitals_gdf['hospital_name'] = hospitals_gdf['name'].fillna('Unknown Hospital')
# We want the centroid for routing if they are polygons
hospitals_gdf['geometry'] = hospitals_gdf.geometry.centroid
hospitals_gdf = hospitals_gdf.to_crs(CRS)
# 3. Routing Analysis
print("Calculating routes...")
# Pre-calculate nearest nodes for incidents and hospitals
inc_nodes = ox.distance.nearest_nodes(G, incidents_gdf.geometry.x, incidents_gdf.geometry.y)
hosp_nodes = ox.distance.nearest_nodes(G, hospitals_gdf.geometry.x, hospitals_gdf.geometry.y)
closest_routes_data = []
distance_matrix_data = []
for idx, inc in incidents_gdf.iterrows():
inc_id = inc['incident_id']
u = inc_nodes[idx]
# Calculate distances to all hospitals for this incident
hosp_distances = []
for h_idx, hosp in hospitals_gdf.iterrows():
v = hosp_nodes[h_idx]
try:
# weight='length' to get network distance
d = nx.shortest_path_length(G, u, v, weight='length')
hosp_distances.append({
'hospital_name': hosp['hospital_name'],
'network_distance_m': d,
'hosp_idx': h_idx,
'v_node': v
})
except nx.NetworkXNoPath:
continue
# Sort by distance
hosp_distances.sort(key=lambda x: x['network_distance_m'])
# Closest hospital for 'closest_hospital' layer
if hosp_distances:
closest = hosp_distances[0]
# Get path for the line geometry
path = nx.shortest_path(G, u, closest['v_node'], weight='length')
# Convert path to LineString
nodes_coords = []
for node in path:
nodes_coords.append((G.nodes[node]['x'], G.nodes[node]['y']))
route_line = LineString(nodes_coords)
closest_routes_data.append({
'incident_id': inc_id,
'hospital_name': closest['hospital_name'],
'network_distance_m': float(closest['network_distance_m']),
'geometry': route_line
})
# Distance matrix (top 3)
for rank, h_info in enumerate(hosp_distances[:3], 1):
distance_matrix_data.append({
'incident_id': inc_id,
'hospital_name': h_info['hospital_name'],
'rank': rank,
'network_distance_m': float(h_info['network_distance_m'])
})
else:
# Fallback if no hospital reachable
distance_matrix_data.append({
'incident_id': inc_id,
'hospital_name': 'None',
'rank': 1,
'network_distance_m': np.nan
})
# Convert to GeoDataFrames
closest_hospital_gdf = gpd.GeoDataFrame(closest_routes_data, crs=CRS)
distance_matrix_gdf = gpd.GeoDataFrame(distance_matrix_data, crs=CRS)
# 4. Isochrones (15 min)
print("Calculating isochrones...")
isochrone_data = []
# 15 minutes = 900 seconds
time_limit = 15 * 60
for idx, hosp in hospitals_gdf.iterrows():
h_name = hosp['hospital_name']
v = hosp_nodes[idx]
# Subgraph of reachable nodes within time limit
try:
subgraph_nodes = nx.single_source_dijkstra_path_length(G, v, cutoff=time_limit, weight='travel_time')
except (nx.NodeNotFound, nx.NetworkXError):
subgraph_nodes = {}
if subgraph_nodes:
# Get geometries of these nodes
node_points = [Point(G.nodes[n]['x'], G.nodes[n]['y']) for n in subgraph_nodes.keys()]
points_gdf = gpd.GeoDataFrame(geometry=node_points, crs=CRS)
try:
isochrone_poly = points_gdf.buffer(25).union_all()
except AttributeError:
isochrone_poly = points_gdf.buffer(25).unary_union
isochrone_data.append({
'hospital_name': h_name,
'travel_time_min': 15.0,
'geometry': isochrone_poly
})
else:
isochrone_data.append({
'hospital_name': h_name,
'travel_time_min': 0.0,
'geometry': Point(hosp.geometry.x, hosp.geometry.y).buffer(10)
})
isochrones_gdf = gpd.GeoDataFrame(isochrone_data, crs=CRS)
# 5. Save to GPKG
print("Saving to GPKG...")
# Prepare incidents layer
incidents_output = incidents_gdf[['incident_id', 'geometry']]
# Write layers
incidents_output.to_file('emergency_routing.gpkg', layer='incidents', driver='GPKG')
closest_hospital_gdf.to_file('emergency_routing.gpkg', layer='closest_hospital', driver='GPKG')
distance_matrix_gdf.to_file('emergency_routing.gpkg', layer='distance_matrix', driver='GPKG')
isochrones_gdf.to_file('emergency_routing.gpkg', layer='isochrones_15min', driver='GPKG')
print("Done!")
if __name__ == "__main__":
solve()
Wrote solve.py (6747 chars).
5:38
Bash
python solve.py
Fetching OSM network and hospitals... Calculating routes...
/work/solve.py:45: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.
hospitals_gdf['geometry'] = hospitals_gdf.geometry.centroid
Traceback (most recent call last):
File "/work/solve.py", line 174, in <module>
solve()
~~~~~^^
File "/work/solve.py", line 52, in solve
inc_nodes = ox.distance.nearest_nodes(G, incidents_gdf.geometry.x, incidents_gdf.geometry.y)
File "/usr/local/lib/python3.14/site-packages/osmnx/distance.py", line 376, in nearest_nodes
raise ImportError(msg)
ImportError: scikit-learn must be installed as an optional dependency to search an unprojected graph.
6:56
Edit
solve.py
Edited solve.py.
7:31
Edit
solve.py
old_string not found in solve.py
7:42
Read
solve.py
1 import osmnx as ox
2 import networkx as nx
3 import geopandas as gpd
4 import pandas as pd
5 from shapely.geometry import Point, LineString, MultiLineString, Polygon, MultiPolygon
6 import numpy as np
7
8 # Constants
9 BBOX = (48.83, 2.30, 48.88, 2.38) # min_lat, min_lon, max_lat, max_lon
10 CRS = "EPSG:2154" # RGF93 / Lambert-93 (Official metric CRS for France/Paris)
11 SPEED_DEFAULT = 30 / 3.6 # 30 km/h to m/s
12
13 def solve():
14 # 1. Load Incidents
15 incidents_df = pd.read_csv('incidents.csv')
16 # Create GeoDataFrame
17 incidents_gdf = gpd.GeoDataFrame(
18 incidents_df,
19 geometry=gpd.points_from_xy(incidents_df.longitude, incidents_df.latitude),
20 crs="EPSG:4326"
21 ).to_crs(CRS)
22
23 # 2. Fetch OSM Data (Network and Hospitals)
24 print("Fetching OSM network and hospitals...")
25 # Use osmnx to get the graph
26 # bbox: (left, bottom, right, top) -> (west, south, east, north)
27 west, south, east, north = BBOX[1], BBOX[0], BBOX[3], BBOX[2]
28 G = ox.graph_from_bbox(bbox=(west, south, east, north), network_type='drive')
29
30 # Apply speeds
31 G = ox.add_edge_speeds(G, fallback=30)
32 G = ox.add_edge_travel_times(G)
33
34 # Fetch Hospitals using Overpass (via osmnx or direct query if needed)
35 tags = {"amenity": "hospital"}
36 hospitals_gdf = ox.features_from_bbox(bbox=(west, south, east, north), tags=tags)
37
38 # Clean hospitals: ensure they are points/polygons and have names
39 hospitals_gdf = hospitals_gdf.copy()
40 if 'name' not in hospitals_gdf.columns:
41 hospitals_gdf['name'] = 'Unknown Hospital'
42 hospitals_gdf['hospital_name'] = hospitals_gdf['name'].fillna('Unknown Hospital')
43
44 # We want the centroid for routing if they are polygons
45 # Project to CRS before centroid for accuracy
46 hospitals_gdf = hospitals_gdf.to_crs(CRS)
47 hospitals_gdf['geometry'] = hospitals_gdf.geometry.centroid
48
49 # 3. Routing Analysis
50 print("Calculating routes...")
51
52 # Pre-calculate nearest nodes for incidents and hospitals
53 # Ensure we use projected coordinates for nearest_nodes if the graph is projected
54 # OSMnx graphs are usually in lat/lon unless specifically projected.
55 # Let's check the graph CRS.
56
57 if G.graph['crs'] is None:
58 # If graph is in lat/lon, we must provide lat/lon to nearest_nodes
59 inc_nodes = ox.distance.nearest_nodes(G, incidents_gdf.to_crs(G.graph['crs']).geometry.x, incidents_gdf.to_crs(G.graph['crs']).geometry.y)
60 else:
61 # If graph is projected, we should use projected coordinates
62 inc_nodes = ox.distance.nearest_nodes(G, incidents_gdf.geometry.x, incidents_gdf.geometry.y)
63
64 closest_routes_data = []
65 distance_matrix_data = []
66
67 for idx, inc in incidents_gdf.iterrows():
68 inc_id = inc['incident_id']
69 u = inc_nodes[idx]
70
71 # Calculate distances to all hospitals for this incident
72 hosp_distances = []
73 for h_idx, hosp in hospitals_gdf.iterrows():
74 v = hosp_nodes[h_idx]
75 try:
76 # weight='length' to get network distance
77 d = nx.shortest_path_length(G, u, v, weight='length')
78 hosp_distances.append({
79 'hospital_name': hosp['hospital_name'],
80 'network_distance_m': d,
81 'hosp_idx': h_idx,
82 'v_node': v
83 })
84 except nx.NetworkXNoPath:
85 continue
86
87 # Sort by distance
88 hosp_distances.sort(key=lambda x: x['network_distance_m'])
89
90 # Closest hospital for 'closest_hospital' layer
91 if hosp_distances:
92 closest = hosp_distances[0]
93 # Get path for the line geometry
94 path = nx.shortest_path(G, u, closest['v_node'], weight='length')
95 # Convert path to LineString
96 nodes_coords = []
97 for node in path:
98 nodes_coords.append((G.nodes[node]['x'], G.nodes[node]['y']))
99
100 route_line = LineString(nodes_coords)
101
102 closest_routes_data.append({
103 'incident_id': inc_id,
104 'hospital_name': closest['hospital_name'],
105 'network_distance_m': float(closest['network_distance_m']),
106 'geometry': route_line
107 })
108
109 # Distance matrix (top 3)
110 for rank, h_info in enumerate(hosp_distances[:3], 1):
111 distance_matrix_data.append({
112 'incident_id': inc_id,
113 'hospital_name': h_info['hospital_name'],
114 'rank': rank,
115 'network_distance_m': float(h_info['network_distance_m'])
116 })
117 else:
118 # Fallback if no hospital reachable
119 distance_matrix_data.append({
120 'incident_id': inc_id,
121 'hospital_name': 'None',
122 'rank': 1,
123 'network_distance_m': np.nan
124 })
125
126 # Convert to GeoDataFrames
127 closest_hospital_gdf = gpd.GeoDataFrame(closest_routes_data, crs=CRS)
128 distance_matrix_gdf = gpd.GeoDataFrame(distance_matrix_data, crs=CRS)
129
130 # 4. Isochrones (15 min)
131 print("Calculating isochrones...")
132 isochrone_data = []
133 # 15 minutes = 900 seconds
134 time_limit = 15 * 60
135
136 for idx, hosp in hospitals_gdf.iterrows():
137 h_name = hosp['hospital_name']
138 v = hosp_nodes[idx]
139
140 # Subgraph of reachable nodes within time limit
141 try:
142 subgraph_nodes = nx.single_source_dijkstra_path_length(G, v, cutoff=time_limit, weight='travel_time')
143 except (nx.NodeNotFound, nx.NetworkXError):
144 subgraph_nodes = {}
145
146 if subgraph_nodes:
147 # Get geometries of these nodes
148 node_points = [Point(G.nodes[n]['x'], G.nodes[n]['y']) for n in subgraph_nodes.keys()]
149 points_gdf = gpd.GeoDataFrame(geometry=node_points, crs=CRS)
150 try:
151 isochrone_poly = points_gdf.buffer(25).union_all()
152 except AttributeError:
153 isochrone_poly = points_gdf.buffer(25).unary_union
154
155 isochrone_data.append({
156 'hospital_name': h_name,
157 'travel_time_min': 15.0,
158 'geometry': isochrone_poly
159 })
160 else:
161 isochrone_data.append({
162 'hospital_name': h_name,
163 'travel_time_min': 0.0,
164 'geometry': Point(hosp.geometry.x, hosp.geometry.y).buffer(10)
165 })
166
167 isochrones_gdf = gpd.GeoDataFrame(isochrone_data, crs=CRS)
168
169 # 5. Save to GPKG
170 print("Saving to GPKG...")
171 # Prepare incidents layer
172 incidents_output = incidents_gdf[['incident_id', 'geometry']]
173
174 # Write layers
175 incidents_output.to_file('emergency_routing.gpkg', layer='incidents', driver='GPKG')
176 closest_hospital_gdf.to_file('emergency_routing.gpkg', layer='closest_hospital', driver='GPKG')
177 distance_matrix_gdf.to_file('emergency_routing.gpkg', layer='distance_matrix', driver='GPKG')
178 isochrones_gdf.to_file('emergency_routing.gpkg', layer='isochrones_15min', driver='GPKG')
179
180 print("Done!")
181
182 if __name__ == "__main__":
183 solve()
8:33
Edit
solve.py
Edited solve.py.
8:41
Bash
python solve.py
Fetching OSM network and hospitals... Calculating routes...
Traceback (most recent call last):
File "/work/solve.py", line 181, in <module>
solve()
~~~~~^^
File "/work/solve.py", line 56, in solve
inc_nodes = ox.distance.nearest_nodes(G_proj, incidents_gdf.geometry.x, incidents_gdf.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.
8:52
Bash
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 64.0 MB/s 0:00:00 Installing collected packages: scipy Successfully installed scipy-1.17.1
9:04
Bash
python solve.py
Fetching OSM network and hospitals... Calculating routes...
Traceback (most recent call last):
File "/work/solve.py", line 181, in <module>
solve()
~~~~~^^
File "/work/solve.py", line 72, in solve
v = hosp_nodes[h_idx]
~~~~~~~~~~^^^^^^^
IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices
9:15
Read
solve.py
1 import osmnx as ox
2 import networkx as nx
3 import geopandas as gpd
4 import pandas as pd
5 from shapely.geometry import Point, LineString, MultiLineString, Polygon, MultiPolygon
6 import numpy as np
7
8 # Constants
9 BBOX = (48.83, 2.30, 48.88, 2.38) # min_lat, min_lon, max_lat, max_lon
10 CRS = "EPSG:2154" # RGF93 / Lambert-93 (Official metric CRS for France/Paris)
11 SPEED_DEFAULT = 30 / 3.6 # 30 km/h to m/s
12
13 def solve():
14 # 1. Load Incidents
15 incidents_df = pd.read_csv('incidents.csv')
16 # Create GeoDataFrame
17 incidents_gdf = gpd.GeoDataFrame(
18 incidents_df,
19 geometry=gpd.points_from_xy(incidents_df.longitude, incidents_df.latitude),
20 crs="EPSG:4326"
21 ).to_crs(CRS)
22
23 # 2. Fetch OSM Data (Network and Hospitals)
24 print("Fetching OSM network and hospitals...")
25 # Use osmnx to get the graph
26 # bbox: (left, bottom, right, top) -> (west, south, east, north)
27 west, south, east, north = BBOX[1], BBOX[0], BBOX[3], BBOX[2]
28 G = ox.graph_from_bbox(bbox=(west, south, east, north), network_type='drive')
29
30 # Apply speeds
31 G = ox.add_edge_speeds(G, fallback=30)
32 G = ox.add_edge_travel_times(G)
33
34 # Fetch Hospitals using Overpass (via osmnx or direct query if needed)
35 tags = {"amenity": "hospital"}
36 hospitals_gdf = ox.features_from_bbox(bbox=(west, south, east, north), tags=tags)
37
38 # Clean hospitals: ensure they are points/polygons and have names
39 hospitals_gdf = hospitals_gdf.copy()
40 if 'name' not in hospitals_gdf.columns:
41 hospitals_gdf['name'] = 'Unknown Hospital'
42 hospitals_gdf['hospital_name'] = hospitals_gdf['name'].fillna('Unknown Hospital')
43
44 # We want the centroid for routing if they are polygons
45 # Project to CRS before centroid for accuracy
46 hospitals_gdf = hospitals_gdf.to_crs(CRS)
47 hospitals_gdf['geometry'] = hospitals_gdf.geometry.centroid
48
49 # 3. Routing Analysis
50 print("Calculating routes...")
51
52 # Project the graph to the target CRS for accurate distance calculations and to avoid scikit-learn dependency in nearest_nodes
53 G_proj = ox.project_graph(G, to_crs=CRS)
54
55 # Now we use projected coordinates for nearest_nodes which avoids the ImportError
56 inc_nodes = ox.distance.nearest_nodes(G_proj, incidents_gdf.geometry.x, incidents_gdf.geometry.y)
57 hosp_nodes = ox.distance.nearest_nodes(G_proj, hospitals_gdf.geometry.x, hospitals_gdf.geometry.y)
58
59 # We'll use the projected graph for all routing
60 G = G_proj
61
62 closest_routes_data = []
63 distance_matrix_data = []
64
65 for idx, inc in incidents_gdf.iterrows():
66 inc_id = inc['incident_id']
67 u = inc_nodes[idx]
68
69 # Calculate distances to all hospitals for this incident
70 hosp_distances = []
71 for h_idx, hosp in hospitals_gdf.iterrows():
72 v = hosp_nodes[h_idx]
73 try:
74 # weight='length' to get network distance
75 d = nx.shortest_path_length(G, u, v, weight='length')
76 hosp_distances.append({
77 'hospital_name': hosp['hospital_name'],
78 'network_distance_m': d,
79 'hosp_idx': h_idx,
80 'v_node': v
81 })
82 except (nx.NetworkXNoPath, nx.NodeNotFound):
83 continue
84
85 # Sort by distance
86 hosp_distances.sort(key=lambda x: x['network_distance_m'])
87
88 # Closest hospital for 'closest_hospital' layer
89 if hosp_distances:
90 closest = hosp_distances[0]
91 # Get path for the line geometry
92 path = nx.shortest_path(G, u, closest['v_node'], weight='length')
93 # Convert path to LineString
94 nodes_coords = []
95 for node in path:
96 nodes_coords.append((G.nodes[node]['x'], G.nodes[node]['y']))
97
98 route_line = LineString(nodes_coords)
99
100 closest_routes_data.append({
101 'incident_id': inc_id,
102 'hospital_name': closest['hospital_name'],
103 'network_distance_m': float(closest['network_distance_m']),
104 'geometry': route_line
105 })
106
107 # Distance matrix (top 3)
108 for rank, h_info in enumerate(hosp_distances[:3], 1):
109 distance_matrix_data.append({
110 'incident_id': inc_id,
111 'hospital_name': h_info['hospital_name'],
112 'rank': rank,
113 'network_distance_m': float(h_info['network_distance_m'])
114 })
115 else:
116 # Fallback if no hospital reachable
117 distance_matrix_data.append({
118 'incident_id': inc_id,
119 'hospital_name': 'None',
120 'rank': 1,
121 'network_distance_m': np.nan
122 })
123
124 # Convert to GeoDataFrames
125 closest_hospital_gdf = gpd.GeoDataFrame(closest_routes_data, crs=CRS)
126 distance_matrix_gdf = gpd.GeoDataFrame(distance_matrix_data, crs=CRS)
127
128 # 4. Isochrones (15 min)
129 print("Calculating isochrones...")
130 isochrone_data = []
131 # 15 minutes = 900 seconds
132 time_limit = 15 * 60
133
134 for idx, hosp in hospitals_gdf.iterrows():
135 h_name = hosp['hospital_name']
136 v = hosp_nodes[idx]
137
138 # Subgraph of reachable nodes within time limit
139 try:
140 subgraph_nodes = nx.single_source_dijkstra_path_length(G, v, cutoff=time_limit, weight='travel_time')
141 except (nx.NodeNotFound, nx.NetworkXError):
142 subgraph_nodes = {}
143
144 if subgraph_nodes:
145 # Get geometries of these nodes
146 node_points = [Point(G.nodes[n]['x'], G.nodes[n]['y']) for n in subgraph_nodes.keys()]
147 points_gdf = gpd.GeoDataFrame(geometry=node_points, crs=CRS)
148 try:
149 isochrone_poly = points_gdf.buffer(25).union_all()
150 except AttributeError:
151 isochrone_poly = points_gdf.buffer(25).unary_union
152
153 isochrone_data.append({
154 'hospital_name': h_name,
155 'travel_time_min': 15.0,
156 'geometry': isochrone_poly
157 })
158 else:
159 isochrone_data.append({
160 'hospital_name': h_name,
161 'travel_time_min': 0.0,
162 'geometry': Point(hosp.geometry.x, hosp.geometry.y).buffer(10)
163 })
164
165 isochrones_gdf = gpd.GeoDataFrame(isochrone_data, crs=CRS)
166
167 # 5. Save to GPKG
168 print("Saving to GPKG...")
169 # Prepare incidents layer
170 incidents_output = incidents_gdf[['incident_id', 'geometry']]
171
172 # Write layers
173 incidents_output.to_file('emergency_routing.gpkg', layer='incidents', driver='GPKG')
174 closest_hospital_gdf.to_file('emergency_routing.gpkg', layer='closest_hospital', driver='GPKG')
175 distance_matrix_gdf.to_file('emergency_routing.gpkg', layer='distance_matrix', driver='GPKG')
176 isochrones_gdf.to_file('emergency_routing.gpkg', layer='isochrones_15min', driver='GPKG')
177
178 print("Done!")
179
180 if __name__ == "__main__":
181 solve()
9:51
Edit
solve.py
old_string not found in solve.py
10:25
Write
solve.py
content · 185 lines · py
import osmnx as ox
import networkx as nx
import geopandas as gpd
import pandas as pd
from shapely.geometry import Point, LineString, MultiLineString, Polygon, MultiPolygon
import numpy as np
# Constants
BBOX = (48.83, 2.30, 48.88, 2.38) # min_lat, min_lon, max_lat, max_lon
CRS = "EPSG:2154" # RGF93 / Lambert-93 (Official metric CRS for France/Paris)
SPEED_DEFAULT = 30 / 3.6 # 30 km/h to m/s
def solve():
# 1. Load Incidents
incidents_df = pd.read_csv('incidents.csv')
# Create GeoDataFrame
incidents_gdf = gpd.GeoDataFrame(
incidents_df,
geometry=gpd.points_from_xy(incidents_df.longitude, incidents_df.latitude),
crs="EPSG:4326"
).to_crs(CRS)
# 2. Fetch OSM Data (Network and Hospitals)
print("Fetching OSM network and hospitals...")
# Use osmnx to get the graph
# bbox: (left, bottom, right, top) -> (west, south, east, north)
west, south, east, north = BBOX[1], BBOX[0], BBOX[3], BBOX[2]
G = ox.graph_from_bbox(bbox=(west, south, east, north), network_type='drive')
# Apply speeds
G = ox.add_edge_speeds(G, fallback=30)
G = ox.add_edge_travel_times(G)
# Fetch Hospitals using Overpass (via osmnx or direct query if needed)
tags = {"amenity": "hospital"}
hospitals_gdf = ox.features_from_bbox(bbox=(west, south, east, north), tags=tags)
# Clean hospitals: ensure they are points/polygons and have names
hospitals_gdf = hospitals_gdf.copy()
if 'name' not in hospitals_gdf.columns:
hospitals_gdf['name'] = 'Unknown Hospital'
hospitals_gdf['hospital_name'] = hospitals_gdf['name'].fillna('Unknown Hospital')
# We want the centroid for routing if they are polygons
hospitals_gdf = hospitals_gdf.to_crs(CRS)
hospitals_gdf['geometry'] = hospitals_gdf.geometry.centroid
# 3. Routing Analysis
print("Calculating routes...")
# Project the graph to the target CRS for accurate distance calculations and to avoid scikit-learn dependency in nearest_nodes
G_proj = ox.project_graph(G, to_crs=CRS)
# Now we use projected coordinates for nearest_nodes which avoids the ImportError
inc_nodes = ox.distance.nearest_nodes(G_proj, incidents_gdf.geometry.x, incidents_gdf.geometry.y)
hosp_nodes = ox.distance.nearest_nodes(G_proj, hospitals_gdf.geometry.x, hospitals_gdf.geometry.y)
# Map nodes to the appropriate index
inc_nodes_series = pd.Series(inc_nodes, index=incidents_gdf.index)
hosp_nodes_series = pd.Series(hosp_nodes, index=hospitals_gdf.index)
# We'll use the projected graph for all routing
G = G_proj
closest_routes_data = []
distance_matrix_data = []
for idx, inc in incidents_gdf.iterrows():
inc_id = inc['incident_id']
u = inc_nodes_series[idx]
# Calculate distances to all hospitals for this incident
hosp_distances = []
for h_idx, hosp in hospitals_gdf.iterrows():
v = hosp_nodes_series[h_idx]
try:
# weight='length' to get network distance
d = nx.shortest_path_length(G, u, v, weight='length')
hosp_distances.append({
'hospital_name': hosp['hospital_name'],
'network_distance_m': d,
'hosp_idx': h_idx,
'v_node': v
})
except (nx.NetworkXNoPath, nx.NodeNotFound):
continue
# Sort by distance
hosp_distances.sort(key=lambda x: x['network_distance_m'])
# Closest hospital for 'closest_hospital' layer
if hosp_distances:
closest = hosp_distances[0]
# Get path for the line geometry
path = nx.shortest_path(G, u, closest['v_node'], weight='length')
# Convert path to LineString
nodes_coords = []
for node in path:
nodes_coords.append((G.nodes[node]['x'], G.nodes[node]['y']))
route_line = LineString(nodes_coords)
closest_routes_data.append({
'incident_id': inc_id,
'hospital_name': closest['hospital_name'],
'network_distance_m': float(closest['network_distance_m']),
'geometry': route_line
})
# Distance matrix (top 3)
for rank, h_info in enumerate(hosp_distances[:3], 1):
distance_matrix_data.append({
'incident_id': inc_id,
'hospital_name': h_info['hospital_name'],
'rank': rank,
'network_distance_m': float(h_info['network_distance_m'])
})
else:
# Fallback if no hospital reachable
distance_matrix_data.append({
'incident_id': inc_id,
'hospital_name': 'None',
'rank': 1,
'network_distance_m': np.nan
})
# Convert to GeoDataFrames
closest_hospital_gdf = gpd.GeoDataFrame(closest_routes_data, crs=CRS)
distance_matrix_gdf = gpd.GeoDataFrame(distance_matrix_data, crs=CRS)
# 4. Isochrones (15 min)
print("Calculating isochrones...")
isochrone_data = []
# 15 minutes = 900 seconds
time_limit = 15 * 60
for idx, hosp in hospitals_gdf.iterrows():
h_name = hosp['hospital_name']
v = hosp_nodes_series[idx]
# Subgraph of reachable nodes within time limit
try:
subgraph_nodes = nx.single_source_dijkstra_path_length(G, v, cutoff=time_limit, weight='travel_time')
except (nx.NodeNotFound, nx.NetworkXError):
subgraph_nodes = {}
if subgraph_nodes:
# Get geometries of these nodes
node_points = [Point(G.nodes[n]['x'], G.nodes[n]['y']) for n in subgraph_nodes.keys()]
points_gdf = gpd.GeoDataFrame(geometry=node_points, crs=CRS)
try:
isochrone_poly = points_gdf.buffer(25).union_all()
except AttributeError:
isochrone_poly = points_gdf.buffer(25).unary_union
isochrone_data.append({
'hospital_name': h_name,
'travel_time_min': 15.0,
'geometry': isochrone_poly
})
else:
isochrone_data.append({
'hospital_name': h_name,
'travel_time_min': 0.0,
'geometry': Point(hosp.geometry.x, hosp.geometry.y).buffer(10)
})
isochrones_gdf = gpd.GeoDataFrame(isochrone_data, crs=CRS)
# 5. Save to GPKG
print("Saving to GPKG...")
# Prepare incidents layer
incidents_output = incidents_gdf[['incident_id', 'geometry']]
# Write layers
incidents_output.to_file('emergency_routing.gpkg', layer='incidents', driver='GPKG')
closest_hospital_gdf.to_file('emergency_routing.gpkg', layer='closest_hospital', driver='GPKG')
distance_matrix_gdf.to_file('emergency_routing.gpkg', layer='distance_matrix', driver='GPKG')
isochrones_gdf.to_file('emergency_routing.gpkg', layer='isochrones_15min', driver='GPKG')
print("Done!")
if __name__ == "__main__":
solve()
Wrote solve.py (7260 chars).
10:31
Bash
python solve.py
Fetching OSM network and hospitals... Calculating routes...
Traceback (most recent call last):
File "/work/solve.py", line 184, in <module>
solve()
~~~~~^^
File "/work/solve.py", line 129, in solve
distance_matrix_gdf = gpd.GeoDataFrame(distance_matrix_data, crs=CRS)
File "/usr/local/lib/python3.14/site-packages/geopandas/geodataframe.py", line 246, in __init__
raise ValueError(
...<3 lines>...
)
ValueError: Assigning CRS to a GeoDataFrame without a geometry column is not supported. Supply geometry using the 'geometry=' keyword argument, or by providing a DataFrame with column name 'geometry'
10:52
Edit
solve.py
Edited solve.py.
10:56
Bash
python solve.py
Fetching OSM network and hospitals... Calculating routes... Calculating isochrones... Saving to GPKG... Done!
11:25
Bash
ls -F
cache/ emergency_routing.gpkg incidents.csv solve.py
11:31
Bash
# Let's verify the layers and their columns in the GPKG
# Let's verify the layers and their columns in the GPKG
python -c "import geopandas as gpd; engine = 'pyogrio'; \
layers = gpd.io.fileio.list_layers('emergency_routing.gpkg'); \
print(layers); \
for layer in layers: \
gdf = gpd.read_file('emergency_routing.gpkg', layer=layer); \
print(f'Layer: {layer}, Columns: {gdf.columns.tolist()}, Geometry: {gdf.geometry.name}')"
File "<string>", line 1
import geopandas as gpd; engine = 'pyogrio'; layers = gpd.io.fileio.list_layers('emergency_routing.gpkg'); print(layers); for layer in layers: gdf = gpd.read_file('emergency_routing.gpkg', layer=layer); print(f'Layer: {layer}, Columns: {gdf.columns.tolist()}, Geometry: {gdf.geometry.name}')
^^^
SyntaxError: invalid syntax
11:37
Bash
python -c "import geopandas as gpd; \
python -c "import geopandas as gpd; \
import fiona; \
layers = fiona.listlayers('emergency_routing.gpkg'); \
print(layers); \
for layer in layers: \
gdf = gpd.read_file('emergency_routing.gpkg', layer=layer); \
print(f'Layer: {layer}, Columns: {gdf.columns.tolist()}, Geometry: {gdf.geometry.name}')"
File "<string>", line 1
import geopandas as gpd; import fiona; layers = fiona.listlayers('emergency_routing.gpkg'); print(layers); for layer in layers: gdf = gpd.read_file('emergency_routing.gpkg', layer=layer); print(f'Layer: {layer}, Columns: {gdf.columns.tolist()}, Geometry: {gdf.geometry.name}')
^^^
SyntaxError: invalid syntax
11:50
Bash
python -c "import geopandas as gpd; \
python -c "import geopandas as gpd; \
import fiona; \
layers = fiona.listlayers('emergency_routing.gpkg'); \
print(layers); \
for layer in layers: \
gdf = gpd.read_file('emergency_routing.gpkg', layer=layer); \
print(f'Layer: {layer}, Columns: {gdf.columns.tolist()}, Geometry: {gdf.geometry.name}')"
File "<string>", line 1
import geopandas as gpd; import fiona; layers = fiona.listlayers('emergency_routing.gpkg'); print(layers); for layer in layers: gdf = gpd.read_file('emergency_routing.gpkg', layer=layer); print(f'Layer: {layer}, Columns: {gdf.columns.tolist()}, Geometry: {gdf.geometry.name}')
^^^
SyntaxError: invalid syntax
12:04
Bash
python -c "import geopandas as gpd;
python -c "import geopandas as gpd;
import fiona;
layers = fiona.listlayers('emergency_routing.gpkg');
print(layers);
for layer in layers:
gdf = gpd.read_file('emergency_routing.gpkg', layer=layer);
print(f'Layer: {layer}, Columns: {gdf.columns.tolist()}, Geometry: {gdf.geometry.name}')"
['incidents', 'closest_hospital', 'distance_matrix', 'isochrones_15min'] Layer: incidents, Columns: ['incident_id', 'geometry'], Geometry: geometry Layer: closest_hospital, Columns: ['incident_id', 'hospital_name', 'network_distance_m', 'geometry'], Geometry: geometry Layer: distance_matrix, Columns: ['incident_id', 'hospital_name', 'rank', 'network_distance_m', 'geometry'], Geometry: geometry Layer: isochrones_15min, Columns: ['hospital_name', 'travel_time_min', 'geometry'], Geometry: geometry
13:06
```python
import osmnx as ox
import networkx as nx
import geopandas as gpd
import pandas as pd
from shapely.geometry import Point, LineString, MultiLineString, Polygon, MultiPolygon
import numpy as np
# Constants
BBOX = (48.83, 2.30, 48.88, 2.38) # min_lat, min_lon, max_lat, max_lon
CRS = "EPSG:2154" # RGF93 / Lambert-93 (Official metric CRS for France/Paris)
SPEED_DEFAULT = 30 / 3.6 # 30 km/h to m/s
def solve():
# 1. Load Incidents
incidents_df = pd.read_csv('incidents.csv')
# Create GeoDataFrame
incidents_gdf = gpd.GeoDataFrame(
incidents_df,
geometry=gpd.points_from_xy(incidents_df.longitude, incidents_df.latitude),
crs="EPSG:4326"
).to_crs(CRS)
# 2. Fetch OSM Data (Network and Hospitals)
print("Fetching OSM network and hospitals...")
# Use osmnx to get the graph
# bbox: (left, bottom, right, top) -> (west, south, east, north)
west, south, east, north = BBOX[1], BBOX[0], BBOX[3], BBOX[2]
G = ox.graph_from_bbox(bbox=(west, south, east, north), network_type='drive')
# Apply speeds
G = ox.add_edge_speeds(G, fallback=30)
G = ox.add_edge_travel_times(G)
# Fetch Hospitals using Overpass (via osmnx or direct query if needed)
tags = {"amenity": "hospital"}
hospitals_gdf = ox.features_from_bbox(bbox=(west, south, east, north), tags=tags)
# Clean hospitals: ensure they are points/polygons and have names
hospitals_gdf = hospitals_gdf.copy()
if 'name' not in hospitals_gdf.columns:
hospitals_gdf['name'] = 'Unknown Hospital'
hospitals_gdf['hospital_name'] = hospitals_gdf['name'].fillna('Unknown Hospital')
# We want the centroid for routing if they are polygons
hospitals_gdf = hospitals_gdf.to_crs(CRS)
hospitals_gdf['geometry'] = hospitals_gdf.geometry.centroid
# 3. Routing Analysis
print("Calculating routes...")
# Project the graph to the target CRS for accurate distance calculations and to avoid scikit-learn dependency in nearest_nodes
G_proj = ox.project_graph(G, to_crs=CRS)
# Now we use projected coordinates for nearest_nodes which avoids the ImportError
inc_nodes = ox.distance.nearest_nodes(G_proj, incidents_gdf.geometry.x, incidents_gdf.geometry.y)
hosp_nodes = ox.distance.nearest_nodes(G_proj, hospitals_gdf.geometry.x, hospitals_gdf.geometry.y)
# Map nodes to the appropriate index
inc_nodes_series = pd.Series(inc_nodes, index=incidents_gdf.index)
hosp_nodes_series = pd.Series(hosp_nodes, index=hospitals_gdf.index)
# We'll use the projected graph for all routing
G = G_proj
closest_routes_data = []
distance_matrix_data = []
for idx, inc in incidents_gdf.iterrows():
inc_id = inc['incident_id']
u = inc_nodes_series[idx]
# Calculate distances to all hospitals for this incident
hosp_distances = []
for h_idx, hosp in hospitals_gdf.iterrows():
v = hosp_nodes_series[h_idx]
try:
# weight='length' to get network distance
d = nx.shortest_path_length(G, u, v, weight='length')
hosp_distances.append({
'hospital_name': hosp['hospital_name'],
'network_distance_m': d,
'hosp_idx': h_idx,
'v_node': v
})
except (nx.NetworkXNoPath, nx.NodeNotFound):
continue
# Sort by distance
hosp_distances.sort(key=lambda x: x['network_distance_m'])
# Closest hospital for 'closest_hospital' layer
if hosp_distances:
closest = hosp_distances[0]
# Get path for the line geometry
path = nx.shortest_path(G, u, closest['v_node'], weight='length')
# Convert path to LineString
nodes_coords = []
for node in path:
nodes_coords.append((G.nodes[node]['x'], G.nodes[node]['y']))
route_line = LineString(nodes_coords)
closest_routes_data.append({
'incident_id': inc_id,
'hospital_name': closest['hospital_name'],
'network_distance_m': float(closest['network_distance_m']),
'geometry': route_line
})
# Distance matrix (top 3)
for rank, h_info in enumerate(hosp_distances[:3], 1):
distance_matrix_data.append({
'incident_id': inc_id,
'hospital_name': h_info['hospital_name'],
'rank': rank,
'network_distance_m': float(h_info['network_distance_m'])
})
else:
# Fallback if no hospital reachable
distance_matrix_data.append({
'incident_id': inc_id,
'hospital_name': 'None',
'rank': 1,
'network_distance_m': np.nan
})
# Convert to GeoDataFrames
closest_hospital_gdf = gpd.GeoDataFrame(closest_routes_data, crs=CRS)
# Distance matrix is tabular, but we must satisfy GeoDataFrame requirements if we want to save as layer
distance_matrix_df = pd.DataFrame(distance_matrix_data)
distance_matrix_gdf = gpd.GeoDataFrame(distance_matrix_df, geometry=[None]*len(distance_matrix_df), crs=CRS)
# 4. Isochrones (15 min)
print("Calculating isochrones...")
isochrone_data = []
# 15 minutes = 900 seconds
time_limit = 15 * 60
for idx, hosp in hospitals_gdf.iterrows():
h_name = hosp['hospital_name']
v = hosp_nodes_series[idx]
# Subgraph of reachable nodes within time limit
try:
subgraph_nodes = nx.single_source_dijkstra_path_length(G, v, cutoff=time_limit, weight='travel_time')
except (nx.NodeNotFound, nx.NetworkXError):
subgraph_nodes = {}
if subgraph_nodes:
# Get geometries of these nodes
node_points = [Point(G.nodes[n]['x'], G.nodes[n]['y']) for n in subgraph_nodes.keys()]
points_gdf = gpd.GeoDataFrame(geometry=node_points, crs=CRS)
try:
isochrone_poly = points_gdf.buffer(25).union_all()
except AttributeError:
isochrone_poly = points_gdf.buffer(25).unary_union
isochrone_data.append({
'hospital_name': h_name,
'travel_time_min': 15.0,
'geometry': isochrone_poly
})
else:
isochrone_data.append({
'hospital_name': h_name,
'travel_time_min': 0.0,
'geometry': Point(hosp.geometry.x, hosp.geometry.y).buffer(10)
})
isochrones_gdf = gpd.GeoDataFrame(isochrone_data, crs=CRS)
# 5. Save to GPKG
print("Saving to GPKG...")
# Prepare incidents layer
incidents_output = incidents_gdf[['incident_id', 'geometry']]
# Write layers
incidents_output.to_file('emergency_routing.gpkg', layer='incidents', driver='GPKG')
closest_hospital_gdf.to_file('emergency_routing.gpkg', layer='closest_hospital', driver='GPKG')
distance_matrix_gdf.to_file('emergency_routing.gpkg', layer='distance_matrix', driver='GPKG')
isochrones_gdf.to_file('emergency_routing.gpkg', layer='isochrones_15min', driver='GPKG')
print("Done!")
if __name__ == "__main__":
solve()
```