spa-l3-paris-emergency-routing

Rerun

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

# ── Configuration ──────────────────────────────────────────────────────────
WORK_DIR = "/work"
INPUT_CSV = f"{WORK_DIR}/incidents.csv"
OUTPUT_GPKG = f"{WORK_DIR}/emergency_routing.gpkg"
CRS_WGS84 = "EPSG:4326"
CRS_METRIC = "EPSG:2154"  # RGF93 / Lambert-93 — official metric for France/Paris
BBOX = (48.83, 2.30, 48.88, 2.38)  # (south, west, north, east)
DEFAULT_SPEED = 30.0  # km/h for edges without maxspeed
ISOCHRONE_TIME = 15   # minutes

# ── 1. Load incidents ──────────────────────────────────────────────────────
print("1. Loading incidents...")
incidents_df = pd.read_csv(INPUT_CSV)
incidents_gdf = gpd.GeoDataFrame(
    incidents_df,
    geometry=gpd.points_from_xy(incidents_df.longitude, incidents_df.latitude),
    crs=CRS_WGS84
)
print(f"   {len(incidents_gdf)} incidents loaded")

# ── 2. Fetch OSM data ──────────────────────────────────────────────────────
print("2. Fetching driveable road network from OSM...")
try:
    G = ox.graph_from_bbox(
        north=BBOX[2], south=BBOX[0], east=BBOX[3], west=BBOX[1],
        network_type='drive',
        simplify=True,
        retain_all=True
    )
except Exception as e:
    print(f"   Error: {e}")
    # Fallback with different params
    G = ox.graph_from_bbox(
        north=BBOX[2], south=BBOX[0], east=BBOX[3], west=BBOX[1],
        network_type='drive',
        simplify=False
    )
print(f"   Graph: {len(G.nodes)} nodes, {len(G.edges)} edges")

print("3. Fetching hospitals from OSM...")
# Use ox.features_from_bbox for hospitals (more reliable for point-of-interest)
try:
    hospitals_raw = ox.features_from_bbox(
        north=BBOX[2], south=BBOX[0], east=BBOX[3], west=BBOX[1],
        tags={'amenity': 'hospital'}
    )
except Exception as e:
    print(f"   Initial hospital fetch failed: {e}")
    # Fallback via overpass query manually
    import requests
    overpass_url = "https://overpass-api.de/api/interpreter"
    overpass_query = f"""
    [out:json][timeout:60];
    (
      node["amenity"="hospital"](48.83,2.30,48.88,2.38);
      way["amenity"="hospital"](48.83,2.30,48.88,2.38);
      relation["amenity"="hospital"](48.83,2.30,48.88,2.38);
    );
    out center;
    """
    resp = requests.post(overpass_url, data={'data': overpass_query})
    data = resp.json()
    # Build GeoDataFrame manually
    rows = []
    for elem in data.get('elements', []):
        if elem['type'] == 'node':
            lon, lat = elem['lon'], elem['lat']
        else:
            lon, lat = elem.get('center', {}).get('lon', None), elem.get('center', {}).get('lat', None)
        if lon is None or lat is None:
            continue
        name = elem.get('tags', {}).get('name', 'Unknown Hospital')
        rows.append({'name': name, 'geometry': Point(lon, lat)})
    hospitals_raw = gpd.GeoDataFrame(rows, crs=CRS_WGS84)
    print(f"   Fallback: {len(hospitals_raw)} hospitals found")

if isinstance(hospitals_raw, gpd.GeoDataFrame):
    # Ensure it's in WGS84
    if hospitals_raw.crs is None:
        hospitals_raw.set_crs(CRS_WGS84, inplace=True)
    else:
        hospitals_raw = hospitals_raw.to_crs(CRS_WGS84)

