Travel Distance and Estimated Travel Time to Nearest Stroke Center¶

In this notebook, we will calculate the distance (in miles) from a county's population center to the nearest primary or acute stroke center. This variable will serve as an indicator of geographic access to basic stroke care. We will also calculate estimated travel time, which is often a more useful metric than travel distance alone, as access to stroke care is very time-sensitive. In addition to this, we will also calculate the travel distance and estimated travel time to the nearest comprehensive or thrombectomy-capable stroke center. This variable will serve as an indicator of geographic access to advanced stroke care, such as surgery. We obtain the data for county population centers in 2020 from U.S. Census Bureau (https://www.census.gov/geographies/reference-files/time-series/geo/centers-population.html). To obtain the data for the geographic locations of the primary stroke centers in the tri-state area, we take their addresses and geocode them into latitude/longitude values. The addresses are obtained from the following links:

  • NY: https://www.health.ny.gov/diseases/cardiovascular/stroke/designation/stroke_designated_centers.htm
  • NJ: https://www.nj.gov/health/healthcarequality/health-care-professionals/cardiac-stroke-services/stroke-services/list.shtml
  • CT: https://portal.ct.gov/dph/emergency-medical-services/ems/certified-stroke-centers?language=en_US

Loading County Population Center Data¶

In [4]:
import pandas as pd
In [5]:
ct_pop = pd.read_csv(    # Importing CT population center data
    "CenPop2020_Mean_CO09.txt"
)

nj_pop = pd.read_csv(    # Importing NJ population center data
    "CenPop2020_Mean_CO34.txt"
)

ny_pop = pd.read_csv(    # Importing NY population center data
    "CenPop2020_Mean_CO36.txt"
)


pop_centers = pd.concat([ct_pop,ny_pop, nj_pop], ignore_index=True) # Combining the two datasets

Creating County FIPS Codes¶

In [7]:
pop_centers["fips"] = (
    pop_centers["STATEFP"].astype(str).str.zfill(2)
    + pop_centers["COUNTYFP"].astype(str).str.zfill(3)
)
In [8]:
pop_centers = pop_centers.rename(
    columns={
        "COUNAME": "county",
        "STNAME": "state"
    }
)

Loading Geocoded Basic Stroke Center Data¶

In [10]:
ct_basic = pd.read_csv(     # CT basic stroke center data
    "ct_basic_geocoded.csv"
)

ny_basic = pd.read_csv(   # NY basic stroke center data
    "ny_primary_stroke_centers_geocoded.csv"
)

nj_basic = pd.read_csv(   # NJ basic stroke center data
    "nj_primary_stroke_centers_geocoded.csv"
)

stroke_basic = pd.concat(
    [ct_basic,ny_basic, nj_basic],
    ignore_index=True
) # Combining the two datasets
In [11]:
# Checking coordinate column names
print(stroke_basic.columns)

Index(['name', 'group', 'latitude', 'longitude', 'address'], dtype='object')

Installing geopy to Calculate Distance to Stroke Center¶

In [13]:
pip install geopy

Requirement already satisfied: geopy in c:\users\janec\anaconda3\lib\site-packages (2.4.1)
Requirement already satisfied: geographiclib<3,>=1.52 in c:\users\janec\anaconda3\lib\site-packages (from geopy) (2.1)
Note: you may need to restart the kernel to use updated packages.
In [14]:
from geopy.distance import geodesic

Calculating Distance From Nearest Stroke Center¶

Defining Function to Calculate Distance to Nearest Stroke Center¶

In [17]:
def nearest_stroke_distance(county_lat,
                            county_lon,
                            stroke_df):

    distances = stroke_df.apply(
        lambda row:
        geodesic(
            (county_lat, county_lon),
            (row["latitude"], row["longitude"])
        ).miles,
        axis=1
    )

    return distances.min()

Computing Minimum Distance For Each County¶

In [19]:
pop_centers["nearest_stroke_distance"] = (
    pop_centers.apply(
        lambda row:
        nearest_stroke_distance(
            row["LATITUDE"],
            row["LONGITUDE"],
            stroke_basic
        ),
        axis=1
    )
)

Keeping Only Neccessary Columns¶

In [21]:
distance_df = pop_centers[
    [
        "fips",
        "county",
        "nearest_stroke_distance"
    ]
]
In [22]:
distance_df.head()
Out[22]:
fips county nearest_stroke_distance
0 09001 Fairfield 14.240214
1 09003 Hartford 4.320017
2 09005 Litchfield 5.375039
3 09007 Middlesex 6.514388
4 09009 New Haven 5.965944
In [23]:
print(pop_centers.columns.tolist())

