Skip To Content ArcGIS for Developers Sign In Dashboard
Download the samples Try it live

Drive time analysis for opioid epidemic

This notebook performs drive time analysis for clinics of opioid epidemic treatment and/or prevention centers in the county of Oakland, Michigan.

Problem

Drug overdoses are a leading cause of injury death in the United States, resulting in approximately 52,000 deaths in 2015. In the same year, the overdose death rate in non metropolitan areas of the country (17 per 100,000) surpassed the metropolitan area rate (16.2 per 100,000), which historically wasn't the case.

While metro areas have more people, and consequently, a greater incidence of drug use and overdose, the one stark disadvantage in the non metro areas is the comparatively limited access to opioid epidemic prevention/treatment facilities and transit.

Data Driven Approach

Here, we demonstrate an analysis of how several parts of Oakland county in Michigan are accessible from their opioid epidemic prevention and/or treatment centers within 5, 10, 15, 20 and 25 minutes of drive time. The areas of the county identified within those time bound service areas are then enriched with factors such as

Health insurance spending

Number of people with a Bachelor's Degree

Total Population in 2017

This provides a better understanding of the population served or under-served by these treatment facilities and the plausible contributing factors.

Steps of analysis

  1. Get data for clinics and census block groups boundary (CBG) in Oakland County
  2. Generate areas served based on time taken to drive to the clinics
  3. Spatially join the service areas generated with the CBG to obtain drive time to closest clinic for each CBG.
  4. Enrich each CBG with additional data to estimate influencing factors.
  5. Map it!

Let's get started by first importing necessary modules and connecting to our GIS

In [1]:
import arcgis
from datetime import datetime
from arcgis.features import FeatureLayer
In [6]:
# Connect to GIS
from arcgis.gis import GIS
gis = GIS(profile="your_online_profile") 

1. Reading in clinics and Oakland county census block group (CBG) boundaries

In [2]:
# Layer URLS for the Opioid Epidemic Clinics in Oakland County and Census Block Groups of Oakland County
clinic_url = 'https://services8.arcgis.com/TWfU0bgvgUkCgHLa/arcgis/rest/services/Oakland_County_Licensed_Program_Listings/FeatureServer/0/'
oakland_cbg = 'https://services.arcgis.com/bkrWlSKcjUDFDtgw/arcgis/rest/services/Oakland_2010_census_bg/FeatureServer/0/'

# create feature layers
clinic_layer = FeatureLayer(clinic_url)
oakland_layer = FeatureLayer(oakland_cbg)
In [3]:
# Query all clinics and census block groups as FeatureSets
clinic_features = clinic_layer.query()
print('Total number of rows in the clinic dataset: ' + str(len(clinic_features.features)))

oakland_features = oakland_layer.query()
print('Total number of rows in the oakland dataset: ' + str(len(oakland_features.features)))
Total number of rows in the clinic dataset: 122
Total number of rows in the oakland dataset: 934

We notice that Oakland County, MI has 122 clinics and 934 CBG

