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:07
Let me start by reading both input files to understand their structure.
0:07
Bash
ls /work/
addresses.parquet hospitals.parquet
0:12
Let me first inspect the two input files to understand their schemas and content.
0:12
Bash
cd /work && python -c "
cd /work && python -c "
import pandas as pd
df = pd.read_parquet('addresses.parquet')
print('=== ADDRESSES ===')
print('Columns:', list(df.columns))
print('Shape:', df.shape)
print(df.dtypes)
print(df.head(3))
print()
"
=== ADDRESSES === Columns: ['address_id', 'geometry'] Shape: (120, 2) address_id str geometry object dtype: object 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...
0:14
Bash
cd /work && python -c "
cd /work && python -c "
import pandas as pd
df = pd.read_parquet('hospitals.parquet')
print('=== HOSPITALS ===')
print('Columns:', list(df.columns))
print('Shape:', df.shape)
print(df.dtypes)
print(df.head(3))
"
=== HOSPITALS === Columns: ['hospital_id', 'name', 'geometry'] Shape: (37, 3) hospital_id str name str geometry object dtype: object 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:32
The geometry columns are stored as WKB (binary) in the parquet files. Let me inspect them more carefully to understand the CRS and structure. Let me read them with geopandas.
0:32
Bash
cd /work && python -c "
cd /work && python -c "
import geopandas as gpd
addr = gpd.read_parquet('addresses.parquet')
hosp = gpd.read_parquet('hospitals.parquet')
print('=== ADDRESSES ===')
print('CRS:', addr.crs)
print('Geometry type:', addr.geometry.dtype)
print(addr.head(2))
print()
print('=== HOSPITALS ===')
print('CRS:', hosp.crs)
print('Geometry type:', hosp.geometry.dtype)
print(hosp.head(2))
print()
print('Bounds addr:', addr.total_bounds)
print('Bounds hosp:', hosp.total_bounds)
"
=== 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}}
Geometry type: geometry
address_id geometry
0 A0001 POINT (18.53441 -33.94184)
1 A0002 POINT (18.50296 -33.93623)
=== 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}}
Geometry type: 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)
Bounds addr: [ 18.4008295 -33.99916655 18.54986077 -33.90055752]
Bounds hosp: [ 18.41153087 -33.99823 18.53345868 -33.89997785]
0:59
Now I have a clear picture: 1. Both datasets are in WGS 84 (EPSG:4326) with point geometries 2. Addresses have columns: `address_id`, `geometry` 3. Hospitals have columns: `hospital_id`, `name`, `geometry` 4. Both are in the Western Cape area (Cape Town region) I need to: 1. For each address, find the nearest hospital by straight-line (Euclidean) distance 2. Since they're in degrees, I should use a metric CRS for accurate distance calculation. A good choice for the Western Cape / Cape Town area would be a UTM zone (Zone 34S, EPSG:32734) or a South African CRS like EPSG:32734 (WGS 84 / UTM zone 34S). Actually, for straight-line distance on a sphere, I could also use the haversine formula, but using a projected CRS is cleaner with geopandas. Let me use the WGS 84 / UTM zone 34S (EPSG:32734) which covers the Western Cape area. Or actually, let me check... Cape Town is around 18.4°E longitude which is in UTM zone 34S (18°E to 24°E). So EPSG:32734 is correct. I'll write the script.
0:59
Write
/work/solve.py
content · 71 lines · py
import geopandas as gpd
import pandas as pd
import numpy as np
# ---------------------------------------------------------------------------
# 1. Read input data
# ---------------------------------------------------------------------------
addr = gpd.read_parquet("/work/addresses.parquet")
hosp = gpd.read_parquet("/work/hospitals.parquet")
# Both are in EPSG:4326 (WGS 84) from inspection
assert addr.crs is not None and hosp.crs is not None, "CRS missing"
# ---------------------------------------------------------------------------
# 2. Reproject to a metric CRS so distances are meaningful (metres)
# Cape Town / Western Cape falls in UTM zone 34S (EPSG:32734)
# ---------------------------------------------------------------------------
UTM_CRS = "EPSG:32734"
addr_utm = addr.to_crs(UTM_CRS)
hosp_utm = hosp.to_crs(UTM_CRS)
# ---------------------------------------------------------------------------
# 3. For each address, find the nearest hospital
# ---------------------------------------------------------------------------
# Extract coordinate arrays for fast vectorised distance computation
addr_coords = np.array([(p.x, p.y) for p in addr_utm.geometry])
hosp_coords = np.array([(p.x, p.y) for p in hosp_utm.geometry])
hosp_names = hosp_utm["name"].values
n_addr = len(addr_coords)
n_hosp = len(hosp_coords)
# Build result arrays
nearest_hospital_names = np.empty(n_addr, dtype=object)
nearest_distances = np.empty(n_addr, dtype=np.float64)
# For each address, compute squared Euclidean distance to every hospital
# (avoid sqrt until we pick the winner)
for i in range(n_addr):
dx = addr_coords[i, 0] - hosp_coords[:, 0]
dy = addr_coords[i, 1] - hosp_coords[:, 1]
dist_sq = dx * dx + dy * dy
idx = np.argmin(dist_sq)
nearest_hospital_names[i] = hosp_names[idx]
nearest_distances[i] = np.sqrt(dist_sq[idx]) # metres
# ---------------------------------------------------------------------------
# 4. Build output GeoDataFrame (keep original geometry + columns)
# ---------------------------------------------------------------------------
out = addr[["address_id", "geometry"]].copy()
out["nearest_hospital_name"] = nearest_hospital_names
out["distance_m"] = nearest_distances
# Sanity checks requested by the prompt
assert out["nearest_hospital_name"].notna().all(), "Null hospital name detected"
assert (out["nearest_hospital_name"].str.len() > 0).all(), "Empty hospital name detected"
assert out["distance_m"].notna().all(), "Null distance detected"
assert np.isfinite(out["distance_m"]).all(), "Non-finite distance detected"
assert (out["distance_m"] >= 0).all(), "Negative distance detected"
print(f"Addresses processed: {len(out)}")
print(f"Distances: min={out['distance_m'].min():.2f} m, "
f"max={out['distance_m'].max():.2f} m, "
f"mean={out['distance_m'].mean():.2f} m")
print(f"Unique hospitals referenced: {out['nearest_hospital_name'].nunique()}")
# ---------------------------------------------------------------------------
# 5. Write to GeoPackage
# ---------------------------------------------------------------------------
out.to_file("/work/nearest_hospital.gpkg", layer="nearest_hospital", driver="GPKG")
print("Written to /work/nearest_hospital.gpkg")Wrote /work/solve.py (3236 chars).
1:00
Let me run the script to see if it works.
1:00
Bash
cd /work && python solve.py
Addresses processed: 120 Distances: min=47.02 m, max=5956.74 m, mean=1841.99 m Unique hospitals referenced: 20 Written to /work/nearest_hospital.gpkg
1:07
The script ran successfully. Let me verify the output file exists and has the right columns.
1:07
Bash
cd /work && python -c "
cd /work && python -c "
import geopandas as gpd
gdf = gpd.read_file('nearest_hospital.gpkg')
print('Columns:', list(gdf.columns))
print('Shape:', gdf.shape)
print('CRS:', gdf.crs)
print(gdf.head(5))
print()
# Verify specific requirements
print('address_id dtype:', gdf['address_id'].dtype)
print('nearest_hospital_name dtype:', gdf['nearest_hospital_name'].dtype)
print('distance_m dtype:', gdf['distance_m'].dtype)
print('All names non-empty:', (gdf['nearest_hospital_name'].str.len() > 0).all())
print('All distances >= 0:', (gdf['distance_m'] >= 0).all())
print('All distances finite:', all(gdf['distance_m'].notna()))
print()
# Verify original address geometry preserved
print('Geometry type:', gdf.geometry.dtype)
print(gdf.geometry.head(3))
"
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 dtype: str nearest_hospital_name dtype: str distance_m dtype: float64 All names non-empty: True All distances >= 0: True All distances finite: True Geometry type: geometry 0 POINT (18.53441 -33.94184) 1 POINT (18.50296 -33.93623) 2 POINT (18.53234 -33.95918) Name: geometry, dtype: geometry
1:08
Bash
ls -lh /work/nearest_hospital.gpkg
-rw-r--r-- 1 agent agent 116K Jun 18 16:22 /work/nearest_hospital.gpkg
1:14
Everything checks out: 1. **File format**: GeoPackage (`.gpkg`) ✓ 2. **Rows**: 120 features (same as input addresses) ✓ 3. **Columns**: `address_id`, `nearest_hospital_name`, `distance_m`, `geometry` ✓ 4. **address_id**: preserved from input, str dtype ✓ 5. **nearest_hospital_name**: all non-empty strings, no nulls ✓ 6. **distance_m**: float64, all ≥ 0, all finite ✓ 7. **Geometry**: POINT geometry preserved from original addresses ✓ 8. **CRS**: EPSG:4326 (WGS 84) - original CRS ✓ The output file `nearest_hospital.gpkg` is ready.
1:14