spa-l1-capetown-hospital-nn

Rerun

import geopandas as gpd
import pandas as pd
from shapely.ops import nearest_points

def solve():
    # Load data
    addresses = gpd.read_parquet('addresses.parquet')
    hospitals = gpd.read_parquet('hospitals.parquet')

    print(f"Addresses: {len(addresses)}")
    print(f"Hospitals: {len(hospitals)}")

    # Check CRS
    if addresses.crs is None or hospitals.crs is None:
        raise ValueError("CRS not found in input files")
    
    if addresses.crs != hospitals.crs:
        print(f"Reprojecting hospitals from {hospitals.crs} to {addresses.crs}")
        hospitals = hospitals.to_crs(addresses.crs)

    # To calculate distance in metres, we need a projected CRS.
    # Let's assume the input might be WGS84 (EPSG:4326).
    # We should project to a metric CRS if it's geographic.
    
    # Check if CRS is geographic
    if addresses.crs.is_geographic:
        print("Projecting to UTM for distance calculations...")
        # Use a suitable UTM zone for Western Cape (likely UTM 34S or 35S)
        # Or just use a general metric CRS like EPSG:3857 for an approximation, 
        # but for better accuracy in Western Cape, let's use a local UTM or similar.
        # Actually, a more robust way is to use a local UTM or a local projection.
        # Let's try to find a suitable UTM zone automatically or just use 3857 for simplicity 
        # if accuracy isn't ultra-critical, but since it's for EMS, accuracy matters.
        # Western Cape is roughly between 18E and 25E and 33S and 35S.
        # UTM zone 34S covers 18E to 24E. UTM zone 35S covers 24E to 30E.
        # Let's use a generic approach: project to a metric CRS.
        target_crs = 'EPSG:3857' # Web Mercator is metric but distorts.
        # Better: find the best UTM zone.
        # Let's just use a simple projection if it's 4326.
        # For Western Cape, EPSG:2056 (South African Lo system) is very accurate.
        target_crs = 'EPSG:2056'
        print(f"Using {target_crs} for distance calculation.")
        addresses_proj = addresses.to_crs(target_crs)
        hospitals_proj = hospitals.to_crs(target_crs)
    else:
        addresses_proj = addresses
        hospitals_proj = hospitals
        target_crs = addresses.crs

    # To find the nearest hospital for each address efficiently
    # We can use a spatial index (sindex)
    
    # Prepare results list
    nearest_hospital_names = []
    distances = []

    # Get spatial index for hospitals
    h_sindex = hospitals_proj.sindex

    for idx, addr in addresses_proj.iterrows():
        # Find nearest hospital
        # sindex.nearest returns indices of the nearest geometry
        # Note: sindex.nearest returns a tuple of (idx_of_addresses, idx_of_hospitals)
        # Since we are iterating over addresses one by one, we just need the hospital index.
        
        # For a single geometry, nearest() is efficient
        # The sindex.nearest approach for a single point:
        nearest_idx_list = list(h_sindex.nearest(addr.geometry, return_distance=False))
        # nearest_idx_list is often a list of arrays if querying multiple, 
        # but here we query one. Let's check documentation or use a different approach.
        
        # Standard way for single point using sindex:
        # 1. Find candidates using bounding box or tree
        # 2. Calculate exact distances
        
        # Using the more robust approach for single geometry:
        # We'll use the spatial index to find the closest one
        # In recent geopandas/pygeos, sindex.nearest is very good.
        
        # Let's try the simplest approach first:
        # For each address, find the nearest hospital
        # Given the likely size of the datasets, we'll use a spatial index.
        
        # For single point 'addr.geometry':
        # nearest_idx = h_sindex.nearest(addr.geometry)[1][0] # This is for bulk.
        # Let's use a more manual approach if sindex.nearest is tricky for single items.
        
        # Let's use the bulk approach for speed
        pass

    # Redoing the logic for efficiency:
    # Use bulk nearest neighbor search
    # This requires geopandas >= 0.10.0 or having pygeos/shapely 2.0
    
    # We will perform a bulk search using sindex.nearest
    # This returns indices of the nearest geometries in hospitals_proj for each geometry in addresses_proj
    
    # addresses_proj.sindex.nearest(hospitals_proj.geometry) is for many-to-many.
    # We want one-to-many (each address -> one hospital).
    
    # Use sindex.query or similar? 
    # In newer geopandas:
    # nearest_indices = addresses_proj.sindex.query(hospitals_proj.geometry, predicate='intersects') -> No.
    
    # Let's use a loop with sindex.query if we want to be safe, but actually
    # addresses_proj.sindex.query(hospitals_proj.geometry, method='nearest')
    # Or simply:
    
    from scipy.spatial import cKDTree
    
    # Extract coordinates for KDTree (only works for Points)
    # If addresses are not points, we'll use the geometry.centroid
    addr_coords = [(p.x, p.y) for p in addresses_proj.geometry.centroid]
    hosp_coords = [(p.x, p.y) for p in hospitals_proj.geometry.centroid]
    
    tree = cKDTree(hosp_coords)
    dist, hosp_indices = tree.query(addr_coords)

    # Mapping indices back to names and distances
    # dist is in the units of the projection (metres)
    
    results_df = pd.DataFrame({
        'address_id': addresses['address_id'],
        'nearest_hospital_name': hospitals.iloc[hosp_indices]['name'].values, # assuming 'name' column exists
        'distance_m': dist
    })
    
    # Wait, I should check the column name for hospital name in hospitals.parquet
    # Let's inspect columns first.

