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.

In [1]:
from IPython.display import HTML
import pandas as pd
from arcgis.gis import GIS

#connect to your GIS
user_name = 'arcgis_python'
password = 'P@ssword123'
my_gis = GIS('https://www.arcgis.com', user_name, password)

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

In [2]:
import arcgis.network as network

Access the analysis URL from the GIS object

In [3]:
analysis_url = my_gis.properties.helperServices.closestFacility.url
analysis_url
Out[3]:
'https://route.arcgis.com/arcgis/rest/services/World/ClosestFacility/NAServer/ClosestFacility_World'

Create a ClosestFacilityLayer object using this URL

In [4]:
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.

In [5]:
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.

In [6]:
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.

In [7]:
from arcgis.features import Feature, FeatureSet
In [8]:
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.

In [9]:
hospitals_fset = FeatureSet(features=hosp_feat_list, 
                            geometry_type='esriGeometryPoint', 
                            spatial_reference={'latestWkid': 4326})

Lets draw our hospitals on a map

In [10]:
map1 = my_gis.map('Rio de Janeiro, Brazil')
map1

closest hosptitals route

In [11]:
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

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

Let us add the incident to the map

In [14]:
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.

In [15]:
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

In [16]:
result.keys()
Out[16]:
dict_keys(['messages', 'facilities', 'routes'])

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

In [76]:
result['routes'].keys()
Out[76]:
dict_keys(['fieldAliases', 'geometryType', 'spatialReference', 'features'])
In [80]:
result['routes']['features'][0].keys()
Out[80]:
dict_keys(['attributes', 'geometry'])

Construct line features out of the routes that are returned.

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

In [19]:
map1.draw(routes_fset)

Analyze the results in a table

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

In [20]:
routes_fset.df
Out[20]:
FacilityCurbApproach FacilityID FacilityRank IncidentCurbApproach IncidentID Name ObjectID Shape_Length Total_Kilometers Total_Miles Total_TravelTime SHAPE
0 1 1 1 2 1 Rua Tangapeme 2-50 - Location 1 1 0.028445 3.012284 1.871747 7.158442 {'paths': [[[-43.28124189499994, -22.865707432...
1 2 4 2 2 1 Rua Tangapeme 2-50 - Location 4 2 0.048777 5.218743 3.242777 11.946830 {'paths': [[[-43.28124189499994, -22.865707432...
2 1 3 3 2 1 Rua Tangapeme 2-50 - Location 3 3 0.066788 7.072069 4.394380 14.159053 {'paths': [[[-43.28124189499994, -22.865707432...
3 2 2 4 2 1 Rua Tangapeme 2-50 - Location 2 4 0.070686 7.467527 4.640106 14.285109 {'paths': [[[-43.28124189499994, -22.865707432...

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

In [23]:
df1 = routes_fset.df
df1['facility_address'] = hospitals_addresses
df1['incident_address'] = [incident_feature.attributes['Match_addr'] for i in range(len(hospitals_addresses))]
In [31]:
df1[['facility_address','incident_address','Total_Miles','Total_TravelTime']]
Out[31]:
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.871747 7.158442
1 R. José dos Reis Engenho de Dentro, Rio de Jan... Rua Tangapeme 2-50, Complexo do Alemão, Rio de... 3.242777 11.946830
2 R. Dezessete, s/n Maré, Rio de Janeiro RJ, 210... Rua Tangapeme 2-50, Complexo do Alemão, Rio de... 4.394380 14.159053
3 Rua Dr. Miguel Vieira Ferreira, 266 Ramos, Rio... Rua Tangapeme 2-50, Complexo do Alemão, Rio de... 4.640106 14.285109

Conclusion

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


Feedback on this topic?