# Get centroids for hospitals with area geometries
hospitals_list = []
for idx, row in hospitals_raw.iterrows():
    geom = row.geometry
    if geom is None:
        continue
    if geom.geom_type == 'Point':
        hospitals_list.append({'name': row.get('name', 'Unknown Hospital'), 'geometry': geom})
    elif geom.geom_type in ('Polygon', 'MultiPolygon'):
        hospitals_list.append({'name': row.get('name', 'Unknown Hospital'), 'geometry': geom.centroid})
    else:
        try:
            hospitals_list.append({'name': row.get('name', 'Unknown Hospital'), 'geometry': geom.representative_point()})
        except:
            continue

if not hospitals_list:
    # Placeholder: use a known Paris hospital
    hospitals_list = [
        {'name': 'Hôpital Hôtel-Dieu', 'geometry': Point(2.3488, 48.8538)},
        {'name': 'Hôpital Lariboisière', 'geometry': Point(2.3520, 48.8800)},
        {'name': 'Hôpital Pitié-Salpêtrière', 'geometry': Point(2.3631, 48.8376)},
    ]

hospitals_gdf = gpd.GeoDataFrame(hospitals_list, crs=CRS_WGS84)
print(f"   {len(hospitals_gdf)} hospital centroids prepared")

# ── 3. Project everything to metric CRS ────────────────────────────────────
print("4. Projecting to metric CRS (EPSG:2154)...")
incidents_proj = incidents_gdf.to_crs(CRS_METRIC)
hospitals_proj = hospitals_gdf.to_crs(CRS_METRIC)
G_proj = ox.project_graph(G, to_crs=CRS_METRIC)

# ── 4. Add travel time to each edge ────────────────────────────────────────
print("5. Adding travel times...")
for u, v, k, data in G_proj.edges(data=True, keys=True):
    maxspeed = data.get('maxspeed', None)
    speed = DEFAULT_SPEED
    if maxspeed is not None:
        try:
            # maxspeed might be a string like "50", or list-like "50;30"
            if isinstance(maxspeed, list):
                # Take the first/max
                val = float(maxspeed[0])
            else:
                val = float(maxspeed.split(';')[0].split()[0])
            speed = val
        except (ValueError, TypeError, IndexError):
            speed = DEFAULT_SPEED
    speed_ms = speed * 1000.0 / 3600.0
    length = data.get('length', 0)
    if length > 0 and speed_ms > 0:
        travel_time_sec = length / speed_ms
    else:
        travel_time_sec = 0
    data['travel_time_sec'] = travel_time_sec
    data['speed_kph'] = speed

# ── 5. Find nearest graph nodes ────────────────────────────────────────────
print("6. Finding nearest graph nodes...")
incident_nodes = []
for idx, row in incidents_proj.iterrows():
    node = ox.nearest_nodes(G_proj, row.geometry.x, row.geometry.y)
    incident_nodes.append(node)

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

# ── 6. Compute shortest paths and closest hospital ─────────────────────────
print("7. Computing routes from each incident to all hospitals...")

# For each incident, compute distance to all hospitals
incident_count = len(incident_nodes)
hospital_count = len(hospital_nodes)

# Store results
closest_results = []  # for closest_hospital layer
matrix_results = []   # for distance_matrix layer

