Drive time analysis for opioid epidemic

import arcgis
from datetime import datetime
from arcgis.features import FeatureLayer

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

# 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

# Convert to SpatialDataFrames
clinic = clinic_features.sdf
clinic.head()

oakland = oakland_features.sdf
oakland.head()

m = gis.map('Oakland, Michigan')
m

m.add_layer(oakland_layer)
m.add_layer(clinic_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

# Get results as a SpatialDataFrame
driveAll = driveAll_features.sdf
driveAll.head()

m_drive_times = gis.map('Oakland, Michigan')
m_drive_times

m_drive_times.add_layer(result.layers[0])

# 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': []}]}

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

# Convert it to a SpatialDataFrame
drive_time_cbg_df = drive_times_features.sdf
drive_time_cbg_df.head()

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

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

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