Find hospitals closest to an incident

The network module of the ArcGIS API for Python can be used to solve different types of network analysis operations. In this sample, we see how to find the hospital that is closest to an incident.

Closest facility

The closest facility solver provides functionality for finding out the closest locations to a particular input point. This solver would be useful in cases when you have an incident and need to find the closest facility or need to get information on the travel time and the distance to each of the facilities from an incident point for reporting purposes.

When finding closest facilities, you can specify how many to find and whether the direction of travel is toward or away from them. The closest facility solver displays the best routes between incidents and facilities, reports their travel costs, and returns driving directions.

Connect to your GIS

As a first step, you would need to establish a connection to your organization which could be an ArcGIS Online organization or an ArcGIS Enterprise.

from IPython.display import HTML
import pandas as pd
from arcgis.gis import GIS

#connect to your GIS
my_gis = GIS("home")

Create a Network Layer

To perform any network analysis (such as finding the closest facility, the best route between multiple stops, or service area around a facility), you would need to create a NetworkLayer object. In this sample, since we are solving for closest facilities, we need to create a ClosestFacilityLayer which is a type of NetworkLayer.

To create any NetworkLayer object, you would need to provide the URL to the appropriate network analysis service. Hence, in this sample, we provide a ClosestFacility URL to create a ClosestFacilityLayer object.

Since all ArcGIS Online organizations already have access to those routing services, you can access this URL through the GIS object's helperServices property. If you have your own ArcGIS Server based map service with network analysis capability enabled, you would need to provide the URL for this service.

Let us start by importing the network module

import arcgis.network as network

Access the analysis URL from the GIS object

analysis_url = my_gis.properties.helperServices.closestFacility.url
analysis_url

'https://route.arcgis.com/arcgis/rest/services/World/ClosestFacility/NAServer/ClosestFacility_World'

Create a ClosestFacilityLayer object using this URL

cf_layer = network.ClosestFacilityLayer(analysis_url, gis=my_gis)

Create hospitals layer

In this sample, we will be looking for the closest hospital (facility) to an incident location. Even though we are interested in finding out the closest one, it would still be helpful to get the information on the distance and travel time to all of them for reference purposes.

In the code below, we need to geocode the hospitals' addresses as well as do the reverse geocode for the incident location which has been supplied in the latitude/longitude format.

To perform the geocode operations, we import the geocoding module of the ArcGIS API.

from arcgis import geocoding

In this sample, we geocode addresses of hospitals to create the facility layer. In your workflows, this could any feature layer. Create a list of hospitals in Rio de Janeiro, Brazil.

hospitals_addresses = ['Estrada Adhemar Bebiano, 339 Del Castilho, Rio de Janeiro RJ, 21051-370, Brazil',
                       'R. José dos Reis Engenho de Dentro, Rio de Janeiro RJ, 20750-000, Brazil',
                       'R. Dezessete, s/n Maré, Rio de Janeiro RJ, 21042-010, Brazil',
                       'Rua Dr. Miguel Vieira Ferreira, 266 Ramos, Rio de Janeiro RJ, Brazil']

Loop through each address and geocode it. The geocode operation returns a list of matches for each address. We pick the first result and extract the coordinates from it and construct a Feature object out of it. Then we combine all the Features representing the hospitals into a FeatureSet object.

from arcgis.features import Feature, FeatureSet

hosp_feat_list = []

for address in hospitals_addresses:
    hit = geocoding.geocode(address)[0]
    hosp_feat = Feature(geometry=hit['location'], attributes=hit['attributes'])

    hosp_feat_list.append(hosp_feat)

Construct a FeatureSet using each hospital Feature.

hospitals_fset = FeatureSet(features=hosp_feat_list, 
                            geometry_type='esriGeometryPoint', 
                            spatial_reference={'wkid' : 4326, 'latestWkid': 4326})

Lets draw our hospitals on a map

from IPython.display import HTML

map1 = my_gis.map('Rio de Janeiro, Brazil')
map1.basemap = 'arcgis-light-gray'
map1

map1.draw(hospitals_fset, symbol={"type": "esriSMS","style": "esriSMSSquare",
                                  "color": [76,115,0,255],"size": 8,})

Create incidents layer

Similarly, let us create the incient layer

incident_coords = '-43.281206,-22.865676'
reverse_geocode = geocoding.reverse_geocode({"x": incident_coords.split(',')[0], 
                                              "y": incident_coords.split(',')[1]})

incident_feature = Feature(geometry=reverse_geocode['location'], 
                           attributes=reverse_geocode['address'])

incident_fset = FeatureSet([incident_feature], geometry_type='esriGeometryPoint',
                          spatial_reference={'wkid' : 4326, 'latestWkid': 4326})

Let us add the incident to the map

map1.draw(incident_fset, symbol={"type": "esriSMS","style": "esriSMSCircle","size": 8})

Solve for closest hospital

By default the closest facility service would return only the closest location, so we need to specify explicitly the default_target_facility_count parameter as well as return_facilities.

result = cf_layer.solve_closest_facility(incidents=incident_fset,
                                        facilities=hospitals_fset,
                                        default_target_facility_count=4,
                                        return_facilities=True,
                                        impedance_attribute_name='TravelTime',
                                        accumulate_attribute_names=['Kilometers','TravelTime'])

Let us inspect the result dictionary

result.keys()

dict_keys(['requestID', 'routes', 'facilities', 'messages'])

Let us use the routes dictionary to construct line features out of the routes to display on the map

result['routes'].keys()

dict_keys(['fieldAliases', 'geometryType', 'spatialReference', 'fields', 'features'])

result['routes']['features'][0].keys()

dict_keys(['attributes', 'geometry'])

Construct line features out of the routes that are returned.

line_feat_list = []
for line_dict in result['routes']['features']:
    f1 = Feature(line_dict['geometry'], line_dict['attributes'])
    line_feat_list.append(f1)

routes_fset = FeatureSet(line_feat_list, 
                         geometry_type=result['routes']['geometryType'],
                         spatial_reference= result['routes']['spatialReference'])

Add the routes back to the map. The route to the closest hospital is in red

allocation_line_symbol_4 = {'type': 'esriSLS', 'style': 'esriSLSSolid',
                                'color': [0,0,255,100], 'width': 6}
map1.draw(routes_fset, symbol = allocation_line_symbol_4)

map1.clear_graphics()

Analyze the results in a table

Since we parsed the routes as a FeatureSet, we can display the attributes easily as a pandas DataFrame.

df1 = routes_fset.sdf
df1

Let us add the hospital addresses and incident address to this table and display only the relevant columns

df1['facility_address'] = hospitals_addresses
df1['incident_address'] = [incident_feature.attributes['Match_addr'] for i in range(len(hospitals_addresses))]

df1[['facility_address','incident_address','Total_Miles','Total_TravelTime']]

Conclusion

Thus using the network module of the ArcGIS API for Python, you can solve for closest facilities from an incident location.