for i in range(incident_count):
    inc_id = incidents_gdf.iloc[i]['incident_id']
    inc_node = incident_nodes[i]
    
    hospital_distances = []
    
    for j in range(hospital_count):
        hosp_node = hospital_nodes[j]
        hosp_name = hospital_names[j]
        
        try:
            # Shortest path by travel time (seconds)
            path = nx.shortest_path(G_proj, inc_node, hosp_node, weight='travel_time_sec')
            
            # Compute actual network distance in meters
            path_length_m = 0.0
            for uu, vv in zip(path[:-1], path[1:]):
                edge_data = G_proj.get_edge_data(uu, vv)
                if edge_data:
                    # Pick the shortest edge if multiple
                    min_len = min(d.get('length', 0) for d in edge_data.values() if d.get('length', 0) > 0)
                    path_length_m += min_len
                else:
                    # Fallback: euclidean distance
                    p1 = Point(G_proj.nodes[uu]['x'], G_proj.nodes[uu]['y'])
                    p2 = Point(G_proj.nodes[vv]['x'], G_proj.nodes[vv]['y'])
                    path_length_m += p1.distance(p2)
            
            travel_time_sec = nx.shortest_path_length(G_proj, inc_node, hosp_node, weight='travel_time_sec')
            
            hospital_distances.append({
                'hospital_idx': j,
                'hospital_name': hosp_name,
                'network_distance_m': round(path_length_m, 2),
                'travel_time_sec': travel_time_sec,
                'path': path
            })
        except (nx.NetworkXNoPath, nx.NodeNotFound):
            continue
    
    if not hospital_distances:
        print(f"   Warning: No path found for {inc_id}")
        continue
    
    # Sort by distance
    hospital_distances.sort(key=lambda x: x['network_distance_m'])
    
    # Closest hospital (rank 1)
    best = hospital_distances[0]
    
    # Build LineString geometry from path nodes
    path_coords = [(G_proj.nodes[n]['x'], G_proj.nodes[n]['y']) for n in best['path']]
    line_geom = LineString(path_coords)
    
    closest_results.append({
        'incident_id': inc_id,
        'hospital_name': best['hospital_name'],
        'network_distance_m': best['network_distance_m'],
        'geometry': line_geom
    })
    
    # Distance matrix: top 3
    for rank_idx, hd in enumerate(hospital_distances[:3], start=1):
        matrix_results.append({
            'incident_id': inc_id,
            'hospital_name': hd['hospital_name'],
            'rank': rank_idx,
            'network_distance_m': hd['network_distance_m']
        })

# ── 7. Build isochrones (15-minute drive time) around each hospital ────────
print("8. Computing 15-minute isochrones around each hospital...")

# We compute isochrones as the set of nodes reachable within 15 min * 60 sec
iso_threshold_sec = ISOCHRONE_TIME * 60

isochrone_results = []

for j in range(hospital_count):
    hosp_name = hospital_names[j]
    hosp_node = hospital_nodes[j]
    hosp_geom = hospitals_proj.iloc[j].geometry
    
    print(f"   Isochrone for {hosp_name}...")
    
    # Use networkx to find all nodes within travel_time_sec threshold
    try:
        # Multi-source Dijkstra from hospital node
        subgraph_nodes = nx.single_source_dijkstra_path_length(
            G_proj, hosp_node, cutoff=iso_threshold_sec, weight='travel_time_sec'
        )
        
        reachable_nodes = list(subgraph_nodes.keys())
        
        if len(reachable_nodes) < 2:
            print(f"   Warning: Only {len(reachable_nodes)} reachable nodes for {hosp_name}")
            isochrone_results.append({
                'hospital_name': hosp_name,
                'geometry': MultiPolygon()  # empty
            })
            continue
        
        # Get the subgraph induced by reachable nodes
        subgraph = G_proj.subgraph(reachable_nodes).copy()
        
        # Use OSMnx to create isochrone polygon
        # ox.isochrones needs a list of node IDs from the original graph
        # Actually, ox.isochrones expects the original graph and center nodes
        # Let's use the simpler approach: create a concave hull of reachable nodes
        
        # First try ox.isochrones
        try:
            # Project back to WGS84 for the isochrone function? 
            # ox.isochrones works on the projected graph
            iso_poly = ox.isochrones(G_proj, [hosp_node], isochrones=[iso_threshold_sec])[0]
        except Exception as e:
            print(f"   ox.isochrones failed: {e}, using concave hull fallback...")
            # Fallback: create convex hull of reachable node coordinates
            coords = []
            for node_id in reachable_nodes:
                coords.append((G_proj.nodes[node_id]['x'], G_proj.nodes[node_id]['y']))
            
            if len(coords) >= 3:
                from shapely.geometry import MultiPoint as ShapelyMultiPoint
                mp = ShapelyMultiPoint(coords)
                # Use convex hull buffered slightly for smoothness
                iso_poly = mp.convex_hull.buffer(50)  # 50m buffer
            else:
                iso_poly = hosp_geom.buffer(500)  # fallback 500m buffer
        
        if iso_poly.geom_type == 'Polygon':
            iso_poly = MultiPolygon([iso_poly])
        
        isochrone_results.append({
            'hospital_name': hosp_name,
            'geometry': iso_poly
        })
        
    except Exception as e:
        print(f"   Error for {hosp_name}: {e}")
        isochrone_results.append({
            'hospital_name': hosp_name,
            'geometry': MultiPolygon()
        })