solve()

import geopandas as gpd
import pandas as pd
import numpy as np
from scipy.spatial import cKDTree

def solve():
    # Load data
    addresses = gpd.read_parquet('addresses.parquet')
    hospitals = gpd.read_parquet('hospitals.parquet')

    print(f"Addresses: {len(addresses)}")
    print(f"Hospitals: {len(hospitals)}")

    if addresses.crs is None or hospitals.crs is None:
        raise ValueError("CRS not found in input files")
    
    # Use a projected CRS for accurate distance calculation in metres
    # For Western Cape, South African Lo (EPSG:2056) is a good choice.
    # If it is already projected, we'll use it.
    
    if addresses.crs.is_geographic:
        print("Projecting to EPSG:2056 for distance calculation...")
        target_crs = 'EPSG:2056'
    else:
        target_crs = addresses.crs
        print(f"Using existing CRS: {target_crs}")

    # Project both to the target metric CRS
    addresses_proj = addresses.to_crs(target_crs)
    hospitals_proj = hospitals.to_crs(target_crs)

    # Use cKDTree for efficient nearest neighbor search
    # Extract coordinates (assuming points)
    # If they are not points, we'll use centroids
    addr_coords = np.array([(p.x, p.y) for p in addresses_proj.geometry.centroid])
    hosp_coords = np.array([(p.x, p.y) for p in hospitals_proj.geometry.centroid])
    
    tree = cKDTree(hosp_coords)
    dist, hosp_indices = tree.query(addr_coords)

    # Prepare the output GeoDataFrame
    # We need to keep original address_id and geometry from the original 'addresses'
    
    # Create a new GeoDataFrame with the requested columns
    # 'nearest_hospital_name' comes from hospitals['name']
    # 'distance_m' comes from 'dist' (which is in metres if using EPSG:2056)
    
    # We must ensure distance_m is non-negative and finite
    dist = np.where(np.isfinite(dist), np.maximum(0, dist), np.nan)
    
    # Build the result DataFrame
    result_df = pd.DataFrame({
        'address_id': addresses['address_id'],
        'nearest_hospital_name': hospitals.iloc[hosp_indices]['name'].astype(str).values,
        'distance_m': dist
    })
    
    # Ensure non-empty strings for nearest_hospital_name (if possible)
    # If a name is empty or NaN, we should handle it. 
    # But typically it shouldn't be for a hospital list.
    # Replacing NaN with empty string as per request if needed, but "non-empty string" 
    # usually implies we should avoid empty strings. 
    # I'll just make sure they are strings.
    mask = result_df['nearest_hospital_name'].str.strip() == ''
    result_df.loc[mask, 'nearest_hospital_name'] = "Unknown"

    # Convert to GeoDataFrame
    # Use the original addresses geometry and CRS
    result_gdf = gpd.GeoDataFrame(result_df, geometry=addresses.geometry, crs=addresses.crs)

    # Write to GeoPackage
    result_gdf.to_file('nearest_hospital.gpkg', driver='GPKG')
    print("Success! Result written to nearest_hospital.gpkg")

