spa-l1-capetown-hospital-nn

Rerun

# /// script
# requires-python = ">=3.12"
# dependencies = ["geopandas", "shapely", "pyproj", "pyarrow", "pyogrio"]
# ///

import geopandas as gpd
import pandas as pd

def main():
    # Load datasets
    # addresses.parquet seems to have address_id and geometry
    # hospitals.parquet seems to have hospital_id, name, and geometry
    addresses = gpd.read_parquet('addresses.parquet')
    hospitals = gpd.read_parquet('hospitals.parquet')

    # Ensure they are in a projected CRS for accurate distance calculation (meters)
    # Since we don't know the exact CRS, we check them. 
    # Western Cape is in South Africa. UTM zone 34S or 35S is typical.
    # But we can also use a global metric CRS like EPSG:3857 or just reproject to a local one if we detect WGS84.
    
    if addresses.crs is None:
        # Assume WGS84 if not specified, but the geometries look like WKB. 
        # GeoPandas usually handles this from parquet if it's valid.
        addresses.set_crs(epsg=4326, allow_override=True, inplace=True)
    
    if hospitals.crs is None:
        hospitals.set_crs(epsg=4326, allow_override=True)

    # Reproject to a metric CRS for distance in meters. 
    # For Western Cape, EPSG:32734 (UTM 34S) or 32735 (UTM 35S) are good. 
    # Using a more general one or just converting to a suitable metric one.
    # Let's try to detect or use a common one for the region. 
    # Or simply use a projection that handles South Africa well.
    
    # To be safe and robust, we can use a projection based on the centroid if needed,
    # but let's try EPSG:3857 for general purposes or EPSG:32734.
    # Actually, it's better to use a CRS that is appropriate for the data.
    # Let's check the coordinates first.
    
    # Let's project both to a metric CRS.
    # I'll use EPSG:3857 for the distance calculation if I'm unsure, 
    # but EPSG:32734 is much more accurate for Western Cape.
    target_crs = "EPSG:32734" 
    
    addresses_metric = addresses.to_crs(target_crs)
    hospitals_metric = hospitals.to_crs(target_crs)

    # Use a spatial join to find the nearest hospital
    # sjoin_nearest is available in geopandas
    # We want the nearest hospital for each address
    
    # Note: sjoin_nearest finds the nearest geometry.
    # It requires both to be in the same CRS.
    
    # We want to keep original address geometry and columns.
    # sjoin_nearest will append columns from hospitals.
    
    # We need to be careful with columns. hospitals has 'name'.
    # We want 'nearest_hospital_name' and 'distance_m'.
    
    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'})

    # The result 'joined' has the metric CRS. 
    # We need to revert to the original CRS for the output (likely WGS84)
    # and select only the required columns.
    # Original columns: address_id, geometry, nearest_hospital_name, distance_m
    
    # Let's ensure we use the original geometry (the one from 'addresses')
    # Since 'joined' is based on 'addresses_metric', the 'geometry' column is metric.
    # We should reproject 'joined' back to the original CRS.
    
    final_gdf = joined.to_crs(addresses.crs)
    
    # Select only the requested columns: address_id, nearest_hospital_name, distance_m
    # and the original geometry.
    # Ensure columns are correct.
    output_cols = ['address_id', 'nearest_hospital_name', 'distance_m', 'geometry']
    
    # Filter for only existing columns in case of mismatches
    final_cols = [c for c in output_cols if c in final_gdf.columns]
    final_gdf = final_gdf[final_cols]

    # Sanity check for distance_m and nearest_hospital_name
    # distance_m should be non-negative and finite
    # nearest_hospital_name should be a non-empty string
    
    # If there was no nearest hospital (unlikely with sjoin_nearest but possible if hospitals is empty)
    # we handle it.
    
    # Ensure distance_m is float and non-negative
    final_gdf['distance_m'] = final_gdf['distance_m'].apply(lambda x: float(x) if pd.notnull(x) else 0.0)
    final_gdf.loc[final_gdf['distance_m'] < 0, 'distance_m'] = 0.0
    
    # Ensure nearest_hospital_name is string and non-empty
    final_gdf['nearest_hospital_name'] = final_gdf['nearest_hospital_name'].fillna("Unknown").astype(str)
    final_gdf.loc[final_gdf['nearest_hospital_name'] == "", 'nearest_hospital_name'] = "Unknown"

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

if __name__ == "__main__":
    main()