<jupyter_start><jupyter_text>Finding suitable spots for placing heart defibrillator equipments in publicIn this sample, we will observe how site [suitability analyses](https://en.wikipedia.org/wiki/Suitability_analysis) 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)](https://en.wikipedia.org/wiki/Automated_external_defibrillator) for public emergencies.The criteria for a suitable place are those that have high incidence of [OHCA (Out of Hospital Cardiac Arrests)](http://www.sca-aware.org/sca-news/aha-releases-latest-statistics-on-out-of-hospital-cardiac-arrest) 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.  Connect to the GIS<jupyter_code>from arcgis.gis import *
gis = GIS("https://python.playground.esri.com/portal","arcgis_python","amazing_arcgis_123")<jupyter_output><empty_output><jupyter_text>Preview the input datasets<jupyter_code>ohca_item = gis.content.search("Pittsburgh_heart_attacks", "Feature Layer")[0]
ohca_item
map1 = gis.map("Pittsburgh, PA", zoomlevel=12)
map1
map1.add_layer(ohca_item)<jupyter_output><empty_output><jupyter_text>Let us take a look at the layers available in this item<jupyter_code>for lyr in ohca_item.layers:
    print(lyr.properties.name)<jupyter_output>Heart attack incidence
Streets
Zoning
Rivers
Pittsburgh boundary<jupyter_text>Outline of the analysisFor the rest of this analysis, we will use the `Zoning` polygon layer and build a `600`-foot buffers around it to represent areas that are accessible to commercial zones. Next, we will use the `Heart attack incidence` point layer to build a density raster. This raster depicts those areas that have a higher incidence of cardiac arrests. Finally, we will overlay the buffers on the density raster to pick places that are suitable to place new AED devices.  Create a 600 feet buffer around commercial areasThe `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 zonesTo 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.<jupyter_code>zoning_flayer = ohca_item.layers[2]
zoning_fset = zoning_flayer.query(result_record_count=10, return_all_records=False)
zoning_fset.df<jupyter_output><empty_output><jupyter_text>The column `zoning_grouped` contains zoning categories. We are intersted in those polygons that correspond to the `Commercial` category.<jupyter_code>zoning_commercial_fset = zoning_flayer.query("zoning_grouped = 'Commercial'")
commercial_zone_df = zoning_commercial_fset.df
commercial_zone_df.head(5)[['name','zoning_grouped']] #display the first 5 results
commercial_zone_df.shape<jupyter_output><empty_output><jupyter_text>Let us draw the selected polygons on a map<jupyter_code>zone_map = gis.map("Pittsburgh, PA")
zone_map
zone_map.draw(zoning_commercial_fset)<jupyter_output><empty_output><jupyter_text>Thus, from `965` zoning polygons, we have narrowed down to `317`. Create buffers on in-memory FeatureCollectionThe ArcGIS API for Python allows you to perform analysis on in-memory objects, meaning, you do not have to create feature layers for intermediate results, such as for the output of the previous feature selection. To execute the `create_buffers` tool, create an in-memory `FeatureCollection` object using these selected features.<jupyter_code>from arcgis.features import FeatureCollection
zoning_commercial_fc = FeatureCollection(zoning_commercial_fset.to_dict())
from arcgis.features import use_proximity<jupyter_output><empty_output><jupyter_text>Use the `zoning_commercial_fc` as an input for the buffer tool<jupyter_code>commercial_buffers = use_proximity.create_buffers(zoning_commercial_fc, 
                                                  distances=[600],units='Feet', 
                                                  dissolve_type='Dissolve',
                                                 output_name='commercial_buffers')<jupyter_output><empty_output><jupyter_text>Draw the results on the `commercial_zone_map` created above<jupyter_code>zone_map.add_layer(commercial_buffers)<jupyter_output><empty_output><jupyter_text>Create a density map to find areas of high heart attack incidenceTo 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.<jupyter_code>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.df.head(10)<jupyter_output><empty_output><jupyter_text>Calculate density<jupyter_code>from arcgis.raster.analytics import calculate_density
from arcgis.raster.functions import *
ha_density = calculate_density(ha_incidence, count_field='num_incidence', 
                               output_cell_size={'distance':150,'units':'feet'},
                               output_name = 'ha_density')
ha_density<jupyter_output><empty_output><jupyter_text>Let us display the density raster on a map<jupyter_code>density_map = gis.map("Pittsburgh, PA", zoomlevel=11)
density_map<jupyter_output><empty_output><jupyter_text>Use the `stretch` [raster function](http://pro.arcgis.com/en/pro-app/help/data/imagery/stretch-function.htm) to enhance the density layer before adding it to the map:<jupyter_code>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')
density_map.add_layer(colormap_rf, {"opacity":0.5})<jupyter_output><empty_output><jupyter_text>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 rasterCalculate 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<jupyter_code>import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>density_hist = density_layer.histograms<jupyter_output><empty_output><jupyter_text>Construct the X axis such that it ranges from min value to max value of the pixel range in the image.<jupyter_code>x = np.linspace(density_hist[0]['min'], density_hist[0]['max'], num=density_hist[0]['size'])
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")<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>conversion_value = (1500*1500)/(5280*5280)
density_5blocks = density_layer * conversion_value #raster arithmetic
density_5blocks<jupyter_output><empty_output><jupyter_text>Let us remap this continuous density raster to a binary layer representing whether a pixel represents high enough density or not.<jupyter_code>density_classified_color = colormap(density_5blocks, colormap_name='Random',astype='u8')
density_classified_color<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>#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<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>density_map2 = gis.map("Pittsburgh, PA")
density_map2
density_map2.add_layer(density_classified_viz)<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>from arcgis.raster import convert_feature_to_raster
buffer_raster = convert_feature_to_raster(commercial_buffers.layers[0],
                                          output_cell_size={'distance':150, 'units':'feet'},
                                          output_name='buffer_raster')<jupyter_output><empty_output><jupyter_text>Query the layer to quickly visualize it as an image<jupyter_code>buffer_raster<jupyter_output><empty_output><jupyter_text>The `raster` module of the Python API provides numerous raster functions. Of which we use the [bitwise_and](http://desktop.arcgis.com/en/arcmap/10.3/manage-data/raster-and-images/local-function.htm) local function which returns an image with pixels that match in both the input rasters.<jupyter_code>bool_overlay = bitwise_and([buffer_raster,density_classified])
bool_overlay<jupyter_output><empty_output><jupyter_text>Let us overlay this final result on a map to visualize the regions that are suitable to locating new AED devices.<jupyter_code>map3 = gis.map("Carnegie Mellon University, PA")
map3<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.07572009188902633, 'fraction_numerical': 0.012546386287329916, 'mean_word_length': 5.944171779141104, 'pattern_counts': {'":': 1, '<': 96, '<?xml version=': 0, '>': 96, 'https://': 3, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '24', 'dir': 'jupyter-structured-clean-dedup', 'id': '342315', 'license': 'Apache-2.0', 'n_tokens_mistral': 3426, 'n_tokens_neox': 3249, 'n_words': 1102, 'path': 'samples/04_gis_analysts_data_scientists/finding_suitable_spots_for_AED_devices_using_raster_analytics.ipynb', 'repo_name': 'merinogeospatial/arcgis-python-api', 'size': '11318'}