if __name__ == "__main__":
    solve()

import geopandas as gpd
import pandas as pd
import numpy as np
from scipy.spatial import cKDTree

def solve():
    # Load data
    addresses = gpd.read_parquet('addresses.parquet')
    hospitals = gpd.read_parquet('hospitals.parquet')

    print(f"Addresses: {len(addresses)}")
    print(f"Hospitals: {len(hospitals)}")

    if addresses.crs is None or hospitals.crs is None:
        raise ValueError("CRS not found in input files")
    
    # Use a projected CRS for accurate distance calculation in metres
    # If it is geographic, project to a metric CRS.
    if addresses.crs.is_geographic:
        print("Projecting to EPSG:2056 for distance calculation...")
        target_crs = 'EPSG:2056'
    else:
        target_crs = addresses.crs
        print(f"Using existing CRS: {target_crs}")

    # Project both to the target metric CRS
    addresses_proj = addresses.to_crs(target_crs)
    hospitals_proj = hospitals.to_crs(target_crs)

    # Use cKDTree for efficient nearest neighbor search
    # Extract coordinates (using centroid in case they are not points)
    addr_coords = np.array([(p.x, p.y) for p in addresses_proj.geometry.centroid])
    hosp_coords = np.array([(p.x, p.y) for p in hospitals_proj.geometry.centroid])
    
    tree = cKDTree(hosp_coords)
    dist, hosp_indices = tree.query(addr_coords)

    # Ensure distance_m is non-negative and finite
    dist = np.where(np.isfinite(dist), np.maximum(0.0, dist), np.nan)
    
    # Prepare the result columns
    # Get the names from hospitals using the indices found
    hospital_names = hospitals.iloc[hosp_indices]['name'].values
    
    # Create a temporary dataframe to build the result
    result_data = {
        'address_id': addresses['address_id'].values,
        'nearest_hospital_name': hospital_names.astype(str),
        'distance_m': dist.astype(float)
    }
    
    # Handle empty strings if any
    # The requirement says: "nearest_hospital_name comes through as a non-empty string"
    # If name is empty or NaN, we replace it.
    for i in range(len(result_data['nearest_hospital_name'])):
        if not result_data['nearest_hospital_name'][i].strip():
            result_data['nearest_hospital_name'][i] = "Unknown"

    # Create the output GeoDataFrame
    # We keep original address geometry and original CRS
    result_gdf = gpd.GeoDataFrame(
        result_data, 
        geometry=addresses.geometry, 
        crs=addresses.crs
    )

    # Write to GeoPackage
    result_gdf.to_file('nearest_hospital.gpkg', driver='GPKG')
    print("Success! Result written to nearest_hospital.gpkg")

if __name__ == "__main__":
    solve()

import geopandas as gpd
import pandas as pd
import numpy as np
from scipy.spatial import cKDTree
import os