['STATEFP', 'COUNTYFP', 'county', 'state', 'POPULATION', 'LATITUDE', 'LONGITUDE', 'fips', 'nearest_stroke_distance']

Saving as a CSV File¶

In [25]:
distance_df.to_csv(
    "county_stroke_distances.csv",
    index=False
)

Calculating Estimated Travel Time to Nearest Stroke Center¶

Installing OpenRouteService¶

In [28]:
pip install openrouteservice

Requirement already satisfied: openrouteservice in c:\users\janec\anaconda3\lib\site-packages (2.3.3)
Requirement already satisfied: requests>=2.0 in c:\users\janec\anaconda3\lib\site-packages (from openrouteservice) (2.32.3)
Requirement already satisfied: charset-normalizer<4,>=2 in c:\users\janec\anaconda3\lib\site-packages (from requests>=2.0->openrouteservice) (3.3.2)
Requirement already satisfied: idna<4,>=2.5 in c:\users\janec\anaconda3\lib\site-packages (from requests>=2.0->openrouteservice) (3.7)
Requirement already satisfied: urllib3<3,>=1.21.1 in c:\users\janec\anaconda3\lib\site-packages (from requests>=2.0->openrouteservice) (2.2.3)
Requirement already satisfied: certifi>=2017.4.17 in c:\users\janec\anaconda3\lib\site-packages (from requests>=2.0->openrouteservice) (2026.6.17)
Note: you may need to restart the kernel to use updated packages.

Connecting to the API¶

In [30]:
import openrouteservice

client = openrouteservice.Client(
    key="eyJvcmciOiI1YjNjZTM1OTc4NTExMTAwMDFjZjYyNDgiLCJpZCI6IjY2MmVhMDc4MGJlOTQyM2VhOTNmZTVkZGM5NzQ0YmY2IiwiaCI6Im11cm11cjY0In0="
)
In [31]:
import numpy as np

Creating County Population Center Coordinates List¶

In [33]:
county_coords = [
    [lon, lat]
    for lat, lon in zip(
        pop_centers["LATITUDE"],
        pop_centers["LONGITUDE"]
    )
]

Creating Basic Stroke Center Coordinates List¶

In [35]:
hospital_coords = [
    [lon, lat]
    for lat, lon in zip(
        stroke_basic["latitude"],
        stroke_basic["longitude"]
    )
]

Initializing an Array of Minimum Travel Times¶

In [37]:
# Creating an array that will store the minimum drive time for each county
# Starting by setting all values to infinity so any real travel time found will be smaller

nearest_minutes = np.full(
    len(pop_centers),
    np.inf
)

Looping Through Hospitals in Groups of 20¶

Since OpenRouteSource has a limit of 3500 routes, we must split the hospitals into chunks.

In [40]:
chunk_size = 20
In [41]:
for i in range(0, len(hospital_coords), chunk_size):


    # Selecting current group of hospitals
    
    hospital_chunk = hospital_coords[i:i+chunk_size]

    # Combining county population center coordinates with hospital coordinates in this chunk

    locations = county_coords + hospital_chunk

    # Defining sources as starting points (county population centers)

    sources = list(range(len(county_coords)))

    # Defining destinations as ending points (hospitals)

    destinations = list(
        range(
            len(county_coords),
            len(locations)
        )
    )

    # Requesting travel time matrix from OpenRouteService to compute driving times from every county to every hospital in the current chunk

    matrix = client.distance_matrix(
        locations=locations,
        profile='driving-car',
        sources=sources,
        destinations=destinations,
        metrics=['duration']
    )

    # Skip chunk if no travel times are returned

    if matrix["durations"] is None:
        print("Failed chunk")
        continue

    # Converting travel time matrix from seconds to minutes

    durations = np.array(
        matrix["durations"],
        dtype=float
    ) / 60

    # Finding minimum travel time for each county 
    
    chunk_min = durations.min(axis=1)

    # Comparing minimum travel times from this chunk with best times found in previous chunks -> keep whichever is smaller

    nearest_minutes = np.minimum(
        nearest_minutes,
        chunk_min
    )

    # Adding final minimum drive times to dataframe

pop_centers["drive_time_min"] = nearest_minutes
In [42]:
print(
    pop_centers[
        ["county", "fips", "drive_time_min"]
    ].head()
)

       county   fips  drive_time_min
