Drive time analysis for opioid epidemic

#Import necessary packages
from arcgis.gis import GIS
from datetime import datetime
from arcgis.features import FeatureLayer

from arcgis.features.use_proximity import create_drive_time_areas
from arcgis.features.analysis import join_features
from arcgis.features.enrich_data import enrich_layer

# Connect to GIS
gis = GIS(profile="your_online_profile") 

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

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

m1.basemap.basemap = 'topo-vector'

m1.content.add(oakland_layer)
m1.content.add(clinic_layer)

m1.zoom_to_layer(clinic_layer)

m1.legend.enabled = True

# Generate service areas
result = create_drive_time_areas(input_layer = clinic_layer, 
                                 travel_mode = 'Driving Time',
                                 break_values = [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.sharing.sharing_level = 'EVERYONE'
result.sharing.shared_with

{'groups': [], 'level': <SharingLevel.EVERYONE: 'EVERYONE'>}

# 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.content.add(result.layers[0])

m_drive_times.zoom_to_layer(result.layers[0])

m_drive_times.legend.enabled = True

# Generate service areas
visible_result = create_drive_time_areas(input_layer = clinic_layer, 
                                         break_values = [5], 
                                         travel_mode = 'Driving Time',
                                         output_name = "DriveTimeToClinics_5mins"+str(datetime.now().microsecond), 
                                         overlap_policy = 'Dissolve')

# Grant 'visible_result' public access
visible_result.sharing.sharing_level = "EVERYONE"
visible_result.sharing.shared_with

{'groups': [], 'level': <SharingLevel.EVERYONE: 'EVERYONE'>}

# Perform spatial join between CBG layer and the service areas created for all time durations
cbg_with_driveTime = join_features(target_layer = oakland_cbg, 
                                   join_layer = 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.sharing.sharing_level = "EVERYONE"
cbg_with_driveTime.sharing.shared_with

{'groups': [], 'level': <SharingLevel.EVERYONE: 'EVERYONE'>}

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

# 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.sharing.sharing_level = "EVERYONE"
enriched_drive_time_cbg.sharing.shared_with

{'groups': [], 'level': <SharingLevel.EVERYONE: 'EVERYONE'>}

# 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(['OBJECTID', 'Join_Count', 'TARGET_FID', 'NAME', 'LABEL', 'TYPE', 'LINK',
       'CNTY_CODE', 'SQKM', 'SQMILES', 'ACRES', 'VER', 'LAYOUT', 'PENINSULA',
       'Name_1', 'FromBreak', 'ToBreak', 'FacilityOID', 'AnalysisArea', 'ID',
       'sourceCountry', 'ENRICH_FID', 'aggregationMethod',
       'populationToPolygonSizeRating', 'apportionmentConfidence', 'HasData',
       'TOTPOP_CY', 'X8002_X', 'BACHDEG_CY', 'SHAPE'],
      dtype='object')

# Create map of Oakland county
map1 = gis.map('Oakland, Michigan')
map1.basemap.basemap = 'arcgis-dark-gray'
map1

map1.content.add(enriched_driveTimes_layer)
#Use smart mapping to render size based classes
sm_manager = map1.content.renderer(0).smart_mapping()
sm_manager.class_breaks_renderer(break_type="size", field="X8002_X", num_classes=5)

#Visualize drive time areas
map1.content.add(visible_result.layers[0], options = {'opacity':0.6})

map1.content.add(clinic_layer, options = {'opacity':0.85})

map1.zoom_to_layer(enriched_driveTimes_layer)