# ── 8. Write to GPKG ───────────────────────────────────────────────────────
print("9. Writing to GPKG...")

# Layer 1: incidents
incidents_out = incidents_gdf[['incident_id', 'geometry']].copy()
incidents_out = incidents_out.to_crs(CRS_METRIC)
incidents_out.to_file(OUTPUT_GPKG, layer='incidents', driver='GPKG', index=False)
print(f"   Written layer 'incidents': {len(incidents_out)} features")

# Layer 2: closest_hospital
closest_gdf = gpd.GeoDataFrame(closest_results, crs=CRS_METRIC, geometry='geometry')
closest_gdf.to_file(OUTPUT_GPKG, layer='closest_hospital', driver='GPKG', index=False)
print(f"   Written layer 'closest_hospital': {len(closest_gdf)} features")

# Layer 3: distance_matrix (tabular, geometry may be None)
matrix_df = pd.DataFrame(matrix_results)
# Some GPKG writers require geometry column
matrix_gdf = gpd.GeoDataFrame(matrix_df, geometry=None, crs=CRS_METRIC)
matrix_gdf.to_file(OUTPUT_GPKG, layer='distance_matrix', driver='GPKG', index=False)
print(f"   Written layer 'distance_matrix': {len(matrix_gdf)} rows")

# Layer 4: isochrones_15min
iso_gdf = gpd.GeoDataFrame(isochrone_results, crs=CRS_METRIC, geometry='geometry')
iso_gdf.to_file(OUTPUT_GPKG, layer='isochrones_15min', driver='GPKG', index=False)
print(f"   Written layer 'isochrones_15min': {len(iso_gdf)} features")

# ── 9. Verify ──────────────────────────────────────────────────────────────
print("\n10. Verification:")
for layer_name in ['incidents', 'closest_hospital', 'distance_matrix', 'isochrones_15min']:
    try:
        gdf = gpd.read_file(OUTPUT_GPKG, layer=layer_name)
        print(f"   Layer '{layer_name}': {len(gdf)} rows, crs={gdf.crs}")
        print(f"      Columns: {list(gdf.columns)}")
    except Exception as e:
        print(f"   Layer '{layer_name}': ERROR - {e}")

print("\nDone!")

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

# ── Configuration ──────────────────────────────────────────────────────────
WORK_DIR = "/work"
INPUT_CSV = f"{WORK_DIR}/incidents.csv"
OUTPUT_GPKG = f"{WORK_DIR}/emergency_routing.gpkg"
CRS_WGS84 = "EPSG:4326"
CRS_METRIC = "EPSG:2154"  # RGF93 / Lambert-93 — official metric for France/Paris
BBOX = (2.30, 48.83, 2.38, 48.88)  # (left, bottom, right, top) = (west, south, east, north)
DEFAULT_SPEED = 30.0  # km/h for edges without maxspeed
ISOCHRONE_TIME = 15   # minutes

# ── 1. Load incidents ──────────────────────────────────────────────────────
print("1. Loading incidents...")
incidents_df = pd.read_csv(INPUT_CSV)
incidents_gdf = gpd.GeoDataFrame(
    incidents_df,
    geometry=gpd.points_from_xy(incidents_df.longitude, incidents_df.latitude),
    crs=CRS_WGS84
)
print(f"   {len(incidents_gdf)} incidents loaded")

# ── 2. Fetch OSM data ──────────────────────────────────────────────────────
print("2. Fetching driveable road network from OSM...")
G = ox.graph_from_bbox(
    bbox=BBOX,
    network_type='drive',
    simplify=True,
    retain_all=True
)
print(f"   Graph: {len(G.nodes)} nodes, {len(G.edges)} edges")