0   Fairfield  09001       34.063000
1    Hartford  09003       14.668000
2  Litchfield  09005        9.997500
3   Middlesex  09007       19.046000
4   New Haven  09009       18.025333
In [43]:
pop_centers[
    pop_centers["county"].isin(
        ["Suffolk", "Nassau", "Kings", "Queens", "New York", "Westchester"]
    )
][["county","drive_time_min"]]
Out[43]:
county drive_time_min
31 Kings 2.752833
37 Nassau 10.227667
38 New York 10.068000
48 Queens 6.769167
59 Suffolk 21.352167
67 Westchester 16.221833
In [44]:
pop_centers.head()
Out[44]:
STATEFP COUNTYFP county state POPULATION LATITUDE LONGITUDE fips nearest_stroke_distance drive_time_min
0 9 1 Fairfield Connecticut 957419 41.206021 -73.367214 09001 14.240214 34.063000
1 9 3 Hartford Connecticut 899498 41.761165 -72.717868 09003 4.320017 14.668000
2 9 5 Litchfield Connecticut 185186 41.726958 -73.190432 09005 5.375039 9.997500
3 9 7 Middlesex Connecticut 164245 41.479018 -72.571631 09007 6.514388 19.046000
4 9 9 New Haven Connecticut 864835 41.396593 -72.942327 09009 5.965944 18.025333

Checking for Missing Values¶

In [46]:
pop_centers.loc[
    pop_centers["drive_time_min"].isna(),
    ["county", "state", "nearest_stroke_distance"]
]
Out[46]:
county state nearest_stroke_distance
23 Essex New York 59.472898
28 Hamilton New York 42.745947

Imputing Missing Values¶

Since Essex County and Hamilton County are both deep in the Adirondacks, OpenRouteService may be struggling to connect very remote points to the road network. Thus, we can impute these values assuming an average speed of 45 mph. We use the following formula:

$$ \text{Time (minutes)} = \frac{\text{Distance (miles)}}{45} \times 60 $$
In [50]:
# Filling in missing travel times using straight-line distance and an assumed average driving speed of 45 mph

pop_centers.loc[
    pop_centers["drive_time_min"].isna(),
    "drive_time_min"
] = (
    pop_centers.loc[
        pop_centers["drive_time_min"].isna(),
        "nearest_stroke_distance"
    ]
    / 45
    * 60
)
In [51]:
# Verifying that there are no missing values remaining
print(
    pop_centers["drive_time_min"].isna().sum()
)

0
In [52]:
# Checking imputed counties
pop_centers.loc[
    pop_centers["county"].isin(["Essex", "Hamilton"]),
    [
        "county",
        "state",
        "nearest_stroke_distance",
        "drive_time_min"
    ]
]
Out[52]:
county state nearest_stroke_distance drive_time_min
23 Essex New York 59.472898 79.297198
28 Hamilton New York 42.745947 56.994596
76 Essex New Jersey 0.955867 3.235667

Loading Geocoded Advanced Stroke Center Data¶

In [54]:
# Loading all data

ct_advanced = pd.read_csv("ct_advanced_geocoded.csv")
nj_all = pd.read_csv("nj_all_stroke_centers_geocoded.csv")
ny_all = pd.read_csv("ny_all_stroke_centers_geocoded.csv")
In [55]:
# Keeping only "advanced" facilities (Comprehensive and Thrombectomy-Capable)

ny_advanced = ny_all[ # NY
    ny_all["designation"].isin([
        "Comprehensive Stroke Center",
        "Thrombectomy Capable Stroke Center"
    ])
]

nj_advanced = nj_all[ # NJ
    nj_all["designation"].isin([
        "Comprehensive"
    ])
]

# NJ only has comprehensive and primary designation, so we don't need to worry about thrombectomy-capable

ct_advanced = ct_advanced.copy()
In [56]:
# Combining all three states

stroke_advanced = pd.concat(
    [
        ny_advanced,
        nj_advanced,
        ct_advanced
    ],
    ignore_index=True
)
In [57]:
# Checking for NA values 
print(
    stroke_advanced[
        stroke_advanced["latitude"].isna()
        |
        stroke_advanced["longitude"].isna()
    ]
)

                                                 name    designation  \
44                 Our Lady of Lourdes Medical Center  Comprehensive   
49  Robert Wood Johnson University Hospital New Br...  Comprehensive   
52                       Morristown Memorial Hospital  Comprehensive   

                                              address  latitude  longitude  \
44                  1600 Haddon Ave, Camden, NJ 08103       NaN        NaN   
49  1 Robert Wood Johnson Pl, New Brunswick, NJ 08901       NaN        NaN   
52              100 Madison Ave, Morristown, NJ 07960       NaN        NaN   

   group  
44   NaN  
49   NaN  
52   NaN
In [58]:
# Manually entering the latitude and longitude for missing values
stroke_advanced.loc[
    stroke_advanced["name"] == "Our Lady of Lourdes Medical Center",
    ["latitude", "longitude"]
] = [
    39.9259,
    -75.0965
]

