How to Measure Distance to the Nearest Hospital in Python

Problem statement

How far people live from a hospital emergency department is one of the most common questions in health geography, and three choices decide the answer before any distance is computed. Which hospitals count: every licensed facility, or those with an emergency department? Where do people live: at tract centroids, which in rural areas can sit in empty farmland, or at population-weighted centres? And where does the study area end: a state boundary that patients cross every day, or the hospitals they actually use?

This page measures the straight-line distance from every Illinois census tract's centre of population to the nearest hospital with emergency services, using the CMS list of Medicare-certified hospitals geocoded with the Census Bureau's geocoder, and shows how much the answer depends on including hospitals across the state line.

Quick answer

Put hospitals and population centres on the unit sphere, find the nearest hospital with a KD-tree, and convert the chord length to a great-circle distance:

def to_xyz(lat, lon):
    la, lo = np.radians(lat), np.radians(lon)
    return np.column_stack([np.cos(la) * np.cos(lo), np.cos(la) * np.sin(lo), np.sin(la)])

tree = cKDTree(to_xyz(HOSP.lat, HOSP.lon))
chord, idx = tree.query(to_xyz(POP.LATITUDE, POP.LONGITUDE))
POP["km"] = 2 * 6371.0 * np.arcsin(chord / 2)                       # chord on the unit sphere → great-circle km
POP["nearest"] = HOSP["Facility Name"].to_numpy()[idx]
w = POP.POPULATION
print(f"{len(POP):,} tracts, {w.sum():,} people; population-weighted median distance {np.median(np.repeat(POP.km, w)):.1f} km")
for d in (10, 20, 30):
    print(f"  more than {d} km from an emergency department: {w[POP.km > d].sum() / w.sum():.1%}")

3,261 tracts, 12,812,508 people; population-weighted median distance 4.0 km
  more than 10 km from an emergency department: 16.6%
  more than 20 km from an emergency department: 4.1%
  more than 30 km from an emergency department: 0.5%

Which hospitals count, checked geocodes, distance from homes

1. Decide which hospitals count

The CMS file has an Emergency Services flag, and it is not always right. Hospitals with an emergency department report its waiting times to CMS (measure OP_18b), which gives a second source to check the flag against:

il = h[h.State == "IL"]
print(f"Illinois hospitals: {len(il)}; flagged with emergency services: {(il['Emergency Services'] == 'Yes').sum()}")
odd = il[(il["Emergency Services"] == "No") & il["Facility ID"].isin(REPORTS_ED)]
print(f"OP_18b period: {op18['Start Date'].mode()[0]} to {op18['End Date'].mode()[0]}")
print("flagged No but reporting emergency department times:", "; ".join(odd["Facility Name"].str.title() + " (" + odd["City/Town"].str.title() + ")"))
print(f"flagged Yes but psychiatric or children's: {((il['Emergency Services'] == 'Yes') & ~GENERAL[il.index]).sum()}; "
      f"used here: {(HOSP.State == 'IL').sum()} in Illinois, {len(HOSP)} with the neighbouring states")

Illinois hospitals: 193; flagged with emergency services: 175
OP_18b period: 10/01/2024 to 09/30/2025
flagged No but reporting emergency department times: West Suburban Medical Center (Oak Park); Osf Saint Anthony'S Health Center (Alton); Saint James Hospital (Pontiac)
flagged Yes but psychiatric or children's: 5; used here: 173 in Illinois, 710 with the neighbouring states

The three hospitals flagged No reported median emergency department times for hundreds of visits in the year to September 2025, so they had working emergency departments then. Whichever source is out of date, a list of facilities should be checked against a second one before distances are computed from it.

2. Geocode the hospitals and check the matches

The CMS file has addresses but no coordinates. The Census Bureau's batch geocoder takes up to 10,000 addresses per request; some will not match, and those must not be silently dropped:

