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')

    # Ensure we are working with a projected CRS for accurate meter distances
    # Western Cape is roughly in UTM zone 34S (EPSG:32734) or 35S (EPSG:32735)
    # However, a more robust way is to check the input CRS and project to a local UTM or use a global metric CRS like World Equidistant Cylindrical or similar, 
    # but for local distance calculations, UTM is best.
    # Let's check the input CRS first.
    
    if addresses.crs is None:
        # If no CRS, assume WGS84 (standard for many datasets)
        addresses.set_crs(epsg=4326, inplace=True)
        hospitals.set_crs(epsg=4326, inplace=True)
    
    # Project to a metric CRS. EPSG:3857 is often used but can be distorted.
    # Since it's Western Cape, let's try to find a suitable UTM.
    # For simplicity and accuracy in straight-line distance in meters, 
    # we will project to a metric system.
    # Let's use EPSG:3857 for a general approach or a local UTM if we want to be precise.
    # But actually, geopandas's distance method on a projected CRS is what we want.
    
    # We will project both to a metric CRS. 
    # A safe bet for distance calculation if we don't want to guess the exact UTM is to 
    # project to a metric CRS like EPSG:3857 or use a local UTM.
    # Let's detect the UTM zone for Western Cape (usually 34S or 35S).
    # Or just use a common one. 32734 is common for South Africa.
    
    target_crs = 'EPSG:3857' # Web Mercator is metric, though distorted at high latitudes. 
    # For Western Cape (lat ~ -34), it's okay-ish, but UTM 34S (32734) is much better.
    # Let's try 32734.
    target_crs = 'EPSG:32734'

    addresses_m = addresses.to_crs(target_crs)
    hospitals_m = hospitals.to_crs(target_crs)

    # Use a spatial index for efficiency
    sindex = hospitals_m.sindex

    nearest_hospital_names = []
    distances = []

    for idx, addr in addresses_m.iterrows():
        # Find the nearest hospital
        # We can use the spatial index to find the nearest neighbor
        # nearest_neighbor = sindex.nearest(addr.geometry, return_distance=True) - this is for newer scipy/geopandas
        
        # Traditional way with spatial index:
        # Find the nearest point in hospitals_m to addr.geometry
        
        # Get the nearest geometry from hospitals
        # Using distance directly on the whole hospitals_m is slow for large datasets,
        # but since we want a robust solution:
        
        # Find index of nearest hospital
        # We'll use the distance method on the hospitals geometry set
        dist_array = hospitals_m.distance(addr.geometry)
        min_dist_idx = dist_array.idxmin()
        
        min_dist = dist_array.min()
        h_name = hospitals.loc[min_dist_idx, 'name']
        
        nearest_hospital_names.append(h_name)
        distances.append(min_dist)

    # Add results to the original GeoDataFrame (to keep original geometry/CRS)
    # But the prompt says "keeping the original address geometry" and "one feature per input address"
    # so we merge back.
    
    addresses['nearest_hospital_name'] = nearest_hospital_names
    addresses['distance_m'] = distances

    # Ensure requirements:
    # nearest_hospital_name: non-empty string
    # distance_m: numeric, finite, non-negative
    
    addresses['nearest_hospital_name'] = addresses['nearest_hospital_name'].fillna('').astype(str)
    addresses['distance_m'] = pd.to_numeric(addresses['distance_m'], errors='coerce').fillna(0.0)
    addresses = addresses[addresses['distance_m'] >= 0]

    # Keep only requested columns: address_id, nearest_hospital_name, distance_m and geometry
    # Note: addresses already has geometry.
    output = addresses[['address_id', 'nearest_hospital_name', 'distance_m', 'geometry']]

    # Write to GeoPackage
    output.to_file('nearest_hospital.gpkg', driver='GPKG')
    print("Successfully wrote nearest_hospital.gpkg")

if __name__ == "__main__":
    solve()

import geopandas as gpd
import pandas as pd

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

    # Ensure CRS is set if not present. Standard for GPS is EPSG:4326.
    if addresses.crs is None:
        addresses.set_crs(epsg=4326, inplace=True)
    if hospitals.crs is None:
        hospitals.set_crs(epsg=4326, inplace=True)

    # Project to a metric CRS for accurate distance calculation in meters.
    # UTM Zone 34S (EPSG:32734) is appropriate for the Western Cape region.
    target_crs = 'EPSG:32734'
    addresses_m = addresses.to_crs(target_crs)
    hospitals_m = hospitals.to_crs(target_crs)

    # Initialize lists to store results
    nearest_hospital_names = []
    distances = []

    # Use the spatial index of hospitals for efficient proximity search.
    # For each address, find the nearest hospital.
    # We use the distance method on the hospitals_m GeoDataFrame.
    # While a full distance calc is O(N*M), it's robust for medium datasets.
    # For much larger datasets, we would use a spatial index (sindex.query).
    
    # To make it faster for potentially larger inputs, let's use the sindex.
    # However, standard distance() on a GeoSeries/GeoDataFrame is well-optimized.
    # Given the likely size of the problem, we'll use an efficient approach.
    
    # Prepare for fast nearest neighbor search
    # We'll iterate through addresses and find the nearest hospital using sindex
    sindex = hospitals_m.sindex

    for idx, addr in addresses_m.iterrows():
        # Query the spatial index for the nearest point
        # sindex.nearest is available in newer versions of geopandas/pygeos/shapely
        # If not, we use the distance method on the whole set or a subset.
        
        # A simple but efficient enough method for this context:
        dist_to_hospitals = hospitals_m.distance(addr.geometry)
        min_dist_idx = dist_to_hospitals.idxmin()
        min_dist = dist_to_hospitals.min()
        
        # The index of hospitals_m might not be the same as the index of hospitals 
        # if it was filtered, but here it's just a projection.
        # We want the 'name' from the original hospitals df using the same index.
        h_name = hospitals.loc[min_dist_idx, 'name']
        
        nearest_hospital_names.append(h_name)
        distances.append(min_dist)

    # Add results to the original GeoDataFrame (to preserve original CRS and geometry)
    addresses['nearest_hospital_name'] = nearest_hospital_names
    addresses['distance_m'] = distances

    # Clean up: ensure non-empty strings and finite, non-negative numeric values
    addresses['nearest_hospital_name'] = addresses['nearest_hospital_name'].fillna('').astype(str)
    addresses['distance_m'] = pd.to_numeric(addresses['distance_m'], errors='coerce').fillna(0.0)
    
    # Filter out any potentially invalid distances (though distance is >= 0 by definition)
    addresses = addresses[addresses['distance_m'] >= 0]

    # Select only the requested columns
    # Note: address_id, nearest_hospital_name, distance_m, and geometry
    output = addresses[['address_id', 'nearest_hospital_name', 'distance_m', 'geometry']]

    # Write to GeoPackage
    output.to_file('nearest_hospital.gpkg', driver='GPKG')
    print("Successfully wrote nearest_hospital.gpkg")

if __name__ == "__main__":
    solve()