def solve():
    # Load data
    addresses = gpd.read_parquet('addresses.parquet')
    hospitals = gpd.read_parquet('hospitals.parquet')

    print(f"Addresses: {len(addresses)}")
    print(f"Hospitals: {len(hospitals)}")

    if addresses.crs is None or hospitals.crs is None:
        raise ValueError("CRS not found in input files")
    
    # Use a projected CRS for accurate distance calculation in metres
    if addresses.crs.is_geographic:
        print("Projecting to EPSG:2056 for distance calculation...")
        target_crs = 'EPSG:2056'
    else:
        target_crs = addresses.crs
        print(f"Using existing CRS: {target_crs}")

    # Project both to the target metric CRS
    addresses_proj = addresses.to_crs(target_crs)
    hospitals_proj = hospitals.to_crs(target_crs)

    # Use cKDTree for efficient nearest neighbor search
    addr_coords = np.array([(p.x, p.y) for p in addresses_proj.geometry.centroid])
    hosp_coords = np.array([(p.x, p.y) for p in hospitals_proj.geometry.centroid])
    
    tree = cKDTree(hosp_coords)
    dist, hosp_indices = tree.query(addr_coords)

    # Ensure distance_m is non-negative and finite
    dist = np.where(np.isfinite(dist), np.maximum(0.0, dist), np.nan)
    
    # Get names from hospitals
    hospital_names = hospitals.iloc[hosp_indices]['name'].values
    
    # Clean hospital names
    hospital_names = [str(name) if (name and str(name).strip()) else "Unknown" for name in hospital_names]

    # Create the output GeoDataFrame
    result_gdf = gpd.GeoDataFrame({
        'address_id': addresses['address_id'].values,
        'nearest_hospital_name': hospital_names,
        'distance_m': dist.astype(float),
        'geometry': addresses.geometry
    }, crs=addresses.crs)

    # Write to GeoPackage
    result_gdf.to_file('nearest_hospital.gpkg', driver='GPKG')
    print("Success! Result written to nearest_hospital.gpkg")

if __name__ == "__main__":
    solve()

import geopandas as gpd
import pandas as pd
import numpy as np
from scipy.spatial import cKDTree

def solve():
    # Load data
    addresses = gpd.read_parquet('addresses.parquet')
    hospitals = gpd.read_parquet('hospitals.parquet')

    print(f"Addresses: {len(addresses)}")
    print(f"Hospitals: {len(hospitals)}")

    if addresses.crs is None or hospitals.crs is None:
        raise ValueError("CRS not found in input files")
    
    # Use a projected CRS for accurate distance calculation in metres
    if addresses.crs.is_geographic:
        print("Projecting to EPSG:2056 for distance calculation...")
        target_crs = 'EPSG:2056'
    else:
        target_crs = addresses.crs
        print(f"Using existing CRS: {target_crs}")

    # Project both to the target metric CRS
    addresses_proj = addresses.to_crs(target_crs)
    hospitals_proj = hospitals.to_crs(target_crs)

    # Use cKDTree for efficient nearest neighbor search
    addr_coords = np.array([(p.x, p.y) for p in addresses_proj.geometry.centroid])
    hosp_coords = np.array([(p.x, p.y) for p in hospitals_proj.geometry.centroid])
    
    tree = cKDTree(hosp_coords)
    dist, hosp_indices = tree.query(addr_coords)

    # Ensure distance_m is non-negative and finite
    dist = np.where(np.isfinite(dist), np.maximum(0.0, dist), np.nan)
    
    # Get names from hospitals
    hospital_names = hospitals.iloc[hosp_indices]['name'].values
    
    # Clean hospital names
    hospital_names = [str(name) if (pd.notna(name) and str(name).strip()) else "Unknown" for name in hospital_names]

    # Create the output GeoDataFrame
    # Keep original address_id, geometry and CRS
    result_gdf = gpd.GeoDataFrame({
        'address_id': addresses['address_id'].values,
        'nearest_hospital_name': hospital_names,
        'distance_m': dist.astype(float),
        'geometry': addresses.geometry.values
    }, crs=addresses.crs)

    # Write to GeoPackage
    result_gdf.to_file('nearest_hospital.gpkg', driver='GPKG')
    print("Success! Result written to nearest_hospital.gpkg")

if __name__ == "__main__":
    solve()

import geopandas as gpd
import pandas as pd
import numpy as np

