spa-l1-capetown-hospital-nn
Map (reference ← swipe → agent)
0:00
I am putting together EMS coverage planning for the Western Cape and need each residential pickup address tagged with its nearest hospital so the ops-research team can join the result against shift schedules. For every address in `addresses.parquet`, find the nearest hospital from `hospitals.parquet` by straight-line distance and record that hospital's name plus the distance in metres. Keep the original `address_id` from the input so we can join back to the address master, and make sure `nearest_hospital_name` comes through as a non-empty string and `distance_m` as a numeric, finite, non-negative value in metres.
Please write the result to `nearest_hospital.gpkg` in GeoPackage format with one feature per input address, keeping the original address geometry and the columns `address_id`, `nearest_hospital_name`, and `distance_m`.
0:00
System info
{
"adapter": "openrouter",
"model": "deepseek/deepseek-v4-flash",
"base_url": "https://openrouter.ai/api/v1"
}
0:04
Let me start by examining the input files to understand their structure.
0:04
Bash
cd /work && python3 -c "
cd /work && python3 -c "
import pandas as pd
addr = pd.read_parquet('addresses.parquet')
hosp = pd.read_parquet('hospitals.parquet')
print('=== ADDRESSES ===')
print('Columns:', list(addr.columns))
print('Count:', len(addr))
print(addr.head(3))
print()
print('=== HOSPITALS ===')
print('Columns:', list(hosp.columns))
print('Count:', len(hosp))
print(hosp.head(3))
"
=== ADDRESSES === Columns: ['address_id', 'geometry'] Count: 120 address_id geometry 0 A0001 b'\x01\x01\x00\x00\x00\xc2\x99$,\xcf\x882@\xcd... 1 A0002 b'\x01\x01\x00\x00\x00\xfcg\xbf\x12\xc2\x802@t... 2 A0003 b'\x01\x01\x00\x00\x00\xb8\xa8\x80\x7fG\x882@r... === HOSPITALS === Columns: ['hospital_id', 'name', 'geometry'] Count: 37 hospital_id ... geometry 0 H001 ... b"\x01\x01\x00\x00\x00\xb7''^\x8ey2@\x12\xf4z\... 1 H002 ... b'\x01\x01\x00\x00\x00\xb5\x01THdw2@\x15\xe1`N... 2 H003 ... b'\x01\x01\x00\x00\x00\x0b\xe8K\x90\x8dk2@\xd6... [3 rows x 3 columns]
0:10
The geometry columns are stored as binary WKB/WKB bytes. Let me inspect more carefully to understand the data.
0:10
Bash
cd /work && python3 -c "
cd /work && python3 -c "
import pandas as pd
import pyogrio
addr = pd.read_parquet('addresses.parquet')
hosp = pd.read_parquet('hospitals.parquet')
print('addresses dtypes:', addr.dtypes)
print('hospitals dtypes:', hosp.dtypes)
print()
print('addresses geometry type:', type(addr['geometry'].iloc[0]))
print('hospitals geometry type:', type(hosp['geometry'].iloc[0]))
print()
# Try loading with geopandas
import geopandas as gpd
addr_gdf = gpd.read_parquet('addresses.parquet')
hosp_gdf = gpd.read_parquet('hospitals.parquet')
print('=== ADDR GDF ===')
print(addr_gdf.head(2))
print('crs:', addr_gdf.crs)
print()
print('=== HOSP GDF ===')
print(hosp_gdf.head(2))
print('crs:', hosp_gdf.crs)
"
addresses dtypes: address_id str
geometry object
dtype: object
hospitals dtypes: hospital_id str
name str
geometry object
dtype: object
addresses geometry type: <class 'bytes'>
hospitals geometry type: <class 'bytes'>
=== ADDR GDF ===
address_id geometry
0 A0001 POINT (18.53441 -33.94184)
1 A0002 POINT (18.50296 -33.93623)
crs: {"$schema": "https://proj.org/schemas/v0.7/projjson.schema.json", "type": "GeographicCRS", "name": "WGS 84", "datum_ensemble": {"name": "World Geodetic System 1984 ensemble", "members": [{"name": "World Geodetic System 1984 (Transit)"}, {"name": "World Geodetic System 1984 (G730)"}, {"name": "World Geodetic System 1984 (G873)"}, {"name": "World Geodetic System 1984 (G1150)"}, {"name": "World Geodetic System 1984 (G1674)"}, {"name": "World Geodetic System 1984 (G1762)"}, {"name": "World Geodetic System 1984 (G2139)"}, {"name": "World Geodetic System 1984 (G2296)"}], "ellipsoid": {"name": "WGS 84", "semi_major_axis": 6378137, "inverse_flattening": 298.257223563}, "accuracy": "2.0", "id": {"authority": "EPSG", "code": 6326}}, "coordinate_system": {"subtype": "ellipsoidal", "axis": [{"name": "Geodetic latitude", "abbreviation": "Lat", "direction": "north", "unit": "degree"}, {"name": "Geodetic longitude", "abbreviation": "Lon", "direction": "east", "unit": "degree"}]}, "scope": "Horizontal component of 3D system.", "area": "World.", "bbox": {"south_latitude": -90, "west_longitude": -180, "north_latitude": 90, "east_longitude": 180}, "id": {"authority": "EPSG", "code": 4326}}
=== HOSP GDF ===
hospital_id name geometry
0 H001 Mowbray Maternity Hospital POINT (18.47483 -33.94921)
1 H002 Life Claremont Hospital POINT (18.46637 -33.98695)
crs: {"$schema": "https://proj.org/schemas/v0.7/projjson.schema.json", "type": "GeographicCRS", "name": "WGS 84", "datum_ensemble": {"name": "World Geodetic System 1984 ensemble", "members": [{"name": "World Geodetic System 1984 (Transit)"}, {"name": "World Geodetic System 1984 (G730)"}, {"name": "World Geodetic System 1984 (G873)"}, {"name": "World Geodetic System 1984 (G1150)"}, {"name": "World Geodetic System 1984 (G1674)"}, {"name": "World Geodetic System 1984 (G1762)"}, {"name": "World Geodetic System 1984 (G2139)"}, {"name": "World Geodetic System 1984 (G2296)"}], "ellipsoid": {"name": "WGS 84", "semi_major_axis": 6378137, "inverse_flattening": 298.257223563}, "accuracy": "2.0", "id": {"authority": "EPSG", "code": 6326}}, "coordinate_system": {"subtype": "ellipsoidal", "axis": [{"name": "Geodetic latitude", "abbreviation": "Lat", "direction": "north", "unit": "degree"}, {"name": "Geodetic longitude", "abbreviation": "Lon", "direction": "east", "unit": "degree"}]}, "scope": "Horizontal component of 3D system.", "area": "World.", "bbox": {"south_latitude": -90, "west_longitude": -180, "north_latitude": 90, "east_longitude": 180}, "id": {"authority": "EPSG", "code": 4326}}
0:22
Now I have a clear picture. Both files are in WGS84 (EPSG:4326). I need to: 1. Read addresses and hospitals as GeoDataFrames 2. For each address, find the nearest hospital by straight-line distance (Euclidean distance in metres - but since coordinates are in degrees, I need to compute great-circle distance or project to a suitable local CRS for distance calculations). Since the data is in the Western Cape, a good local projected CRS would be something like EPSG:32734 (UTM zone 34S) which covers the Western Cape. Or I could use the haversine formula for great-circle distance. Let me use pyproj to project to a suitable CRS and then compute Euclidean distances. Actually, the simplest approach: use `geopandas.sjoin_nearest` with a projected CRS. Let me reproject to EPSG:32734 (UTM 34S) which covers the Western Cape. Let me write the script.
0:22
Write
/work/solve.py
content · 47 lines · py
import geopandas as gpd
import pandas as pd
import numpy as np
# Read input files
addresses = gpd.read_parquet('/work/addresses.parquet')
hospitals = gpd.read_parquet('/work/hospitals.parquet')
# Both are in WGS84 (EPSG:4326) based on inspection.
# Western Cape is in UTM zone 34S (EPSG:32734). We'll project for accurate
# straight-line distance calculations in metres.
utm_crs = "EPSG:32734"
addresses_utm = addresses.to_crs(utm_crs)
hospitals_utm = hospitals.to_crs(utm_crs)
# For each address, find the nearest hospital using spatial join
# sjoin_nearest returns the nearest neighbour
nearest = gpd.sjoin_nearest(
addresses_utm,
hospitals_utm[['hospital_id', 'name', 'geometry']],
how='left',
distance_col='distance_m'
)
# distance_m is already in metres since both layers are in a projected CRS
# Build result: keep original address geometry (in WGS84 to match input),
# plus address_id, nearest_hospital_name, distance_m
result = addresses[['address_id', 'geometry']].copy()
result['nearest_hospital_name'] = nearest['name'].values
result['distance_m'] = nearest['distance_m'].values
# Validate constraints
assert result['address_id'].notna().all(), "address_id must not be null"
assert result['nearest_hospital_name'].notna().all(), "nearest_hospital_name must not be null"
assert (result['nearest_hospital_name'].str.len() > 0).all(), "nearest_hospital_name must be non-empty"
assert result['distance_m'].notna().all(), "distance_m must not be null"
assert np.isfinite(result['distance_m']).all(), "distance_m must be finite"
assert (result['distance_m'] >= 0).all(), "distance_m must be non-negative"
# Write to GeoPackage
result.to_file('/work/nearest_hospital.gpkg', layer='nearest_hospital', driver='GPKG')
print(f"Wrote {len(result)} features to /work/nearest_hospital.gpkg")
print(result[['address_id', 'nearest_hospital_name', 'distance_m']].head(10))
print(f"\nDistance range: {result['distance_m'].min():.2f} m – {result['distance_m'].max():.2f} m")Wrote /work/solve.py (1991 chars).
0:25
Let me run the script.
0:25
Bash
cd /work && python solve.py
Wrote 120 features to /work/nearest_hospital.gpkg address_id nearest_hospital_name distance_m 0 A0001 Melomed Gatesville Private Hospital 3192.958737 1 A0002 Life Vincent Pallotti Hospital 1468.601852 2 A0003 Melomed Gatesville Private Hospital 1271.941611 3 A0004 Maitland Cottage Children's Orthopaedic Hospital 984.047189 4 A0005 Kingsbury Hospital 3021.737824 5 A0006 Intercare Physical Rehabilitation Hospital - N... 1777.422286 6 A0007 Melomed Gatesville Private Hospital 1268.558960 7 A0008 Melomed Gatesville Private Hospital 1264.240167 8 A0009 Melomed Gatesville Private Hospital 2020.478965 9 A0010 Red Cross Children's Hospital 1811.017934 Distance range: 47.02 m – 5956.74 m
0:30
The script ran successfully. Let me verify the output by reading back the GeoPackage to make sure everything looks good.
0:30
Bash
cd /work && python3 -c "
cd /work && python3 -c "
import geopandas as gpd
result = gpd.read_file('/work/nearest_hospital.gpkg')
print('Columns:', list(result.columns))
print('Count:', len(result))
print('CRS:', result.crs)
print()
print(result.head())
print()
print('dtypes:')
print(result.dtypes)
print()
print('Null check:')
print('address_id nulls:', result['address_id'].isna().sum())
print('nearest_hospital_name nulls:', result['nearest_hospital_name'].isna().sum())
print('distance_m nulls:', result['distance_m'].isna().sum())
print()
print('distance_m range:', result['distance_m'].min(), '-', result['distance_m'].max())
print('All finite:', result['distance_m'].apply(lambda x: __import__('numpy').isfinite(x)).all())
print('All non-neg:', (result['distance_m'] >= 0).all())
print('name non-empty:', (result['nearest_hospital_name'].str.len() > 0).all())
"
Columns: ['address_id', 'nearest_hospital_name', 'distance_m', 'geometry'] Count: 120 CRS: EPSG:4326 address_id ... geometry 0 A0001 ... POINT (18.53441 -33.94184) 1 A0002 ... POINT (18.50296 -33.93623) 2 A0003 ... POINT (18.53234 -33.95918) 3 A0004 ... POINT (18.44857 -33.98103) 4 A0005 ... POINT (18.49928 -33.99649) [5 rows x 4 columns] dtypes: address_id str nearest_hospital_name str distance_m float64 geometry geometry dtype: object Null check: address_id nulls: 0 nearest_hospital_name nulls: 0 distance_m nulls: 0 distance_m range: 47.01874251111449 - 5956.744689131489 All finite: True All non-neg: True name non-empty: True
0:40
All checks pass. Let me summarize the results.
0:40