print(HOSP.how.value_counts().to_string())
print("unmatched in Illinois, placed at their ZIP code:", ", ".join(HOSP[(HOSP.how == "ZIP centre") & (HOSP.State == "IL")]["Facility Name"].str.title().head(5)), "...")
ok = HOSP[HOSP.how == "address"].assign(zlat=lambda d: d["ZIP Code"].str[:5].map(zip_pt.INTPTLAT),
                                        zlon=lambda d: d["ZIP Code"].str[:5].map(zip_pt.INTPTLONG)).dropna(subset=["zlat"])
err = Geod(ellps="WGS84").inv(ok.lon.to_numpy(), ok.lat.to_numpy(), ok.zlon.to_numpy(), ok.zlat.to_numpy())[2] / 1000
print(f"for matched hospitals, the ZIP centre is {np.median(err):.1f} km from the address (median), {np.percentile(err, 90):.1f} km at the 90th percentile")

how
address       625
ZIP centre     85
unmatched in Illinois, placed at their ZIP code: Osf Saint Anthony'S Health Center, Palos Community Hospital, Uchicago Medicine Adventhealth La Grange, Hines Va Medical Center, Swedish Hospital ...
for matched hospitals, the ZIP centre is 2.6 km from the address (median), 6.4 km at the 90th percentile

The matched hospitals show how far off a ZIP centre is likely to be: small against rural distances, not against urban ones. Keep the flag, so the approximate positions can be found again.

3. Measure from where people live

Tract centroids are the centres of the tract's area. The Census Bureau's centres of population are the averages of where the tract's residents live, and in large rural tracts the two can be far apart:

tr = TRACTS.merge(POP[["GEOID", "LATITUDE", "LONGITUDE"]], on="GEOID")
cen = tr.geometry.to_crs(5070).centroid.to_crs(4326)
gap = Geod(ellps="WGS84").inv(cen.x.to_numpy(), cen.y.to_numpy(), tr.LONGITUDE.to_numpy(), tr.LATITUDE.to_numpy())[2] / 1000
print(f"distance between tract centroid and centre of population: median {np.median(gap):.2f} km, "
      f"90th percentile {np.percentile(gap, 90):.1f} km, largest {gap.max():.1f} km")

distance between tract centroid and centre of population: median 0.24 km, 90th percentile 2.2 km, largest 9.7 km

4. Include the hospitals across the border

Illinois residents along the Mississippi, the Ohio and the Indiana line use hospitals in neighbouring states. Measuring only to Illinois hospitals overstates their distances:

il_only = HOSP[HOSP.State == "IL"].reset_index(drop=True)
chord_il, _ = cKDTree(to_xyz(il_only.lat, il_only.lon)).query(to_xyz(POP.LATITUDE, POP.LONGITUDE))
POP["km_il_only"] = 2 * 6371.0 * np.arcsin(chord_il / 2)
moved = POP.km_il_only - POP.km
print(f"tracts whose nearest emergency department is out of state: {(moved > 0).sum()} ({w[moved > 0].sum():,} people)")
print(f"share of people over 30 km away: {w[POP.km_il_only > 30].sum() / w.sum():.2%} with Illinois hospitals only, "
      f"{w[POP.km > 30].sum() / w.sum():.2%} with neighbours'")

tracts whose nearest emergency department is out of state: 57 (212,358 people)
share of people over 30 km away: 0.63% with Illinois hospitals only, 0.54% with neighbours'

How often the nearest is a critical access hospital

The data used on this page

The CMS Hospital General Information file (Medicare-certified hospitals, with an Emergency Services flag), filtered to Illinois and its five neighbouring states; the Census Bureau's batch geocoder for coordinates, with unmatched addresses placed at the internal point of their ZIP code tabulation area from the 2023 Gazetteer; the 2020 centres of population of Illinois tracts; and the 2023 tract boundary file.

import io
import zipfile
from pathlib import Path

import geopandas as gpd
import numpy as np
import pandas as pd
import requests
from pyproj import Geod
from scipy.spatial import cKDTree