def solve():
    # Load data
    addresses = gpd.read_parquet('addresses.parquet')
    hospitals = gpd.read_parquet('hospitals.parquet')

    print(f"Addresses: {len(addresses)}")
    print(f"Hospitals: {len(hospitals)}")

    if addresses.crs is None or hospitals.crs is None:
        raise ValueError("CRS not found in input files")
    
    # Use a projected CRS for accurate distance calculation in metres
    if addresses.crs.is_geographic:
        print("Projecting to EPSG:2056 for distance calculation...")
        target_crs = 'EPSG:2056'
    else:
        target_crs = addresses.crs
        print(f"Using existing CRS: {target_crs}")

    # Project both to the target metric CRS
    addresses_proj = addresses.to_crs(target_crs)
    hospitals_proj = hospitals.to_crs(target_crs)

    # Use spatial index to find nearest neighbors efficiently
    # We'll use the sindex.query with method='nearest' if possible, 
    # but to be safe and compatible, we use a more manual approach with sindex.
    
    # Create a list to store results
    nearest_hospital_names = []
    distances = []
    
    # Get hospital geometries and names
    h_geom = hospitals_proj.geometry
    h_names = hospitals['name'].values

    # For each address, find the nearest hospital
    # We can use the sindex of the hospitals to speed up search
    h_sindex = hospitals_proj.sindex

    for idx, addr in addresses_proj.iterrows():
        # Use the spatial index to find the nearest hospital
        # sindex.nearest is available in newer geopandas (with pygeos/shapely 2.0)
        # It returns indices of the nearest geometries.
        try:
            # This is the most efficient way in modern geopandas
            # It returns a tuple of (indices_of_addresses, indices_of_hospitals)
            # Since we are looking at one address at a time, we provide it as a single-item collection
            # However, sindex.nearest is usually a batch operation. 
            # Let's try a more compatible approach for single items:
            
            # For small-ish number of hospitals, we can use a spatial query.
            # But let's assume we need speed.
            
            # We'll use the batch approach for all addresses at once for performance.
            # sindex.query(other_geometries, method='nearest') is quite efficient.
            pass
        except Exception:
            pass

    # Let's use the batch approach:
    # sindex.query returns the indices of neighbors. 
    # We want to find the one with the minimum distance for each address.
    
    # Using the batch nearest approach from geopandas:
    # This returns a 2D array of indices of hospitals near each address.
    # Actually, we can use the 'nearest' method on sindex.
    # Since we don't know the exact geopandas/shapely version, let's use a loop with sindex.query
    
    # BUT, let's check if we can use a simple distance-based approach with sindex.
    # If the datasets aren't massive, a loop with sindex is fine.
    
    # Efficient approach:
    # For each address, find the nearest hospital using the spatial index.
    
    print("Calculating nearest hospitals...")
    
    # To avoid the loop being too slow, we'll use the sindex for batching.
    # If sindex.query(..., method='nearest') is not available, we'll use a different approach.
    
    # Let's try to use the property that sindex.nearest returns the nearest indices for a set of geometries
    # if the version is recent enough.
    
    try:
        # Try the modern batch approach
        # This returns indices of hospitals that are nearest to the input geometries (addresses)
        # It's highly efficient.
        nearest_indices = addresses_proj.sindex.query(hospitals_proj.geometry, predicate='intersects') # This is for overlap. 
        # We need nearest.
        
        # Let's use a fallback for finding nearest neighbor if sindex.query doesn't support 'nearest'
        # Actually, the most reliable way for 'nearest' in modern geopandas is:
        # idxs = addresses_proj.sindex.nearest(hospitals_proj.geometry)
        # But 'nearest' in sindex is for the other way around? 
        # Let's use a loop with sindex.query and a large enough radius, or just iterate.
        
        # If I can't use scipy, I'll use a simple loop with sindex.query which is still much faster than no index.
        
        # Let's just do the loop. For typical EMS data (thousands of points), it's fine.
        for idx, addr in addresses_proj.iterrows():
            # Find index of nearest hospital
            # We query the spatial index for the nearest hospital to the current address
            # We can use the distance to the nearest object.
            
            # A trick to find the nearest using sindex:
            # Query with a very small radius or just find the nearest in the tree.
            # Since I can't use scipy, I'll use the fact that hospitals are points (usually)
            # and use the sindex.
            
            # For each address, we search for the nearest hospital.
            # We'll use the sindex.query with a "nearest" approach if available or 
            # simply use the sindex to get potential candidates.
            
            # Actually, the simplest way is to use the geometry's distance to all hospitals? 
            # No, that's O(N*M).
            
            # Let's try:
            # 1. Get the bounding box of the address (itself)
            # 2. Use sindex.query to find hospitals in a growing radius? No.
            
            # Let's just use the loop and for each address find the nearest hospital.
            # To make it faster, we'll use sindex.query to get a small set of neighbors.
            
            # If there are no scipy/pygeos, we'll just iterate. 
            # Let's see if it's too slow.
            
            pass
    except:
        pass

    # Let's reconsider. The most robust way without scipy:
    # Use a loop and for each address, find the closest hospital.
    # We can use the spatial index to find candidates.
    
    # Actually, let's try to use the sindex.query for all addresses against all hospitals to find candidates,
    # but that's for intersection.
    
    # Let's use this:
    # For each address, use sindex.query(addr.geometry) to get nearby hospitals, then find the min.
    # If no hospitals found in a small box, expand the box.
    
    # BUT, there's a simpler way:
    # For each address, we find the nearest hospital.
    
    # Let's try the most efficient way possible in pure geopandas/shapely:
    # Pre-calculate hospital coordinates.
    # Use a loop.
    
    h_coords = np.array([(p.x, p.y) for p in hospitals_proj.geometry])
    a_coords = np.array([(p.x, p.y) for p in addresses_proj.geometry])
    
    # Since I can't use scipy.spatial.cKDTree, I'll use the sindex to find candidates.
    # Or, if the number of addresses and hospitals is small (e.g. < 10,000 each), 
    # a vectorized distance calculation with numpy might be acceptable but memory intensive.
    # Let's try to process in chunks.

    print("Starting distance calculation...")
    
    # Result lists
    res_names = []
    res_dists = []
    res_ids = []

    # Process in chunks to manage memory
    chunk_size = 1000
    for i in range(0, len(addresses), chunk_size):
        chunk = addresses_proj.iloc[i : i + chunk_size]
        chunk_orig = addresses.iloc[i : i + chunk_size]
        
        # Current chunk coordinates
        c_coords = np.array([(p.x, p.y) for p in chunk.geometry])
        
        # For each address in chunk, find the nearest hospital.
        # We'll use the sindex to find the nearest one.
        for j, (idx_a, addr_a) in enumerate(chunk.iterrows()):
            # Use sindex to find the nearest hospital.
            # In modern geopandas, sindex.nearest is the way.
            # Let's check if it works.
            try:
                # This returns the index of the nearest hospital
                nearest_h_idx = list(h_sindex.nearest(addr_a.geometry, return_distance=False))[0]
                # In some versions, sindex.nearest(geom) returns a tuple or array.
                # Let's try to be robust.
                if isinstance(nearest_h_idx, (list, np.ndarray)):
                    nearest_h_idx = nearest_h_idx[0]
            except:
                # Fallback: if sindex.nearest is not available, use distance to all hospitals
                # This is slow but will work if the number of hospitals is not huge.
                # We'll use the chunk's coordinates to do a vectorized distance calculation.
                # For a single address:
                dists_to_all = chunk_orig.geometry.iloc[j].distance(h_proj_geoms) # This is still O(M)
                # Wait, let's just use a simple distance calculation for each address.
                
                # Let's try a different way for the fallback:
                # For each address, find the hospital with minimum distance.
                # We'll use a loop. 
                
                # To speed up, we'll use the sindex to get a candidate list.
                # For simplicity, let's just use the distance to all hospitals.
                # If it's too slow, we'll optimize.
                pass

            # Actually, let's use the most compatible way:
            # For each address: 
            #   1. Get the nearest hospital index using sindex.query
            #   2. If not found, use distance to all.
            
            # Let's refine the loop.
            pass

    # Actually, let's just use the pure Python/Geopandas approach. 
    # It's usually fast enough for several thousand points.
    
    print("Looping through addresses...")
    for idx, addr in addresses_proj.iterrows():
        # This is the simplest way to find the nearest geometry in a GeoSeries
        # using only geopandas/shapely.
        # We find the index of the minimum distance.
        
        # Note: addr.distance(hospitals_proj.geometry) is a Series of distances.
        dists = addr.distance(hospitals_proj.geometry)
        min_idx = dists.idxmin()
        min_dist = dists.min()
        
        res_ids.append(addresses.loc[idx, 'address_id'])
        res_names.append(hospitals.loc[min_idx, 'name'])
        res_dists.append(min_dist)
        
    # Clean names and build result
    res_names = [str(n) if (pd.notna(n) and str(n).strip()) else "Unknown" for n in res_names]
    
    result_gdf = gpd.GeoDataFrame({
        'address_id': res_ids,
        'nearest_hospital_name': res_names,
        'distance_m': res_dists,
        'geometry': addresses.geometry
    }, crs=addresses.crs)

    result_gdf.to_file('nearest_hospital.gpkg', driver='GPKG')
    print("Success! Result written to nearest_hospital.gpkg")

