Drive time analysis for opioid epidemic

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

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

In [1]:
import arcgis
from arcgis.features import FeatureLayer
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 clonic and census block groups as FeatureSets
clinic_features = clinic_layer.query()
print('Total number of rows in the clinic dataset: ')
print(len(clinic_features.features))

oakland_features = oakland_layer.query()
print('Total number of rows in the oakland dataset: ')
print(len(oakland_features.features))
Total number of rows in the clinic dataset: 
122
Total number of rows in the oakland dataset: 
934
In [4]:
# Convert to SpatialDataFrames
clinic = clinic_features.sdf
clinic.head()
Out[4]:
Address City Contact_Name FID Key_ LicenseNumber Name Phone_Number State Type Zip SHAPE
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 MI Treatment 48342 {'x': -9268869.4566, 'y': 5259550.238600001, '...
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 MI Prevention 48301 {'x': -9271532.8795, 'y': 5242995.079599999, '...
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 MI Both 48323 {'x': -9278211.6833, 'y': 5246510.740000002, '...
3 2710 W. 12 Mile Rd. Berkley Kathy Smaller, Program Director 4 Case Management; Early Intervention; Integrate... # SA0631165 RECOVERY CONSULTANTS, INC. 248-543-1090 MI Treatment 48072 {'x': -9260261.7618, 'y': 5236694.745800003, '...
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 MI Treatment 48328 {'x': -9278492.0901, 'y': 5256952.089100003, '...
In [5]:
oakland = oakland_features.sdf
oakland.head()
Out[5]:
ACRES CNTY_CODE FID LABEL LAYOUT LINK NAME PENINSULA SQKM SQMILES Shape__Area Shape__Length TYPE VER SHAPE
0 4075.858559 125 1 BG 1 portrait 1251200001 1200001 lower 16.494592 6.368568 3.073481e+07 22611.540412 BG 14a {'rings': [[[-9250818.94177608, 5289778.006450...
1 8797.887541 125 2 BG 1 portrait 1251203001 1203001 lower 35.604171 13.746784 6.633129e+07 35417.658796 BG 14a {'rings': [[[-9264026.53064651, 5290079.339554...
2 6693.987915 125 3 BG 1 landscape 1251210001 1210001 lower 27.089900 10.459421 5.047787e+07 29332.713639 BG 14a {'rings': [[[-9270503.93750677, 5288614.193867...
3 5399.750207 125 4 BG 1 portrait 1251214001 1214001 lower 21.852249 8.437162 4.070574e+07 29772.396132 BG 14a {'rings': [[[-9269511.1796278, 5286236.4642479...
4 3719.286317 125 5 BG 2 portrait 1251222002 1222002 lower 15.051580 5.811421 2.804448e+07 23245.466116 BG 14a {'rings': [[[-9281344.03654134, 5290365.890503...

2. Drive time analysis

In [6]:
# Connect to GIS
from arcgis.gis import GIS
gis = GIS('http://dcdev.maps.arcgis.com/', 'username', 'password')

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
  2. [5] that is needed purely for visualization purposes

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

In [7]:
# Generate service areas
result = arcgis.features.use_proximity.create_drive_time_areas(clinic_url, [5, 10, 15, 20, 25], output_name="DriveTimeToClinics_All", overlap_policy='Dissolve')
In [8]:
# Check type of result
type(result)
Out[8]:
arcgis.gis.Item
In [ ]:
# Share results for public consumption
result.share(everyone=True)
In [10]:
# Convert it to a FeatureSet
result_url = str(result.url)+'/0/'
drive_times_all_layer = FeatureLayer(result_url)
driveAll_features = drive_times_all_layer.query()
print('Total number of rows in the service area (upto 25 mins) dataset: ')
print(len(driveAll_features.features))
Total number of rows in the service area (upto 25 mins) dataset: 
5
In [11]:
# Get results as a SpatialDataFrame
driveAll = driveAll_features.sdf
driveAll.head()
Out[11]:
AnalysisArea FacilityOID FromBreak Name OBJECTID ToBreak SHAPE
0 503.992431 None 20 20 - 25 1 25 {'rings': [[[-9247241.4852, 5209081.8363], [-9...
1 449.082348 None 15 15 - 20 2 20 {'rings': [[[-9262016.5505, 5298939.5558], [-9...
2 412.954912 None 10 10 - 15 3 15 {'rings': [[[-9244516.5052, 5228407.9967], [-9...
3 392.044252 None 5 5 - 10 4 10 {'rings': [[[-9301666.348, 5279365.8795], [-93...
4 332.783551 None 0 0 - 5 5 5 {'rings': [[[-9269666.4107, 5286863.1463], [-9...

Service area generation for the visible [5] minutes layer

In [12]:
# Generate service areas
visible_result = arcgis.features.use_proximity.create_drive_time_areas(clinic_url, [5], output_name="DriveTimeToClinics_5mins", overlap_policy='Dissolve')
In [13]:
# Grant 'visible_result' public access
visible_result.share(everyone=True)
Out[13]:
{'itemId': 'b4d0dd9ad9544990aab465d8063102ab', 'notSharedWith': []}
In [14]:
# Generate the URL for this layer (later used for visualization)
visible_result_url = str(visible_result.url)+'/0'
visible_result_url
Out[14]:
'https://services.arcgis.com/bkrWlSKcjUDFDtgw/arcgis/rest/services/DriveTimeToClinics_5mins/FeatureServer/0'

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

In [15]:
# Perform spatial join between CBG layer and the service areas created for all time durations
cbg_with_driveTime = arcgis.features.analysis.join_features(oakland_cbg, result_url, spatial_relationship='Intersects', output_name='Oakland_CBG_allDriveTimes')
In [16]:
# Grant 'cbg_with_driveTime' public access
cbg_with_driveTime.share(everyone=True)
Out[16]:
{'itemId': '7b3da4f6386a426ea83608032626d55f', 'notSharedWith': []}
In [17]:
# Converting it to a FeatureSet
cbg_with_driveTime_url = str(cbg_with_driveTime.url)+'/0/'
drive_times_cbg = FeatureLayer(cbg_with_driveTime_url)
drive_times_features = drive_times_cbg.query()
print('Total number of rows in the drive time (upto 25 mins) dataset: ')
print(len(drive_times_features.features))
Total number of rows in the drive time (upto 25 mins) dataset: 
934
In [18]:
# Convert it to a SpatialDataFrame
drive_time_cbg_df = drive_times_features.sdf
drive_time_cbg_df.head()
Out[18]:
ACRES AnalysisArea CNTY_CODE FacilityOID FromBreak Join_Count LABEL LAYOUT LINK NAME ... OBJECTID PENINSULA SQKM SQMILES Shape__Area Shape__Length TYPE ToBreak VER SHAPE
0 4075.858559 503.992431 125 None 20 2 BG 1 portrait 1251200001 1200001 ... 1 lower 16.494592 6.368568 3.073481e+07 22611.540412 BG 25 14a {'rings': [[[-9250818.9418, 5289778.0065], [-9...
1 8797.887541 503.992431 125 None 20 4 BG 1 portrait 1251203001 1203001 ... 2 lower 35.604171 13.746784 6.633129e+07 35417.658796 BG 25 14a {'rings': [[[-9264026.5306, 5290079.3396], [-9...
2 6693.987915 449.082348 125 None 15 3 BG 1 landscape 1251210001 1210001 ... 3 lower 27.089900 10.459421 5.047787e+07 29332.713639 BG 20 14a {'rings': [[[-9270503.9375, 5288614.1939], [-9...
3 5399.750207 412.954912 125 None 10 3 BG 1 portrait 1251214001 1214001 ... 4 lower 21.852249 8.437162 4.070574e+07 29772.396132 BG 15 14a {'rings': [[[-9269511.1796, 5286236.4643], [-9...
4 3719.286317 449.082348 125 None 15 3 BG 2 portrait 1251222002 1222002 ... 5 lower 15.051580 5.811421 2.804448e+07 23245.466116 BG 20 14a {'rings': [[[-9281344.0365, 5290365.8905], [-9...

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 [19]:
# Perform data enrichment
enriched_drive_time_cbg = arcgis.features.enrich_data.enrich_layer(drive_times_cbg, analysis_variables=['TOTPOP_CY', 'X8002_X', 'BACHDEG_CY'], output_name="Enriched_DriveTime_CBG")
In [20]:
# Check type of 'enriched_drive_time_cbg'
type(enriched_drive_time_cbg)
Out[20]:
arcgis.gis.Item
In [21]:
# Grant 'enriched_drive_time_cbg' public access
enriched_drive_time_cbg.share(everyone=True)
Out[21]:
{'itemId': 'dd47f3d9754f4559b641ca946d1e2c42', 'notSharedWith': []}
In [22]:
# Convert it to a Feature Set
enriched_url = str(enriched_drive_time_cbg.url)+'/0/'
enriched_driveTimes_layer = FeatureLayer(enriched_url)
enriched_driveTimes_features = enriched_driveTimes_layer.query()
print('Total number of rows in the drive time (upto 25 mins) dataset: ')
print(len(enriched_driveTimes_features.features))
Total number of rows in the drive time (upto 25 mins) dataset: 
934
In [23]:
# Convert it to a SpatialDataFrame
enriched_df = enriched_driveTimes_features.sdf
enriched_df.head()
Out[23]:
ACRES AnalysisArea BACHDEG_CY CNTY_CODE ENRICH_FID FacilityOID FromBreak HasData ID Join_Count ... Shape__Area Shape__Length TOTPOP_CY TYPE ToBreak VER X8002_X aggregationMethod sourceCountry SHAPE
0 4075.858559 503.992431 115 125 1 None 20 1 0 2 ... 3.073481e+07 22611.540412 856 BG 25 14a 1616323 BlockApportionment:US.BlockGroups US {'rings': [[[-9250818.9418, 5289778.0065], [-9...
1 8797.887541 503.992431 228 125 2 None 20 1 1 4 ... 6.633129e+07 35417.658796 1599 BG 25 14a 2979135 BlockApportionment:US.BlockGroups US {'rings': [[[-9264026.5306, 5290079.3396], [-9...
2 6693.987915 449.082348 66 125 3 None 15 1 2 3 ... 5.047787e+07 29332.713639 649 BG 20 14a 1284444 BlockApportionment:US.BlockGroups US {'rings': [[[-9270503.9375, 5288614.1939], [-9...
3 5399.750207 412.954912 858 125 4 None 10 1 3 3 ... 4.070574e+07 29772.396132 4185 BG 15 14a 8096459 BlockApportionment:US.BlockGroups US {'rings': [[[-9269511.1796, 5286236.4643], [-9...
4 3719.286317 449.082348 196 125 5 None 15 1 4 3 ... 2.804448e+07 23245.466116 1686 BG 20 14a 2423423 BlockApportionment:US.BlockGroups US {'rings': [[[-9281344.0365, 5290365.8905], [-9...

5 rows × 29 columns

In [24]:
# Print column names for ease of use during visualization
enriched_df.columns
Out[24]:
Index(['ACRES', 'AnalysisArea', 'BACHDEG_CY', 'CNTY_CODE', 'ENRICH_FID',
       'FacilityOID', 'FromBreak', 'HasData', 'ID', 'Join_Count', 'LABEL',
       'LAYOUT', 'LINK', 'NAME', 'Name_1', 'OBJECTID', 'PENINSULA', 'SQKM',
       'SQMILES', 'Shape__Area', 'Shape__Length', 'TOTPOP_CY', 'TYPE',
       'ToBreak', 'VER', 'X8002_X', 'aggregationMethod', 'sourceCountry',
       'SHAPE'],
      dtype='object')
In [25]:
# Format URL for visualization
enriched_lyr_url = enriched_url[:-1]
enriched_lyr_url
Out[25]:
'https://services.arcgis.com/bkrWlSKcjUDFDtgw/arcgis/rest/services/Enriched_DriveTime_CBG/FeatureServer/0'

5. Visualize results

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

pythonresultmap

In [27]:
# Add next layer for each CBG viualized based on its health insurance spending (X8002_2)
map1.add_layer({"type":"FeatureLayer",
                "url": enriched_lyr_url,
                "renderer":"ClassedSizeRenderer",
                "field_name":"X8002_X",
                "opacity":0.6
               })
In [28]:
# Add layer for service areas reachable within a 5 min drive time of the clinics
map1.add_layer({"type":"FeatureLayer",
                "url": visible_result_url,
                "opacity":0.75
               })
In [29]:
# Add layer of clinics
map1.add_layer({"type":"FeatureLayer",
                "url": clinic_url,
                "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)

Feedback on this topic?