Path("data").mkdir(exist_ok=True)
meta = requests.get("https://data.cms.gov/provider-data/api/1/metastore/schemas/dataset/items/xubh-q36u", timeout=60).json()
meta_ed = requests.get("https://data.cms.gov/provider-data/api/1/metastore/schemas/dataset/items/yv7e-xc69", timeout=60).json()
SRC = {"Hospital_General_Information.csv": meta["distribution"][0]["downloadURL"],
       "Timely_and_Effective_Care-Hospital.csv": meta_ed["distribution"][0]["downloadURL"],
       "CenPop2020_Mean_TR17.txt": "https://www2.census.gov/geo/docs/reference/cenpop2020/tract/CenPop2020_Mean_TR17.txt",
       "gaz_zcta.zip": "https://www2.census.gov/geo/docs/maps-data/data/gazetteer/2023_Gazetteer/2023_Gaz_zcta_national.zip",
       "cb_2023_17_tract_500k.zip": "https://www2.census.gov/geo/tiger/GENZ2023/shp/cb_2023_17_tract_500k.zip"}
for name, url in SRC.items():
    if not Path("data", name).exists():
        Path("data", name).write_bytes(requests.get(url, timeout=300).content)

h = pd.read_csv("data/Hospital_General_Information.csv", dtype=str)
h = h[h.State.isin(["IL", "IN", "IA", "MO", "KY", "WI"])]
tec = pd.read_csv("data/Timely_and_Effective_Care-Hospital.csv", dtype=str)
op18 = tec[(tec["Measure ID"] == "OP_18b") & (tec.Score != "Not Available")]    # median time in the emergency department
REPORTS_ED = set(op18["Facility ID"])
h["has_ed"] = (h["Emergency Services"] == "Yes") | h["Facility ID"].isin(REPORTS_ED)
GENERAL = ~h["Hospital Type"].isin(["Psychiatric", "Childrens"])            # not general emergency departments
HOSP = h[h.has_ed & GENERAL].reset_index(drop=True)
G = Path("data/geocoded.csv")
if not G.exists():                                                          # one request for every address
    batch = HOSP[["Facility ID", "Address", "City/Town", "State", "ZIP Code"]].to_csv(index=False, header=False)
    r = requests.post("https://geocoding.geo.census.gov/geocoder/locations/addressbatch", timeout=600,
                      files={"addressFile": ("hospitals.csv", batch)}, data={"benchmark": "Public_AR_Current"})
    G.write_text(r.text)
geo = pd.read_csv(G, header=None, dtype=str, names=["Facility ID", "input", "match", "exact", "matched", "coords", "tiger", "side"])
geo = geo[geo.match == "Match"]
geo[["lon", "lat"]] = geo.coords.str.split(",", expand=True).astype(float)
HOSP = HOSP.merge(geo[["Facility ID", "lat", "lon"]], on="Facility ID", how="left")
z = pd.read_csv(zipfile.ZipFile("data/gaz_zcta.zip").open("2023_Gaz_zcta_national.txt"), sep="\t", dtype={"GEOID": str})
z.columns = z.columns.str.strip()                                           # the last header carries trailing spaces
zip_pt = z.set_index("GEOID")[["INTPTLAT", "INTPTLONG"]]
miss = HOSP.lat.isna()
HOSP.loc[miss, "lat"] = HOSP.loc[miss, "ZIP Code"].str[:5].map(zip_pt.INTPTLAT)
HOSP.loc[miss, "lon"] = HOSP.loc[miss, "ZIP Code"].str[:5].map(zip_pt.INTPTLONG)
HOSP["how"] = np.where(miss, "ZIP centre", "address")
lost = HOSP[HOSP.lat.isna()]                                                # ZIP codes with no tabulation area
HOSP = HOSP.dropna(subset=["lat"]).reset_index(drop=True)

POP = pd.read_csv("data/CenPop2020_Mean_TR17.txt", dtype={"STATEFP": str, "COUNTYFP": str, "TRACTCE": str}, encoding="utf-8-sig")
POP["GEOID"] = POP.STATEFP + POP.COUNTYFP + POP.TRACTCE
POP = POP[POP.POPULATION > 0].reset_index(drop=True)
TRACTS = gpd.read_file("data/cb_2023_17_tract_500k.zip")
print(f"{len(HOSP)} hospitals with emergency services in IL and neighbours ({(HOSP.State == 'IL').sum()} in Illinois); "
      f"{len(POP):,} populated Illinois tracts")
