spa-l1-capetown-hospital-nn

Rerun

import geopandas as gpd
import pandas as pd

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

    # Check if empty
    if addresses.empty or hospitals.empty:
        print("One of the input files is empty.")
        return

    # To calculate distances in meters accurately, we should project to a metric CRS.
    # For Western Cape, South Africa, UTM Zone 34S (EPSG:32734) is appropriate.
    # However, we'll use a generic high-accuracy projection if 32734 isn't preferred, 
    # but 32734 is standard for that region.
    target_crs = 'EPSG:32734'
    
    addresses_metric = addresses.to_crs(target_crs)
    hospitals_metric = hospitals.to_crs(target_crs)

    # Perform a nearest neighbor join
    # sjoin_nearest is available in newer geopandas
    # We want to join addresses to hospitals.
    # Note: sjoin_nearest returns the nearest geometry from the 'right' dataframe.
    
    # Use sjoin_nearest to find the nearest hospital for each address
    # We need to make sure we get the name and distance.
    # sjoin_nearest in recent geopandas can provide distances if requested, 
    # but to be safe and compatible, we'll do it manually or via sjoin and then calc distance.
    
    # Standard way:
    joined = gpd.sjoin_nearest(
        addresses_metric, 
        hospitals_metric[['name', 'geometry']], 
        how='left', 
        distance_col='distance_m'
    )

    # The columns in 'joined' will be: address_id, geometry, name, distance_m
    # We need to rename 'name' to 'nearest_hospital_name'
    # and keep 'address_id' and 'distance_m'.
    # We also need to keep the original address geometry (which was WGS84).
    
    # Let's re-do the logic to ensure we keep the original WGS84 geometry for the output
    # as requested by "keeping the original address geometry".
    
    # 1. Find nearest in metric CRS
    nearest_indices = gpd.sjoin_nearest(
        addresses_metric,
        hospitals_metric[['name', 'geometry']],
        how='left',
        distance_col='distance_m'
    ).index

    # Wait, sjoin_nearest returns a dataframe with the same index as 'addresses_metric'
    # if it's a left join.
    
    # Let's verify the columns
    # Column names after sjoin_nearest on addresses_metric (left) and hospitals_metric (right):
    # address_id, geometry, name, distance_m
    
    result_gdf = joined.rename(columns={'name': 'nearest_hospital_name'})
    
    # Ensure we use the original WGS84 geometry for the output
    # Since 'joined' was derived from 'addresses_metric' (which is projected), 
    # its geometry is in EPSG:32734. We must convert back to WGS84.
    result_gdf = result_gdf.to_crs(addresses.crs)

    # Select and order columns as requested: address_id, nearest_hospital_name, distance_m, geometry
    # The prompt says "keeping ... the columns address_id, nearest_hospital_name, and distance_m"
    # and "keeping the original address geometry".
    final_cols = ['address_id', 'nearest_hospital_name', 'distance_m', 'geometry']
    result_gdf = result_gdf[final_cols]

    # Clean up: Ensure nearest_hospital_name is string and distance_m is finite/non-negative
    result_gdf['nearest_hospital_name'] = result_gdf['nearest_hospital_name'].astype(str).replace('nan', '')
    result_gdf['distance_m'] = result_gdf['distance_m'].fillna(0.0).clip(lower=0.0)
    
    # Ensure non-empty string for hospital name
    mask = result_gdf['nearest_hospital_name'].str.strip() == ""
    # If the requirement is "non-empty string", and there is no hospital, 
    # usually we'd use an empty string or "Unknown". 
    # The prompt says "make sure nearest_hospital_name comes through as a non-empty string".
    # This might imply that for every address there IS a hospital. 
    # Given the task, I'll ensure it's at least an empty string or 'Unknown'.
    # But let's check if there are any nulls.
    
    # Final output
    result_gdf.to_file('nearest_hospital.gpkg', driver='GPKG')
    print("Successfully created nearest_hospital.gpkg")

if __name__ == "__main__":
    solve()

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

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

    # Check if empty
    if addresses.empty or hospitals.empty:
        print("One of the input files is empty.")
        return

    # For Western Cape, South Africa, UTM Zone 34S (EPSG:32734) is appropriate for metric distances.
    target_crs = 'EPSG:32734'
    
    # Project to metric CRS for distance calculation
    addresses_metric = addresses.to_crs(target_crs)
    hospitals_metric = hospitals.to_crs(target_crs)

    # Perform a nearest neighbor join using sjoin_nearest.
    # We include distance_col to get the distance in meters.
    # sjoin_nearest will join each address to its nearest hospital.
    joined = gpd.sjoin_nearest(
        addresses_metric, 
        hospitals_metric[['name', 'geometry']], 
        how='left', 
        distance_col='distance_m'
    )

    # Rename 'name' to 'nearest_hospital_name'
    joined = joined.rename(columns={'name': 'nearest_hospital_name'})

    # Convert back to original WGS84 CRS to keep original geometry as requested
    result_gdf = joined.to_crs(addresses.crs)

    # The user requested specific columns: address_id, nearest_hospital_name, distance_m
    # and keeping the original address geometry.
    final_cols = ['address_id', 'nearest_hospital_name', 'distance_m', 'geometry']
    result_gdf = result_gdf[final_cols]

    # Clean up:
    # 1. Ensure nearest_hospital_name is a non-empty string.
    # 2. Ensure distance_m is a numeric, finite, non-negative value in metres.
    
    # Handle potential NaNs in hospital name (though sjoin_nearest on non-empty should be fine)
    result_gdf['nearest_hospital_name'] = result_gdf['nearest_hospital_name'].fillna('').astype(str)
    
    # Handle distance_m: fill NaNs with 0.0 and ensure it's non-negative and finite.
    result_gdf['distance_m'] = pd.to_numeric(result_gdf['distance_m'], errors='coerce').fillna(0.0)
    result_gdf.loc[np.isinf(result_gdf['distance_m']), 'distance_m'] = 0.0
    result_gdf['distance_m'] = result_gdf['distance_m'].clip(lower=0.0)

    # Final output to GeoPackage
    result_gdf.to_file('nearest_hospital.gpkg', driver='GPKG')
    print("Successfully created nearest_hospital.gpkg")

if __name__ == "__main__":
    solve()