if __name__ == "__main__":
    solve()

import geopandas as gpd
import pandas as pd
import numpy as np

def solve():
    # Load data
    addresses = gpd.read_parquet('addresses.parquet')
    hospitals = gpd.read_parquet('hospitals.parquet')

    print(f"Addresses: {len(addresses)}")
    print(f"Hospitals: {len(hospitals)}")

    if addresses.crs is None or hospitals.crs is None:
        raise ValueError("CRS not found in input files")
    
    # Use a projected CRS for accurate distance calculation in metres
    if addresses.crs.is_geographic:
        print("Projecting to EPSG:2056 for distance calculation...")
        target_crs = 'EPSG:2056'
    else:
        target_crs = addresses.crs
        print(f"Using existing CRS: {target_crs}")

    # Project both to the target metric CRS
    addresses_proj = addresses.to_crs(target_crs)
    hospitals_proj = hospitals.to_crs(target_crs)

    # Result lists
    res_ids = []
    res_names = []
    res_dists = []

    print("Looping through addresses...")
    # We'll iterate through the addresses and for each address, find the nearest hospital.
    # Since there are only 120 addresses and 37 hospitals, a nested loop is fine.
    # O(N*M) is 120 * 37 = 4440 operations, which is negligible.

    for idx, addr in addresses_proj.iterrows():
        # addr is a GeoSeries (a single row from addresses_proj)
        # We want its geometry
        geom = addr.geometry
        
        # Calculate distance to all hospital geometries
        # hospitals_proj.geometry is a GeoSeries
        # In geopandas, you can call distance on the geometry object against a GeoSeries
        dists = geom.distance(hospitals_proj.geometry)
        
        # Find the index of the minimum distance
        min_dist = dists.min()
        min_idx_in_hospitals_proj = dists.idxmin()
        
        # Get the hospital name using the index from the original hospitals dataframe
        # Since hospitals_proj is a projection of hospitals, it shares the same index.
        hospital_name = hospitals.loc[min_idx_in_hospitals_proj, 'name']
        
        # Store results
        res_ids.append(addresses.loc[idx, 'address_id'])
        res_names.append(str(hospital_name) if (pd.notna(hospital_name) and str(hospital_name).strip()) else "Unknown")
        res_dists.append(float(min_dist))

    # Create the output GeoDataFrame
    # Use original addresses geometry and CRS
    result_gdf = gpd.GeoDataFrame({
        'address_id': res_ids,
        'nearest_hospital_name': res_names,
        'distance_m': res_dists,
        'geometry': addresses.geometry
    }, crs=addresses.crs)

    # Write to GeoPackage
    result_gdf.to_file('nearest_hospital.gpkg', driver='GPKG')
    print("Success! Result written to nearest_hospital.gpkg")

if __name__ == "__main__":
    solve()