In [4]:
# Convert to SpatialDataFrames
clinic = clinic_features.sdf
clinic.head()
Out[4]:
Address City Contact_Name FID Key_ LicenseNumber Name Phone_Number SHAPE State Type Zip
0 989 University Dr., Ste. 108 Pontiac Dr. Carmen Cardona, Program Director 1 Case Management; Integrated Treatment; Outpa... # SA0631356 SUNRISE TREATMENT CENTER 248-481-2531 {"x": -9268869.4566, "y": 5259550.238600001, "... MI Treatment 48342
1 4120 W. Maple Rd., Ste. 201 Bloomfield Hills Stacy Gutowsky, Program Director 2 Alcohol Hwy Safety Education; Community Change... # SA0631330 EIP, INC. (Education Intervention Programs) 248-693-0336 {"x": -9271532.8795, "y": 5242995.079599999, "... MI Prevention 48301
2 3283 Bloomfield Park Dr West Bloomfield Karen Fraiberg, Program Director 3 Case Management; Community Change, Alternative... # SA0631364 THE ART OF HOPE, HEALING & HAPPINESS 248-229-4717 {"x": -9278211.6833, "y": 5246510.740000002, "... MI Both 48323
3 2710 W. 12 Mile Rd. Berkley Kathy Smaller, Program Director 4 Case Management; Early Intervention; Integrate... # SA0631165 RECOVERY CONSULTANTS, INC. 248-543-1090 {"x": -9260261.7618, "y": 5236694.745800003, "... MI Treatment 48072
4 3139 W. Huron St. Waterford, MI 48328 Waterford Kathy Smaller, Program Director 5 Case Management; Early Intervention; Integrate... # SA0631166 RECOVERY CONSULTANTS, INC. 248-738-8400 {"x": -9278492.0901, "y": 5256952.089100003, "... MI Treatment 48328
In [5]:
oakland = oakland_features.sdf
oakland.head()
Out[5]:
ACRES CNTY_CODE FID LABEL LAYOUT LINK NAME PENINSULA SHAPE SQKM SQMILES Shape__Area Shape__Length TYPE VER
0 4075.858559 125 1 BG 1 portrait 1251200001 1200001 lower {"rings": [[[-9250818.94177608, 5289778.006450... 16.494592 6.368568 3.073481e+07 22611.540412 BG 14a
1 8797.887541 125 2 BG 1 portrait 1251203001 1203001 lower {"rings": [[[-9264026.53064651, 5290079.339554... 35.604171 13.746784 6.633129e+07 35417.658796 BG 14a
2 6693.987915 125 3 BG 1 landscape 1251210001 1210001 lower {"rings": [[[-9270503.93750677, 5288614.193867... 27.089900 10.459421 5.047787e+07 29332.713639 BG 14a
3 5399.750207 125 4 BG 1 portrait 1251214001 1214001 lower {"rings": [[[-9269511.1796278, 5286236.4642479... 21.852249 8.437162 4.070574e+07 29772.396132 BG 14a
4 3719.286317 125 5 BG 2 portrait 1251222002 1222002 lower {"rings": [[[-9281344.03654134, 5290365.890503... 15.051580 5.811421 2.804448e+07 23245.466116 BG 14a
In [19]:
m = gis.map('Oakland, Michigan')
m
In [8]:
m.add_layer(oakland_layer)
m.add_layer(clinic_layer)

2. Drive time analysis

Service areas around the clinics need to be generated twice for two different sets of drive time inputs:

  1. [5, 10, 15, 20, 25] that is needed to get the drive time required from each CBG to the closest clinic. While adding more/less intervals is certainly an option, the choice of time intervals here is meant to span the range (closest-farthest clinic) as well as assume a range of safety, i.e. time within which a person can be driven/drive to a clinic to avoid delays.

  2. [5] that is needed purely for visualization purposes (as seen towards the end of the notebook). Feel free to modify this function parameter based on how you would like to visualize the map.

Service area generation for the [5, 10, 15, 20, 25] minutes layer

In [9]:
from arcgis.features.use_proximity import create_drive_time_areas
In [10]:
# Generate service areas
result = create_drive_time_areas(clinic_layer, [5, 10, 15, 20, 25], 
                                 output_name="DriveTimeToClinics_All_2"+ str(datetime.now().microsecond), 
                                 overlap_policy='Dissolve')
In [11]:
# Check type of result
type(result)
Out[11]:
arcgis.gis.Item
In [12]:
# Share results for public consumption
result.share(everyone=True)
Out[12]:
{'results': [{'itemId': 'b31a041a220e40e6a47539a98e2d208f',
   'success': True,
   'notSharedWith': []}]}
In [13]:
# Convert it to a FeatureSet
drive_times_all_layer = result.layers[0]
driveAll_features = drive_times_all_layer.query()
print('Total number of rows in the service area (upto 25 mins) dataset: '+ str(len(driveAll_features.features)))
Total number of rows in the service area (upto 25 mins) dataset: 5

This implies that only 5 clinics in this county can be driven to within 5 minutes.

In [14]:
# Get results as a SpatialDataFrame
driveAll = driveAll_features.sdf
driveAll.head()
Out[14]:
AnalysisArea FacilityOID FromBreak Name OBJECTID SHAPE ToBreak
0 508.239318 None 20 20 - 25 1 {"rings": [[[-9332991.0914, 5246110.1884], [-9... 25
1 452.186719 None 15 15 - 20 2 {"rings": [[[-9320991.3273, 5233355.6918], [-9... 20
2 426.955043 None 10 10 - 15 3 {"rings": [[[-9317466.3064, 5285772.2205], [-9... 15
3 392.647341 None 5 5 - 10 4 {"rings": [[[-9293641.3093, 5272963.9065], [-9... 10
4 345.389051 None 0 0 - 5 5 {"rings": [[[-9249366.647, 5243768.4208], [-92... 5
In [20]:
m_drive_times = gis.map('Oakland, Michigan')
m_drive_times
In [16]:
m_drive_times.add_layer(result.layers[0])

Service area generation for the visible [5] minutes layer

In [29]:
# Generate service areas
visible_result = create_drive_time_areas(clinic_layer, 
                                         [5], 
                                         output_name="DriveTimeToClinics_5mins"+str(datetime.now().microsecond), 
                                         overlap_policy='Dissolve')
In [30]:
# Grant 'visible_result' public access
visible_result.share(everyone=True)
Out[30]:
{'results': [{'itemId': 'b99b686a85a749e09ca46c480c2b787f',
   'success': True,
   'notSharedWith': []}]}

3. Spatially join the drive_times_all_layer layer with oakland_cbg layer to obtain drive times for each census block group

In [31]:
from arcgis.features.analysis import join_features
In [32]:
# Perform spatial join between CBG layer and the service areas created for all time durations
cbg_with_driveTime = join_features(oakland_cbg, 
                                   drive_times_all_layer, 
                                   spatial_relationship='Intersects', 
                                   output_name='Oakland_CBG_allDriveTimes'+str(datetime.now().microsecond))
In [33]:
# Grant 'cbg_with_driveTime' public access
cbg_with_driveTime.share(everyone=True)
Out[33]:
{'results': [{'itemId': '09dbd34a0c144ef7bc692d284b00784a',
   'success': True,
   'notSharedWith': []}]}
In [34]:
# Converting it to a FeatureSet
drive_times_cbg = cbg_with_driveTime.layers[0]
drive_times_features = drive_times_cbg.query()
print('Total number of rows in the drive time (upto 25 mins) dataset: '+ str(len(drive_times_features.features)))
Total number of rows in the drive time (upto 25 mins) dataset: 934

This implies that the farthest clinic in the county is within a 25 minute drive time window

In [37]:
# Convert it to a SpatialDataFrame
drive_time_cbg_df = drive_times_features.sdf
drive_time_cbg_df.head()
Out[37]:
ACRES AnalysisArea CNTY_CODE FacilityOID FromBreak Join_Count LABEL LAYOUT LINK NAME ... OBJECTID PENINSULA SHAPE SQKM SQMILES Shape__Area Shape__Length TYPE ToBreak VER
0 4075.858559 508.239318 125 None 20 2 BG 1 portrait 1251200001 1200001 ... 1 lower {"rings": [[[-9250818.9418, 5289778.0065], [-9... 16.494592 6.368568 3.073481e+07 22611.540412 BG 25 14a
1 8797.887541 508.239318 125 None 20 4 BG 1 portrait 1251203001 1203001 ... 2 lower {"rings": [[[-9264026.5306, 5290079.3396], [-9... 35.604171 13.746784 6.633129e+07 35417.658796 BG 25 14a
2 6693.987915 452.186719 125 None 15 3 BG 1 landscape 1251210001 1210001 ... 3 lower {"rings": [[[-9270503.9375, 5288614.1939], [-9... 27.089900 10.459421 5.047787e+07 29332.713639 BG 20 14a
3 5399.750207 426.955043 125 None 10 3 BG 1 portrait 1251214001 1214001 ... 4 lower {"rings": [[[-9269511.1796, 5286236.4643], [-9... 21.852249 8.437162 4.070574e+07 29772.396132 BG 15 14a
4 3719.286317 452.186719 125 None 15 3 BG 2 portrait 1251222002 1222002 ... 5 lower {"rings": [[[-9281344.0365, 5290365.8905], [-9... 15.051580 5.811421 2.804448e+07 23245.466116 BG 20 14a

5 rows × 21 columns

4. Enrich census block groups with suitable data

The 3 geo-enrichment factors that we will be including in this study are:

  1. 2017 Total Population (Analysis Variable - TOTPOP_CY)
  2. 2017 Health Insurance (Analysis Variable - X8002_X)
  3. 2017 Education : Bachelor's Degree (Analysis Variable - BACHDEG_CY)
In [38]:
from arcgis.features.enrich_data import enrich_layer
In [39]:
# Perform data enrichment
enriched_drive_time_cbg = enrich_layer(drive_times_cbg, 
                                       analysis_variables=['TOTPOP_CY', 'X8002_X', 'BACHDEG_CY'], 
                                       output_name="Enriched_DriveTime_CBG"+str(datetime.now().microsecond))
In [40]:
# Check type of 'enriched_drive_time_cbg'
type(enriched_drive_time_cbg)
Out[40]:
arcgis.gis.Item
In [41]:
# Grant 'enriched_drive_time_cbg' public access
enriched_drive_time_cbg.share(everyone=True)
Out[41]:
{'results': [{'itemId': '024a0e8ac00047899b0dc7677cdf0e30',
   'success': True,
   'notSharedWith': []}]}
In [42]:
# Convert it to a Feature Set
enriched_driveTimes_layer = enriched_drive_time_cbg.layers[0]
enriched_driveTimes_features = enriched_driveTimes_layer.query()
print('Total number of rows in the drive time (upto 25 mins) dataset: '+ str(len(enriched_driveTimes_features.features)))
Total number of rows in the drive time (upto 25 mins) dataset: 
934
In [43]:
# Convert it to a SpatialDataFrame
enriched_df = enriched_driveTimes_features.sdf
enriched_df.head()
Out[43]:
ACRES AnalysisArea BACHDEG_CY CNTY_CODE ENRICH_FID FacilityOID FromBreak HasData ID Join_Count ... Shape__Length TOTPOP_CY TYPE ToBreak VER X8002_X aggregationMethod apportionmentConfidence populationToPolygonSizeRating sourceCountry
0 4075.858559 508.239318 98 125 1 None 20 1 0 2 ... 22611.540412 862 BG 25 14a 1535591 BlockApportionment:US.BlockGroups 2.576 2.191 US
1 8797.887541 508.239318 258 125 2 None 20 1 1 4 ... 35417.658796 1617 BG 25 14a 3183423 BlockApportionment:US.BlockGroups 2.576 2.191 US
2 6693.987915 452.186719 74 125 3 None 15 1 2 3 ... 29332.713639 653 BG 20 14a 1323720 BlockApportionment:US.BlockGroups 2.576 2.191 US
3 5399.750207 426.955043 909 125 4 None 10 1 3 3 ... 29772.396132 4224 BG 15 14a 8447741 BlockApportionment:US.BlockGroups 2.576 2.191 US
4 3719.286317 452.186719 124 125 5 None 15 1 4 3 ... 23245.466116 1691 BG 20 14a 2553739 BlockApportionment:US.BlockGroups 2.576 2.191 US

5 rows × 31 columns

In [46]:
# Print column names for ease of use during visualization
enriched_df.columns
Out[46]:
Index(['ACRES', 'AnalysisArea', 'BACHDEG_CY', 'CNTY_CODE', 'ENRICH_FID',
       'FacilityOID', 'FromBreak', 'HasData', 'ID', 'Join_Count', 'LABEL',
       'LAYOUT', 'LINK', 'NAME', 'Name_1', 'OBJECTID', 'PENINSULA', 'SHAPE',
       'SQKM', 'SQMILES', 'Shape__Area', 'Shape__Length', 'TOTPOP_CY', 'TYPE',
       'ToBreak', 'VER', 'X8002_X', 'aggregationMethod',
       'apportionmentConfidence', 'populationToPolygonSizeRating',
       'sourceCountry'],
      dtype='object')

5. Visualize results

In [2]:
# Create map of Oakland county
map1 = gis.map('Oakland, Michigan', zoomlevel=10)
map1
Out[2]: