Finding suitable spots for placing heart defibrillator equipments in public

In this sample, we will observe how site suitability analyses can be performed using the ArcGIS API for Python. The objective of this sample is to find locations in the city of Philadelphia that are suitable for placing AED (Automated External Defibrillator) for public emergencies.

Image of an AED device attached to a wall at San Diego Convention Center during Esri UC

The criteria for suitable places are those that have high incidence of OHCA (Out of Hospital Cardiac Arrests) and be accessible to public, such as commercial areas.

As inputs, we start with geocoded OCHA (Out-of-Hospital Cardiac Arrest) point data, along with a few base layers for the city of Pittsburgh published as feature layers. As output, we need to generate a list of locations that have a high incidence of heart-attacks and located within commercial areas, allowing easy access at times of emergencies.

Getting set up

Input
from arcgis.gis import GIS
from arcgis.mapping import WebMap
from arcgis.widgets import MapView
from arcgis.features import FeatureCollection, use_proximity
from datetime import datetime
Input
gis = GIS(url='https://pythonapi.playground.esri.com/portal', username='arcgis_python', password='amazing_arcgis_123')

Load input datasets

Input
ohca_item = gis.content.get('a5719916dff4442789a59680c25a4284')
ohca_item
Output
Pittsburgh heart attacks
OHCA in PittsburghFeature Layer Collection by api_data_owner
Last Modified: June 21, 2018
0 comments, 102 views
Input
ohca_map_item = gis.content.get('b8b6cf2bcbeb4903a5372b7f4cbfb252')
ohca_map_item
Output
Pittsburgh heart attacks
Map showing cardiac arrest information in PittsburghWeb Map by api_data_owner
Last Modified: March 11, 2020
0 comments, 44 views

Let us take a look at the layers available in this item

Input
for lyr in ohca_item.layers:
    print(lyr.properties.name)
Heart attack incidence
Streets
Zoning
Boundary

Let us display the Web Map item to view these layers on a map.

Input
map1 = MapView(item=ohca_map_item)
map1.legend=True
map1

Outline of the analysis

The idea of this analysis is to find places suitable for placing the AED devices. Based on prior knowledge we happen to know areas that are commercial, accessible to public and showing a high incidence of out-of-hospital cardiac arrests are good candidates. We will build the suitability model by performing these steps:

  • use Zoning polygon layer to identify commercial areas and build a 600 feet buffer around them
  • perform density analysis on Heart attack incidence point layer
  • perform spatial overlay to find overlapping locations

Create a 600 feet buffer around commercial areas

The Zoning feature layer contains polygon features that represent different zones such as commercial, residential etc. We need to select those features that correspond to commercial zones and create a buffer of 600 feet around them. The 600 feet area roughly corresponds to two-blocks, a walk able distance in case of an emergency.

Select commercial zones

To select the commercial zones using a query, we need to know what columns and values are available. Hence, let us construct a small query that gives the first few rows / features.

Input
zoning_flayer = ohca_item.layers[2]
zoning_sdf = zoning_flayer.query(result_record_count=5, return_all_records=False, as_df=True)
Input
zoning_sdf
Output
objectid area perimeter zoning_ zoning_id zon_new symbol acres sqmiles code name x y code_2 zoning_grouped SHAPE
0 1 1.591398e+07 52563.176 2 5521 R1D-L 80 370.108 0.576 R1D-L Single-Unit Detached Residential/Low Density 1.340082e+06 429618.202356 R1D-L Residential {'rings': [[[-80.01772016299998, 40.4977767210...
1 2 8.107219e+04 1524.910 3 5522 R2-L 80 370.108 0.576 R2-L Two-Unit Residential/Low Density 1.338226e+06 432380.904671 R2-L Non Commercial {'rings': [[[-80.01757926499994, 40.4981592720...
2 3 2.591469e+06 14071.084 4 5526 R2-L 83 60.461 0.094 R2-L Two-Unit Residential/Low Density 1.335954e+06 430498.618340 R2-L Non Commercial {'rings': [[[-80.02201914599993, 40.4938213470...
3 4 8.370962e+05 5395.625 5 5517 RP 85 19.306 0.030 RP Residential Planned Unit Development 1.340037e+06 431409.546743 RP Residential {'rings': [[[-80.01117829099985, 40.4967724760...
4 5 9.650762e+03 441.323 6 5519 RM-M 84 0.223 0.000 RM-M Multi-Unit Residential/Low Density 1.339642e+06 431798.332613 RM-M Residential {'rings': [[[-80.01130567899997, 40.4967769370...

The column zoning_grouped contains zoning categories. We are intersted in those polygons that correspond to the Commercial category.

Input
zoning_commercial_fset = zoning_flayer.query(where="zoning_grouped = 'Commercial'")
commercial_zone_df = zoning_commercial_fset.sdf
commercial_zone_df.head(5)[['name','zoning_grouped']] #display the first 5 results
Output
name zoning_grouped
0 Urban Industrial District Commercial
1 Local Neighborhood Commercial Commercial
2 Local Neighborhood Commercial Commercial
3 General Industrial Commercial
4 Educational/Medical Institutional Commercial
Input
commercial_zone_df.shape
Output
(317, 16)

Let us draw the selected polygons on a map

Input
zone_map = gis.map("Pittsburgh, PA")
zone_map
Input
zone_map.draw(zoning_commercial_fset)

Thus, from 965 zoning polygons, we have narrowed down to 317.

Create buffers around commercial zones

The ArcGIS API for Python allows you to define definition queries or filters on Feature Layers. When you run a spatial analysis on those layers, only the features that fit the filter criteria you specified will be used. Thus, you can use the 'where' clause you used earlier (to get commercial zones) to set as a filter on the zoning_flayer and pass that as the input to the create_buffers tool. The advantage of this workflow is, you are not sending the features from the local FeatureSet object to the tool, instead, you are asking to the tool to get the features directly from the feature layer which is colocated with the tool. This paradigm of colocating the compute with the data is highly preferred to improve efficiency and scalability of your analyses.

Input
# create a filter using the where clause from earlier
zoning_flayer.filter = "zoning_grouped = 'Commercial'"
Input
# create a timestamp to create a unique output
timestamp=datetime.now().strftime('%d_%m_%Y_%H_%M_%S')

# create buffers
commercial_buffers = use_proximity.create_buffers(input_layer=zoning_flayer,
                                                  distances=[600],units='Feet', 
                                                  dissolve_type='Dissolve',
                                                 output_name=f'commercial_buffers_{timestamp}')
commercial_buffers
Output
commercial_buffers_28_05_2021_09_36_44
Feature Layer Collection by arcgis_python
Last Modified: May 28, 2021
0 comments, 0 views

Draw the results on the commercial_zone_map created above

Input
zone_map.add_layer(commercial_buffers)

Create a density map to find areas of high heart attack incidence

To calculate the density, we use calculate_density tool available under the raster module and provide the Heart attack incidence feature layer as its input. This layer has a column named num_incidence that additionally specifies the number of heart attacks that happened at each point location. Below we bring up a few of the features to get an idea.

Input
ha_incidence = ohca_item.layers[0] #the first layer in the input feature layer collection
ha_incidence_fset = ha_incidence.query(result_record_count=10, return_all_records=False)
ha_incidence_fset.sdf.head(10)
Output
objectid_1 fid_1 id pop2000 no yes num_incidence SHAPE
0 1 0 1 0 1 0 0 {"x": -79.97274830899988, "y": 40.437756305000...
1 2 8 9 96 1 0 0 {"x": -79.97639852099996, "y": 40.437202953000...
2 3 12 13 5 1 0 0 {"x": -79.98023401899997, "y": 40.438334899000...
3 4 13 14 3 1 0 0 {"x": -79.9818761219999, "y": 40.4383995900001...
4 5 25 26 15 1 0 0 {"x": -79.98428402499985, "y": 40.437456611000...
5 6 28 29 8 1 0 0 {"x": -79.98319929899998, "y": 40.436778390000...
6 7 30 31 42 1 0 0 {"x": -79.98183133499987, "y": 40.437055132000...
7 8 51 52 0 4 0 0 {"x": -79.98873209699991, "y": 40.435164507000...
8 9 70 71 0 2 0 0 {"x": -79.98717537399989, "y": 40.437418760000...
9 10 72 73 536 1 0 0 {"x": -79.9918783149999, "y": 40.4381269750001...

Calculate density

Input
from arcgis.raster.analytics import calculate_density
from arcgis.raster.functions import *
Input
# create a timestamp to create a unique output
timestamp=datetime.now().strftime('%d_%m_%Y_%H_%M_%S')

ha_density = calculate_density(ha_incidence, count_field='num_incidence', 
                               output_cell_size={'distance':150,'units':'feet'},
                               output_name = f'ha_density_{timestamp}')
print(ha_density)
<Item title:"ha_density_28_05_2021_09_38_24" type:Imagery Layer owner:arcgis_python>

Let us display the density raster on a map

Input
density_map = gis.map("Pittsburgh, PA", zoomlevel=11)
density_map

Use the stretch raster function to enhance the density layer before adding it to the map:

Input
density_layer = ha_density.layers[0]

stretch_rf = stretch(density_layer, stretch_type='StdDev',num_stddev=2)
colormap_rf = colormap(stretch_rf, colormap_name='Gray')
Input
density_map.add_layer(colormap_rf, {"opacity":0.5})

From the density_map, we see certain regions (in shades of white) have a higher density of heart attack incidences compared to the rest.

Reclassify the density raster

Calculate density tool returns the number of incidences per sq.mile. We are interested in the number of heart attacks at a larger scale of about 5 square blocks. In Pittsburgh, each block spans about 300 ft in length, thus 5 sq. blocks cover an area of 1500 x 1500 sq.feet. We apply remap raster function to convert the density from sq. miles to that in 5 block area

Input
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np

Plot the histogram to view actual density values and its distribution. The histograms property of the ImageryLayer object returns you histogram of each of its bands.

Input
density_hist = density_layer.histograms

Construct the X axis such that it ranges from min value to max value of the pixel range in the image.

Input
x = np.linspace(density_hist[0]['min'], density_hist[0]['max'], num=density_hist[0]['size'])
Input
fig = plt.figure(figsize=(10,8))
ax = fig.add_axes([0,0,1,1])
ax.bar(x,density_hist[0]['counts'])
ax.set_title("Histogram of heart attack density layer")
ax.set_xlabel("Heart attacks per sq. mile")
ax.set_ylabel("Number of pixels")

ax2 = fig.add_axes([0.25,0.2,0.7,0.7])
ax2.bar(x[-200:], density_hist[0]['counts'][-200:])
ax2.set_title("Histogram of heart attack density layer - zoomed")
ax2.set_xlabel("Heart attacks per sq. mile")
ax2.set_ylabel("Number of pixels")
Output
Text(0, 0.5, 'Number of pixels')

Convert units from sqmile to city blocks

The inset histogram chart has the histogram zoomed to view the distribution in the upper end of the density spectrum. We are interested in selecting those regions that have a heart attack of at least 5 per 5 block area. To achieve this, we need to convert the density from square miles to 5 square blocks.

Input
conversion_value = (1500*1500)/(5280*5280)
density_5blocks = density_layer * conversion_value #raster arithmetic
density_5blocks
Output

Let us remap this continuous density raster to a binary layer representing whether a pixel represents high enough density or not.

Input
density_classified_color = colormap(density_5blocks, colormap_name='Random',astype='u8')
density_classified_color
Output

Next, we classify the density raster such that pixels that have heart attacks greater than 5 get value 1 and rest become 'no data' pixels.

Input
#remap pixel values to create a binary raster
density_classified = remap(density_5blocks, input_ranges=[5,16], output_values=[1],astype='u8',no_data_ranges=[0,5])
density_classified_viz = colormap(density_classified, colormap_name='Random', astype='u8')
density_classified_viz
Output

Through classification, we have determined there are 3 hotspots in our density raster. Let us overlay this on a map to see which areas these hotspots correspond to.

Input
density_map2 = gis.map("Pittsburgh, PA")
density_map2
Input
density_map2.add_layer(density_classified_viz)

Perform overlay analysis

The site selection condition requires two inputs, the heart attack density layer (which we created earlier) and the accessibility layer (from the buffer analysis). To perform overlay, we need to convert the buffers layer to a raster layer of matching cell size as that of the density raster layer. To perform this conversion we use the convert_feature_to_raster method.

Input
from arcgis.raster.analytics import convert_feature_to_raster
Input
# create a timestamp to create a unique output
timestamp=datetime.now().strftime('%d_%m_%Y_%H_%M_%S')

# convert zoning buffer polygon to a raster layer of matching cell size
buffer_raster = convert_feature_to_raster(commercial_buffers.layers[0],
                                          output_cell_size={'distance':150, 'units':'feet'},
                                          output_name=f'buffer_raster_{timestamp}')

print(buffer_raster)
<Item title:"buffer_raster_28_05_2021_09_44_37" type:Imagery Layer owner:arcgis_python>

Query the layer to quickly visualize it as an image

Input
buffer_raster = buffer_raster.layers[0]
buffer_raster
Output

The raster module of the Python API provides numerous raster functions. Of which we use the bitwise_and local function which returns an image with pixels that match in both the input rasters.

Input
bool_overlay = bitwise_and([buffer_raster,density_classified])
bool_overlay
Output

Let us overlay this final result on a map to visualize the regions that are suitable to locating new AED devices.

Input
map3 = gis.map("Carnegie Mellon University, PA")
map3
Output
Input
map3.add_layer(bool_overlay)
 

Conclusion

Thus, in this sample, we observed how site-suitability analyses can be performed using ArcGIS and the ArcGIS API for Python. We started with the requirements for placing new AED devices as -- high intensity of cardiac arrests and proximity to commercial areas. Using a combination of feature analysis and raster analysis, we were able to process and extract the suitable sites. The analyst could convert the results from raster to vector, perform a centroid operation on the polygons, followed by reverse geocode to get the addresses of these 3 suitable locations for reporting and further action.

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