print("dropped, no address match and no ZIP area:", "; ".join(lost["Facility Name"].str.title() + ", " + lost.State))

710 hospitals with emergency services in IL and neighbours (173 in Illinois); 3,261 populated Illinois tracts
dropped, no address match and no ZIP area: The Medical Center At Franklin, KY; Cox Medical Center Branson, MO; University Of Missouri Health Care, MO

Example 1: how often the nearest hospital is a critical access hospital

Critical access hospitals are small rural hospitals with at most 25 beds:

POP["type"] = HOSP["Hospital Type"].to_numpy()[idx]
share = POP.groupby("type").POPULATION.sum() / w.sum()
print(share.sort_values(ascending=False).map("{:.1%}".format).to_string())
far = POP[POP.km > 20]
print(f"of people more than 20 km away, nearest is a critical access hospital for "
      f"{far[far.type == 'Critical Access Hospitals'].POPULATION.sum() / far.POPULATION.sum():.0%}")

type
Acute Care Hospitals                    90.1%
Critical Access Hospitals                8.7%
Acute Care - Veterans Administration     1.2%
of people more than 20 km away, nearest is a critical access hospital for 34%

Example 2: the counties farthest from care

POP["county"] = POP.COUNTYFP
c = POP.groupby("county").apply(lambda d: pd.Series({"people": d.POPULATION.sum(),
                                                     "median_km": np.median(np.repeat(d.km, d.POPULATION))}), include_groups=False)
names = TRACTS.drop_duplicates("COUNTYFP").set_index("COUNTYFP").NAMELSADCO
print(c.sort_values("median_km", ascending=False).head(5).assign(name=lambda d: d.index.map(names)).round(1).to_string())

         people  median_km             name
county                                     
193     13877.0       32.3     White County
123     11742.0       31.9  Marshall County
041     19740.0       30.6   Douglas County
153      5193.0       29.8   Pulaski County
079      9287.0       29.1    Jasper County

Straight-line distance is a floor

Straight-line distance is a floor, not a travel time. Roads, rivers and bridges make real journeys longer, sometimes much longer where a river has few crossings. Drive times from a road network are the better measure where they matter, and the same three choices apply to them: which facilities, measured from where, and across which boundaries. The straight-line version is still the right first step, because it is quick, needs no routing data, and identifies where the more expensive analysis is needed.

The emergency services flag in the CMS file means a hospital has an emergency department; it says nothing about trauma level, stroke or heart attack capability, or obstetrics. A study of access to a particular kind of care needs a list of the facilities that provide it, and the nearest hospital on this map may not be the one an ambulance takes a stroke patient to.

List coverage, geocoding and ZIP fallbacks

• Coverage of the list: CMS lists Medicare-certified hospitals; freestanding emergency departments and some facilities are not in it.
• Geocoding: check the match rate and the unmatched list; large, well-known hospitals fail too, such as Hines VA Medical Center here.
• ZIP fallback: a ZIP tabulation area's internal point can be several kilometres from the hospital in rural ZIPs.
• Population: 2020 centres of population; tracts with no residents are dropped.
• Distances: the chord-to-arc formula is exact on a sphere; the difference from the ellipsoid is well under 1% at these distances.

Internal links

• Accessibility measures explained — distance, time and opportunity measures
• Why do accessibility scores fail at the study area edge? — the border problem in general

FAQ

How do I find the nearest hospital to each census tract in Python?

Put hospital and tract coordinates on the unit sphere, query a scipy cKDTree, and convert the chord distance to kilometres with 2R·asin(c/2).

Where can I get a list of US hospitals with coordinates?

The CMS Hospital General Information file has every Medicare-certified hospital with its address and emergency services flag; geocode the addresses with the Census Bureau's free batch geocoder.

Should I use tract centroids or population centres?

Population centres, published by the Census Bureau for every tract, which place each tract where its residents live.

Why include hospitals in other states?

People near a state line often use hospitals across it; leaving them out overstates distances for border communities.