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

import arcgis
from datetime import datetime
from arcgis.features import FeatureLayer
# Connect to GIS
from arcgis.gis import GIS
gis = GIS('home') 

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

# 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)
# 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

# Convert to SpatialDataFrames
clinic = clinic_features.sdf
clinic.head()
AddressCityContact_NameFIDKey_LicenseNumberNamePhone_NumberSHAPEStateTypeZip
0989 University Dr., Ste. 108PontiacDr. Carmen Cardona, Program Director1Case Management; Integrated Treatment; Outpa...# SA0631356SUNRISE TREATMENT CENTER248-481-2531{"x": -9268869.4566, "y": 5259550.238600001, "...MITreatment48342
14120 W. Maple Rd., Ste. 201Bloomfield HillsStacy Gutowsky, Program Director2Alcohol Hwy Safety Education; Community Change...# SA0631330EIP, INC. (Education Intervention Programs)248-693-0336{"x": -9271532.8795, "y": 5242995.079599999, "...MIPrevention48301
23283 Bloomfield Park DrWest BloomfieldKaren Fraiberg, Program Director3Case Management; Community Change, Alternative...# SA0631364THE ART OF HOPE, HEALING & HAPPINESS248-229-4717{"x": -9278211.6833, "y": 5246510.740000002, "...MIBoth48323
32710 W. 12 Mile Rd.BerkleyKathy Smaller, Program Director4Case Management; Early Intervention; Integrate...# SA0631165RECOVERY CONSULTANTS, INC.248-543-1090{"x": -9260261.7618, "y": 5236694.745800003, "...MITreatment48072
43139 W. Huron St. Waterford, MI 48328WaterfordKathy Smaller, Program Director5Case Management; Early Intervention; Integrate...# SA0631166RECOVERY CONSULTANTS, INC.248-738-8400{"x": -9278492.0901, "y": 5256952.089100003, "...MITreatment48328
oakland = oakland_features.sdf
oakland.head()
ACRESCNTY_CODEFIDLABELLAYOUTLINKNAMEPENINSULASHAPESQKMSQMILESShape__AreaShape__LengthTYPEVER
04075.8585591251BG 1portrait12512000011200001lower{"rings": [[[-9250818.94177608, 5289778.006450...16.4945926.3685683.073481e+0722611.540412BG14a
18797.8875411252BG 1portrait12512030011203001lower{"rings": [[[-9264026.53064651, 5290079.339554...35.60417113.7467846.633129e+0735417.658796BG14a
26693.9879151253BG 1landscape12512100011210001lower{"rings": [[[-9270503.93750677, 5288614.193867...27.08990010.4594215.047787e+0729332.713639BG14a
35399.7502071254BG 1portrait12512140011214001lower{"rings": [[[-9269511.1796278, 5286236.4642479...21.8522498.4371624.070574e+0729772.396132BG14a
43719.2863171255BG 2portrait12512220021222002lower{"rings": [[[-9281344.03654134, 5290365.890503...15.0515805.8114212.804448e+0723245.466116BG14a
m = gis.map('Oakland, Michigan')
m
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

from arcgis.features.use_proximity import create_drive_time_areas
# 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')
# Check type of result
type(result)
arcgis.gis.Item
# Share results for public consumption
result.share(everyone=True)
{'results': [{'itemId': 'b31a041a220e40e6a47539a98e2d208f',
   'success': True,
   'notSharedWith': []}]}
# 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.

# Get results as a SpatialDataFrame
driveAll = driveAll_features.sdf
driveAll.head()
AnalysisAreaFacilityOIDFromBreakNameOBJECTIDSHAPEToBreak
0508.239318None2020 - 251{"rings": [[[-9332991.0914, 5246110.1884], [-9...25
1452.186719None1515 - 202{"rings": [[[-9320991.3273, 5233355.6918], [-9...20
2426.955043None1010 - 153{"rings": [[[-9317466.3064, 5285772.2205], [-9...15
3392.647341None55 - 104{"rings": [[[-9293641.3093, 5272963.9065], [-9...10
4345.389051None00 - 55{"rings": [[[-9249366.647, 5243768.4208], [-92...5
m_drive_times = gis.map('Oakland, Michigan')
m_drive_times
m_drive_times.add_layer(result.layers[0])

Service area generation for the visible [5] minutes layer

# Generate service areas
visible_result = create_drive_time_areas(clinic_layer, 
                                         [5], 
                                         output_name="DriveTimeToClinics_5mins"+str(datetime.now().microsecond), 
                                         overlap_policy='Dissolve')
# Grant 'visible_result' public access
visible_result.share(everyone=True)
{'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

from arcgis.features.analysis import join_features
# 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))
# Grant 'cbg_with_driveTime' public access
cbg_with_driveTime.share(everyone=True)
{'results': [{'itemId': '09dbd34a0c144ef7bc692d284b00784a',
   'success': True,
   'notSharedWith': []}]}
# 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

# Convert it to a SpatialDataFrame
drive_time_cbg_df = drive_times_features.sdf
drive_time_cbg_df.head()
ACRESAnalysisAreaCNTY_CODEFacilityOIDFromBreakJoin_CountLABELLAYOUTLINKNAME...OBJECTIDPENINSULASHAPESQKMSQMILESShape__AreaShape__LengthTYPEToBreakVER
04075.858559508.239318125None202BG 1portrait12512000011200001...1lower{"rings": [[[-9250818.9418, 5289778.0065], [-9...16.4945926.3685683.073481e+0722611.540412BG2514a
18797.887541508.239318125None204BG 1portrait12512030011203001...2lower{"rings": [[[-9264026.5306, 5290079.3396], [-9...35.60417113.7467846.633129e+0735417.658796BG2514a
26693.987915452.186719125None153BG 1landscape12512100011210001...3lower{"rings": [[[-9270503.9375, 5288614.1939], [-9...27.08990010.4594215.047787e+0729332.713639BG2014a
35399.750207426.955043125None103BG 1portrait12512140011214001...4lower{"rings": [[[-9269511.1796, 5286236.4643], [-9...21.8522498.4371624.070574e+0729772.396132BG1514a
43719.286317452.186719125None153BG 2portrait12512220021222002...5lower{"rings": [[[-9281344.0365, 5290365.8905], [-9...15.0515805.8114212.804448e+0723245.466116BG2014a

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)
from arcgis.features.enrich_data import enrich_layer
# 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))
# Check type of 'enriched_drive_time_cbg'
type(enriched_drive_time_cbg)
arcgis.gis.Item
# Grant 'enriched_drive_time_cbg' public access
enriched_drive_time_cbg.share(everyone=True)
{'results': [{'itemId': '024a0e8ac00047899b0dc7677cdf0e30',
   'success': True,
   'notSharedWith': []}]}
# 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
# Convert it to a SpatialDataFrame
enriched_df = enriched_driveTimes_features.sdf
enriched_df.head()
ACRESAnalysisAreaBACHDEG_CYCNTY_CODEENRICH_FIDFacilityOIDFromBreakHasDataIDJoin_Count...Shape__LengthTOTPOP_CYTYPEToBreakVERX8002_XaggregationMethodapportionmentConfidencepopulationToPolygonSizeRatingsourceCountry
04075.858559508.239318981251None20102...22611.540412862BG2514a1535591BlockApportionment:US.BlockGroups2.5762.191US
18797.887541508.2393182581252None20114...35417.6587961617BG2514a3183423BlockApportionment:US.BlockGroups2.5762.191US
26693.987915452.186719741253None15123...29332.713639653BG2014a1323720BlockApportionment:US.BlockGroups2.5762.191US
35399.750207426.9550439091254None10133...29772.3961324224BG1514a8447741BlockApportionment:US.BlockGroups2.5762.191US
43719.286317452.1867191241255None15143...23245.4661161691BG2014a2553739BlockApportionment:US.BlockGroups2.5762.191US

5 rows × 31 columns

# Print column names for ease of use during visualization
enriched_df.columns
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

# Create map of Oakland county
map1 = gis.map('Oakland, Michigan', zoomlevel=10)
map1
map1.add_layer(enriched_driveTimes_layer, {'renderer':'ClassedSizeRenderer',
                                           'field_name':'X8002_X',
                                           'opacity':0.6})
map1.add_layer(visible_result.layers[0], {'opacity':0.6})
map1.add_layer(clinic_layer, {'opacity':0.75})

Similar maps can be generated for the following variables too

  1. Travel Time (ToBreak)
  2. Population (TOTPOP_CY)
  3. Citizens with Bachelor's Degree (BACHDEG_CY)

Conclusion

As we can see from the map, clinics are concentrated in a certain part of the county. A few counties have relatively larger population, as compared to their health care spending and population count with Bachelor's degrees and don't have clinics very close by (drime time greater than 10 mins). This subtly hints at disparity in location of clinics with population/affluence level of the county.

This analysis also sets stage for bringing in data about patients of this epidemic, to do a socio-economic study of this epidemic.

Your browser is no longer supported. Please upgrade your browser for the best experience. See our browser deprecation post for more details.