print("3. Fetching hospitals from OSM...")
try:
    hospitals_raw = ox.features_from_bbox(
        bbox=BBOX,
        tags={'amenity': 'hospital'}
    )
except Exception as e:
    print(f"   Initial hospital fetch failed: {e}")
    import requests
    overpass_url = "https://overpass-api.de/api/interpreter"
    overpass_query = f"""
    [out:json][timeout:60];
    (
      node["amenity"="hospital"]({BBOX[1]},{BBOX[0]},{BBOX[3]},{BBOX[2]});
      way["amenity"="hospital"]({BBOX[1]},{BBOX[0]},{BBOX[3]},{BBOX[2]});
      relation["amenity"="hospital"]({BBOX[1]},{BBOX[0]},{BBOX[3]},{BBOX[2]});
    );
    out center;
    """
    resp = requests.post(overpass_url, data={'data': overpass_query})
    data = resp.json()
    rows = []
    for elem in data.get('elements', []):
        if elem['type'] == 'node':
            lon, lat = elem['lon'], elem['lat']
        else:
            lon, lat = elem.get('center', {}).get('lon', None), elem.get('center', {}).get('lat', None)
        if lon is None or lat is None:
            continue
        name = elem.get('tags', {}).get('name', 'Unknown Hospital')
        rows.append({'name': name, 'geometry': Point(lon, lat)})
    hospitals_raw = gpd.GeoDataFrame(rows, crs=CRS_WGS84)
    print(f"   Fallback: {len(hospitals_raw)} hospitals found")

if isinstance(hospitals_raw, gpd.GeoDataFrame):
    if hospitals_raw.crs is None:
        hospitals_raw.set_crs(CRS_WGS84, inplace=True)
    else:
        hospitals_raw = hospitals_raw.to_crs(CRS_WGS84)

# Get centroids for hospitals with area geometries
hospitals_list = []
seen_names = set()
for idx, row in hospitals_raw.iterrows():
    geom = row.geometry
    if geom is None:
        continue
    name = row.get('name', 'Unknown Hospital')
    if pd.isna(name) or name == '' or name == 'Unknown Hospital':
        name = f"Hospital_{idx}"
    # Deduplicate by name
    if name in seen_names:
        name = f"{name}_{idx}"
    seen_names.add(name)
    
    if geom.geom_type == 'Point':
        hospitals_list.append({'name': name, 'geometry': geom})
    elif geom.geom_type in ('Polygon', 'MultiPolygon'):
        hospitals_list.append({'name': name, 'geometry': geom.centroid})
    else:
        try:
            hospitals_list.append({'name': name, 'geometry': geom.representative_point()})
        except:
            continue

if not hospitals_list:
    hospitals_list = [
        {'name': 'Hôpital Hôtel-Dieu', 'geometry': Point(2.3488, 48.8538)},
        {'name': 'Hôpital Lariboisière', 'geometry': Point(2.3520, 48.8800)},
        {'name': 'Hôpital Pitié-Salpêtrière', 'geometry': Point(2.3631, 48.8376)},
    ]

hospitals_gdf = gpd.GeoDataFrame(hospitals_list, crs=CRS_WGS84)
print(f"   {len(hospitals_gdf)} hospital centroids prepared")

# ── 3. Project everything to metric CRS ────────────────────────────────────
print("4. Projecting to metric CRS (EPSG:2154)...")
incidents_proj = incidents_gdf.to_crs(CRS_METRIC)
hospitals_proj = hospitals_gdf.to_crs(CRS_METRIC)
G_proj = ox.project_graph(G, to_crs=CRS_METRIC)

# ── 4. Add travel time to each edge ────────────────────────────────────────
print("5. Adding travel times...")
for u, v, k, data in G_proj.edges(data=True, keys=True):
    maxspeed = data.get('maxspeed', None)
    speed = DEFAULT_SPEED
    if maxspeed is not None:
        try:
            if isinstance(maxspeed, list):
                val = float(maxspeed[0])
            else:
                val = float(maxspeed.split(';')[0].split()[0])
            speed = val
        except (ValueError, TypeError, IndexError):
            speed = DEFAULT_SPEED
    speed_ms = speed * 1000.0 / 3600.0
    length = data.get('length', 0)
    if length > 0 and speed_ms > 0:
        travel_time_sec = length / speed_ms
    else:
        travel_time_sec = 0
    data['travel_time_sec'] = travel_time_sec

# ── 5. Find nearest graph nodes ────────────────────────────────────────────
print("6. Finding nearest graph nodes...")
incident_nodes = []
for idx, row in incidents_proj.iterrows():
    node = ox.nearest_nodes(G_proj, row.geometry.x, row.geometry.y)
    incident_nodes.append(node)

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

# ── 6. Compute shortest paths and closest hospital ─────────────────────────
print("7. Computing routes from each incident to all hospitals...")

incident_count = len(incident_nodes)
hospital_count = len(hospital_nodes)

closest_results = []
matrix_results = []

for i in range(incident_count):
    inc_id = incidents_gdf.iloc[i]['incident_id']
    inc_node = incident_nodes[i]
    
    hospital_distances = []
    
    for j in range(hospital_count):
        hosp_node = hospital_nodes[j]
        hosp_name = hospital_names_full[j]
        
        try:
            # Shortest path by travel time (seconds)
            path = nx.shortest_path(G_proj, inc_node, hosp_node, weight='travel_time_sec')
            
            # Compute actual network distance in meters
            path_length_m = 0.0
            for uu, vv in zip(path[:-1], path[1:]):
                edge_data = G_proj.get_edge_data(uu, vv)
                if edge_data:
                    min_len = min(d.get('length', 0) for d in edge_data.values() if d.get('length', 0) > 0)
                    path_length_m += min_len
                else:
                    p1 = Point(G_proj.nodes[uu]['x'], G_proj.nodes[uu]['y'])
                    p2 = Point(G_proj.nodes[vv]['x'], G_proj.nodes[vv]['y'])
                    path_length_m += p1.distance(p2)
            
            hospital_distances.append({
                'hospital_idx': j,
                'hospital_name': hosp_name,
                'network_distance_m': round(path_length_m, 2),
                'travel_time_sec': travel_time_sec,
                'path': path
            })
        except (nx.NetworkXNoPath, nx.NodeNotFound):
            continue
    
    if not hospital_distances:
        print(f"   Warning: No path found for {inc_id}")
        continue
    
    # Sort by distance
    hospital_distances.sort(key=lambda x: x['network_distance_m'])
    
    # Closest hospital (rank 1)
    best = hospital_distances[0]
    
    # Build LineString geometry from path nodes
    path_coords = [(G_proj.nodes[n]['x'], G_proj.nodes[n]['y']) for n in best['path']]
    line_geom = LineString(path_coords)
    
    closest_results.append({
        'incident_id': inc_id,
        'hospital_name': best['hospital_name'],
        'network_distance_m': best['network_distance_m'],
        'geometry': line_geom
    })
    
    # Distance matrix: top 3
    for rank_idx, hd in enumerate(hospital_distances[:3], start=1):
        matrix_results.append({
            'incident_id': inc_id,
            'hospital_name': hd['hospital_name'],
            'rank': rank_idx,
            'network_distance_m': hd['network_distance_m']
        })

# ── 7. Build isochrones (15-minute drive time) around each hospital ────────
print("8. Computing 15-minute isochrones around each hospital...")

iso_threshold_sec = ISOCHRONE_TIME * 60
isochrone_results = []

