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:09
Let me start by examining the input files to understand their structure.
0:09
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('Shape:', addr.shape)
print(addr.head(3))
print('dtypes:', addr.dtypes)
print()
print('=== HOSPITALS ===')
print('Columns:', list(hosp.columns))
print('Shape:', hosp.shape)
print(hosp.head(3))
print('dtypes:', hosp.dtypes)
"
=== ADDRESSES === Columns: ['address_id', 'geometry'] Shape: (120, 2) 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... dtypes: address_id str geometry object dtype: object === HOSPITALS === Columns: ['hospital_id', 'name', 'geometry'] Shape: (37, 3) 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] dtypes: hospital_id str name str geometry object dtype: object
0:15
The geometry columns are stored as WKB bytes (well-known binary). I need to convert them to proper geometry objects using shapely. Let me load them with geopandas instead.
0:15
Bash
cd /work && python3 -c "
cd /work && python3 -c "
import geopandas as gpd
addr = gpd.read_parquet('addresses.parquet')
hosp = gpd.read_parquet('hospitals.parquet')
print('=== ADDRESSES ===')
print('CRS:', addr.crs)
print('Columns:', list(addr.columns))
print(addr.head(3))
print()
print('=== HOSPITALS ===')
print('CRS:', hosp.crs)
print('Columns:', list(hosp.columns))
print(hosp.head(3))
"
=== ADDRESSES ===
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}}
Columns: ['address_id', 'geometry']
address_id geometry
0 A0001 POINT (18.53441 -33.94184)
1 A0002 POINT (18.50296 -33.93623)
2 A0003 POINT (18.53234 -33.95918)
=== HOSPITALS ===
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}}
Columns: ['hospital_id', 'name', 'geometry']
hospital_id name geometry
0 H001 Mowbray Maternity Hospital POINT (18.47483 -33.94921)
1 H002 Life Claremont Hospital POINT (18.46637 -33.98695)
2 H003 Old Somerset Hospital POINT (18.42013 -33.9158)
0:30
Now I understand the data. Both are in WGS 84 (EPSG:4326). I need to: 1. For each address, find the nearest hospital by straight-line (Euclidean) distance. 2. Since it's lat/lon, I should compute proper geodesic distances or reproject to a suitable projected CRS for Western Cape, South Africa. Actually, for "straight-line distance" in metres, using WGS 84 lat/lon directly would give degrees, not metres. Better to use a projected CRS suitable for the Western Cape area. The Western Cape is around lat -33 to -34, lon 18 to 23. A good choice would be UTM zone 34S (EPSG:32734) which covers from 18°E to 24°E, or I could use an Equidistant projection or just use the Haversine formula / geodesic distance. Let me think... The easiest robust approach: use `pyproj` to compute geodesic distances, or reproject to a suitable local CRS. Since the area is small enough, I could also use `shapely`'s `distance` method in a projected CRS. Let me use a geographic distance calculation approach. I'll reproject both to a projected CRS (UTM 34S - EPSG:32734) which covers the Western Cape, then compute Euclidean distances in metres. Actually, the simplest approach for "straight-line distance" in metres is to use the `distance` method from geopandas after projecting to a suitable metric CRS. Let me write a script.
0:30
Write
/work/solve.py
content · 64 lines · py
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely import wkb
# Read input data
addresses = gpd.read_parquet('/work/addresses.parquet')
hospitals = gpd.read_parquet('/work/hospitals.parquet')
print(f"Addresses: {len(addresses)} features")
print(f"Hospitals: {len(hospitals)} features")
# Both are in WGS 84 (EPSG:4326). Project to UTM zone 34S (EPSG:32734)
# which covers the Western Cape area and gives distances in metres.
addr_proj = addresses.to_crs('EPSG:32734')
hosp_proj = hospitals.to_crs('EPSG:32734')
# Extract coordinates for efficient nearest-neighbour search
addr_pts = addr_proj.geometry.values # shapely Point objects
hosp_pts = hosp_proj.geometry.values
hosp_names = hospitals['name'].values
hosp_indices = np.arange(len(hospitals))
results = []
for i, addr_pt in enumerate(addr_pts):
# Compute Euclidean distances to all hospitals
# shapely Point has .x and .y attributes
ax = addr_pt.x
ay = addr_pt.y
# Vectorised distance calculation
distances = np.sqrt((np.array([p.x for p in hosp_pts]) - ax)**2 +
(np.array([p.y for p in hosp_pts]) - ay)**2)
nearest_idx = np.argmin(distances)
nearest_dist = distances[nearest_idx]
nearest_name = str(hosp_names[nearest_idx])
results.append({
'address_id': addresses.iloc[i]['address_id'],
'nearest_hospital_name': nearest_name,
'distance_m': float(nearest_dist),
'geometry': addresses.iloc[i].geometry # Keep original WGS84 geometry
})
# Build output GeoDataFrame
out_gdf = gpd.GeoDataFrame(results, crs='EPSG:4326')
out_gdf = out_gdf[['address_id', 'nearest_hospital_name', 'distance_m', 'geometry']]
# Validate
assert out_gdf['address_id'].notna().all()
assert (out_gdf['nearest_hospital_name'] != '').all()
assert out_gdf['nearest_hospital_name'].notna().all()
assert np.isfinite(out_gdf['distance_m']).all()
assert (out_gdf['distance_m'] >= 0).all()
print(f"\nOutput shape: {out_gdf.shape}")
print(out_gdf.head(10))
print(f"\nDistance stats: min={out_gdf['distance_m'].min():.2f}m, max={out_gdf['distance_m'].max():.2f}m, mean={out_gdf['distance_m'].mean():.2f}m")
# Write to GeoPackage
out_gdf.to_file('/work/nearest_hospital.gpkg', layer='nearest_hospital', driver='GPKG')
print("\nWritten to /work/nearest_hospital.gpkg")Wrote /work/solve.py (2336 chars).
0:44
Let me run this to see if it works.
0:44
Bash
cd /work && python solve.py
Addresses: 120 features Hospitals: 37 features Output shape: (120, 4) 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 A0006 ... POINT (18.48709 -33.97114) 6 A0007 ... POINT (18.54461 -33.96394) 7 A0008 ... POINT (18.52017 -33.97329) 8 A0009 ... POINT (18.51179 -33.97298) 9 A0010 ... POINT (18.49876 -33.9682) [10 rows x 4 columns] Distance stats: min=47.02m, max=5956.74m, mean=1841.99m Written to /work/nearest_hospital.gpkg
0:49
The script ran successfully. Let me verify the output by reading the GeoPackage to make sure everything is correct.
0:49
Bash
cd /work && python3 -c "
cd /work && python3 -c "
import geopandas as gpd
gdf = gpd.read_file('/work/nearest_hospital.gpkg')
print('Columns:', list(gdf.columns))
print('Shape:', gdf.shape)
print('CRS:', gdf.crs)
print()
print(gdf.head())
print()
print('address_id nulls:', gdf['address_id'].isna().sum())
print('name empty strings:', (gdf['nearest_hospital_name'] == '').sum())
print('name nulls:', gdf['nearest_hospital_name'].isna().sum())
print('distance_m nulls:', gdf['distance_m'].isna().sum())
print('distance_m < 0:', (gdf['distance_m'] < 0).sum())
print('distance_m infinite:', (gdf['distance_m'] == float('inf')).sum())
print()
import numpy as np
print('All finite:', np.isfinite(gdf['distance_m']).all())
print('All non-negative:', (gdf['distance_m'] >= 0).all())
"
Columns: ['address_id', 'nearest_hospital_name', 'distance_m', 'geometry'] Shape: (120, 4) 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] address_id nulls: 0 name empty strings: 0 name nulls: 0 distance_m nulls: 0 distance_m < 0: 0 distance_m infinite: 0 All finite: True All non-negative: True
1:01
All validations pass. The output file `/work/nearest_hospital.gpkg` is ready with: - 120 features (one per input address) - Columns: `address_id`, `nearest_hospital_name`, `distance_m`, `geometry` (original WGS84 geometry) - No nulls, no empty strings, no negative or infinite distances - Distances in metres computed via UTM projection (EPSG:32734) for accurate Euclidean measurement in the Western Cape region
1:01