stroke_advanced.loc[
    stroke_advanced["name"] ==
    "Robert Wood Johnson University Hospital New Brunswick",
    ["latitude", "longitude"]
] = [
    40.4959,
    -74.4518
]

stroke_advanced.loc[
    stroke_advanced["name"] ==
    "Morristown Memorial Hospital",
    ["latitude", "longitude"]
] = [
    40.7968,
    -74.4815
]

Creating Advanced Stroke Center Coordinates List¶

In [60]:
hospital_coords_advanced = [
    [lon, lat]
    for lat, lon in zip(
        stroke_advanced["latitude"],
        stroke_advanced["longitude"]
    )
]

Calculating Minimum Travel Distance¶

In [62]:
nearest_distances_advanced = []
In [63]:
for county in county_coords:

    distances = [

        geodesic(
            (county[1], county[0]),
            (hospital[1], hospital[0])
        ).miles

        for hospital in hospital_coords_advanced

    ]

    nearest_distances_advanced.append(
        min(distances)
    )
In [64]:
pop_centers[
    "nearest_stroke_distance_advanced"
] = nearest_distances_advanced

Initializing Array of Minimum Travel Time¶

In [66]:
nearest_minutes_advanced = np.full(
    len(county_coords),
    np.inf
)

Looping Through Hospitals in Groups of 20¶

In [68]:
chunk_size = 20
In [69]:
for i in range(0, len(hospital_coords_advanced), chunk_size):


    # Selecting current group of hospitals
    
    hospital_chunk = hospital_coords_advanced[i:i+chunk_size]

    # Combining county population center coordinates with hospital coordinates in this chunk

    locations = county_coords + hospital_chunk

    # Defining sources as starting points (county population centers)

    sources = list(range(len(county_coords)))

    # Defining destinations as ending points (hospitals)

    destinations = list(
        range(
            len(county_coords),
            len(locations)
        )
    )

    # Requesting travel time matrix from OpenRouteService to compute driving times from every county to every hospital in the current chunk

    matrix = client.distance_matrix(
        locations=locations,
        profile='driving-car',
        sources=sources,
        destinations=destinations,
        metrics=['duration']
    )

    # Skip chunk if no travel times are returned

    if matrix["durations"] is None:
        print("Failed chunk")
        continue

    # Converting travel time matrix from seconds to minutes

    durations = np.array(
        matrix["durations"],
        dtype=float
    ) / 60

    # Finding minimum travel time for each county 
    
    chunk_min = durations.min(axis=1)

    # Comparing minimum travel times from this chunk with best times found in previous chunks -> keep whichever is smaller

    nearest_minutes_advanced = np.minimum(
        nearest_minutes_advanced,
        chunk_min
    )

    # Adding final minimum drive times to dataframe

pop_centers["drive_time_advanced"] = nearest_minutes_advanced
In [70]:
pop_centers.head()
Out[70]:
STATEFP COUNTYFP county state POPULATION LATITUDE LONGITUDE fips nearest_stroke_distance drive_time_min nearest_stroke_distance_advanced drive_time_advanced
0 9 1 Fairfield Connecticut 957419 41.206021 -73.367214 09001 14.240214 34.063000 7.162505 19.296333
1 9 3 Hartford Connecticut 899498 41.761165 -72.717868 09003 4.320017 14.668000 1.359024 6.081167
2 9 5 Litchfield Connecticut 185186 41.726958 -73.190432 09005 5.375039 9.997500 25.652852 48.925167
3 9 7 Middlesex Connecticut 164245 41.479018 -72.571631 09007 6.514388 19.046000 19.785476 41.052000
4 9 9 New Haven Connecticut 864835 41.396593 -72.942327 09009 5.965944 18.025333 6.357835 19.217333

Again, we must impute the values for Essex and Hamilton county, as they are also missing for advanced stroke centers as well. We use the same formula as the one that we used to impute the values for the basic stroke center data.

In [72]:
# Filling in missing travel times using straight-line distance and an assumed average driving speed of 45 mph

pop_centers.loc[
    pop_centers["drive_time_advanced"].isna(),
    "drive_time_advanced"
] = (
    pop_centers.loc[
        pop_centers["drive_time_advanced"].isna(),
        "nearest_stroke_distance_advanced"
    ]
    / 45
    * 60
)
In [73]:
# Verifying that there are no missing values remaining
print(
    pop_centers["drive_time_advanced"].isna().sum()
)

0

Saving as a CSV file¶

In [75]:
distance_df = pop_centers[
    [
        "fips",
        "county",
        "state",
        "drive_time_min",
        "drive_time_advanced",
        "nearest_stroke_distance",
        "nearest_stroke_distance_advanced"
    ]
]

distance_df.to_csv(
    "geographic_stroke_accessibility.csv",
    index=False
)