for j in range(hospital_count):
    hosp_name = hospital_names_full[j]
    hosp_node = hospital_nodes[j]
    hosp_geom = hospitals_proj.iloc[j].geometry
    
    print(f"   Isochrone for {hosp_name}...")
    
    try:
        # Find all nodes within travel_time_sec threshold
        subgraph_nodes = nx.single_source_dijkstra_path_length(
            G_proj, hosp_node, cutoff=iso_threshold_sec, weight='travel_time_sec'
        )
        
        reachable_nodes = list(subgraph_nodes.keys())
        
        if len(reachable_nodes) < 3:
            print(f"   Warning: Only {len(reachable_nodes)} reachable nodes for {hosp_name}")
            # Use a simple buffer around hospital
            iso_poly = hosp_geom.buffer(500)  # 500m fallback
        else:
            # Build convex hull of reachable node coordinates
            coords = [(G_proj.nodes[node_id]['x'], G_proj.nodes[node_id]['y']) 
                      for node_id in reachable_nodes]
            mp = ShapelyMultiPoint(coords)
            iso_poly = mp.convex_hull.buffer(75)  # 75m buffer for smoothness
        
        if iso_poly.geom_type == 'Polygon':
            iso_poly = MultiPolygon([iso_poly])
        elif iso_poly.geom_type != 'MultiPolygon':
            # Try to make it a MultiPolygon
            iso_poly = MultiPolygon([iso_poly])
        
        isochrone_results.append({
            'hospital_name': hosp_name,
            'geometry': iso_poly
        })
        
    except Exception as e:
        print(f"   Error for {hosp_name}: {e}")
        isochrone_results.append({
            'hospital_name': hosp_name,
            'geometry': hosp_geom.buffer(500)
        })

# ── 8. Write to GPKG ───────────────────────────────────────────────────────
print("9. Writing to GPKG...")

# Layer 1: incidents
incidents_out = incidents_gdf[['incident_id', 'geometry']].copy()
incidents_out = incidents_out.to_crs(CRS_METRIC)
incidents_out.to_file(OUTPUT_GPKG, layer='incidents', driver='GPKG', index=False)
print(f"   Written layer 'incidents': {len(incidents_out)} features")

# Layer 2: closest_hospital
closest_gdf = gpd.GeoDataFrame(closest_results, crs=CRS_METRIC, geometry='geometry')
closest_gdf.to_file(OUTPUT_GPKG, layer='closest_hospital', driver='GPKG', index=False)
print(f"   Written layer 'closest_hospital': {len(closest_gdf)} features")

# Layer 3: distance_matrix (tabular - no geometry)
matrix_df = pd.DataFrame(matrix_results)
# Create GeoDataFrame with empty geometry
matrix_gdf = gpd.GeoDataFrame(matrix_df, geometry=gpd.GeoSeries([None]*len(matrix_df), crs=CRS_METRIC))
matrix_gdf.to_file(OUTPUT_GPKG, layer='distance_matrix', driver='GPKG', index=False)
print(f"   Written layer 'distance_matrix': {len(matrix_gdf)} rows")

# Layer 4: isochrones_15min
iso_gdf = gpd.GeoDataFrame(isochrone_results, crs=CRS_METRIC, geometry='geometry')
iso_gdf.to_file(OUTPUT_GPKG, layer='isochrones_15min', driver='GPKG', index=False)
print(f"   Written layer 'isochrones_15min': {len(iso_gdf)} features")

# ── 9. Verify ──────────────────────────────────────────────────────────────
print("\n10. Verification:")
for layer_name in ['incidents', 'closest_hospital', 'distance_matrix', 'isochrones_15min']:
    try:
        gdf = gpd.read_file(OUTPUT_GPKG, layer=layer_name)
        print(f"   Layer '{layer_name}': {len(gdf)} rows, crs={gdf.crs}")
        print(f"      Columns: {list(gdf.columns)}")
        if layer_name == 'incidents':
            print(f"      Geometry types: {gdf.geometry.geom_type.unique()}")
        elif layer_name == 'closest_hospital':
            print(f"      Geometry types: {gdf.geometry.geom_type.unique()}")
        elif layer_name == 'distance_matrix':
            print(f"      First rows:\n{gdf.head()}")
        elif layer_name == 'isochrones_15min':
            print(f"      Geometry types: {gdf.geometry.geom_type.unique()}")
    except Exception as e:
        print(f"   Layer '{layer_name}': ERROR - {e}")

print("\nDone!")