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.

Input
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

Input
import arcgis.network as network

Access the analysis URL from the GIS object

Input
analysis_url = my_gis.properties.helperServices.closestFacility.url
analysis_url
Output
'https://route.arcgis.com/arcgis/rest/services/World/ClosestFacility/NAServer/ClosestFacility_World'

Create a ClosestFacilityLayer object using this URL

Input
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.

Input
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.

Input
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.

Input
from arcgis.features import Feature, FeatureSet
Input
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.

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

Lets draw our hospitals on a map

Input
from IPython.display import HTML
Input
map1 = my_gis.map('Rio de Janeiro, Brazil')
map1
Output
Input
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

Input
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'])
Input
incident_fset = FeatureSet([incident_feature], geometry_type='esriGeometryPoint',
                          spatial_reference={'latestWkid': 4326})

Let us add the incident to the map

Input
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.

Input
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

Input
result.keys()
Output
dict_keys(['messages', 'routes', 'facilities'])

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

Input
result['routes'].keys()
Output
dict_keys(['fieldAliases', 'geometryType', 'spatialReference', 'features'])
Input
result['routes']['features'][0].keys()
Output
dict_keys(['attributes', 'geometry'])

Construct line features out of the routes that are returned.

Input
line_feat_list = []
for line_dict in result['routes']['features']:
    f1 = Feature(line_dict['geometry'], line_dict['attributes'])
    line_feat_list.append(f1)
Input
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

Input
allocation_line_symbol_4 = {'type': 'esriSLS', 'style': 'esriSLSSolid',
                                'color': [0,0,255,100], 'width': 6}
map1.draw(routes_fset, symbol = allocation_line_symbol_4)
Input
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.

Input
df1 = routes_fset.sdf
df1
Output
FacilityCurbApproach FacilityID FacilityRank IncidentCurbApproach IncidentID Name ObjectID SHAPE Shape_Length Total_Kilometers Total_Miles Total_TravelTime
0 1 1 1 2 1 Rua Tangapeme 36 - Location 1 1 {"paths": [[[-43.28124189499994, -22.865707432... 0.028461 3.016318 1.874253 7.554929
1 2 4 2 2 1 Rua Tangapeme 36 - Location 4 2 {"paths": [[[-43.28124189499994, -22.865707432... 0.048805 5.222468 3.245091 12.202391
2 2 3 3 2 1 Rua Tangapeme 36 - Location 3 3 {"paths": [[[-43.28124189499994, -22.865707432... 0.061555 6.464249 4.016698 14.337250
3 2 2 4 2 1 Rua Tangapeme 36 - Location 2 4 {"paths": [[[-43.28124189499994, -22.865707432... 0.086735 9.143542 5.681533 18.793156

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

Input
df1['facility_address'] = hospitals_addresses
df1['incident_address'] = [incident_feature.attributes['Match_addr'] for i in range(len(hospitals_addresses))]
Input
df1[['facility_address','incident_address','Total_Miles','Total_TravelTime']]
Output
facility_address incident_address Total_Miles Total_TravelTime
0 Estrada Adhemar Bebiano, 339 Del Castilho, Rio... Rua Tangapeme 2-50, Complexo do Alemão, Rio de... 1.874253 7.554929
1 R. José dos Reis Engenho de Dentro, Rio de Jan... Rua Tangapeme 2-50, Complexo do Alemão, Rio de... 3.245091 12.202391
2 R. Dezessete, s/n Maré, Rio de Janeiro RJ, 210... Rua Tangapeme 2-50, Complexo do Alemão, Rio de... 4.016698 14.337250
3 Rua Dr. Miguel Vieira Ferreira, 266 Ramos, Rio... Rua Tangapeme 2-50, Complexo do Alemão, Rio de... 5.681533 18.793156

Conclusion

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

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