2017

July 2017

# Raster Functions and the ArcGIS API for Python 1.2

Posted by david_johnson-esristaff Jul 11, 2017

One of the great things about working in the Lab is you get to experiment with the new goodies from our core software developers before they are released.  When I heard that version 1.2 of the ArcGIS API for Python would include a new module for raster functions, I could not wait to give it a try.  Now that v.1.2 of the API is released, I can finally show you a Jupyter Notebook I built which has an example of a weighted overlay analysis implemented with raster functions.   The following is a non-interactive version of that notebook which I exported to HTML.  I hope it will give you some ideas for how you could use the ArcGIS API for Python to perform your own raster analysis.

Finding Natural and Accessible Areas in the State of Washington, USA

The weighted overlay is a standard GIS analysis technique for site-suitability and travel cost studies. This notebook leverages the new "arcgis.raster.functions" module in the ArcGIS API for Python 1.2 to demonstrate an example of a weighted overlay analysis.  This example attempts to identify areas in the State of Washington that are "natural" while also being easy to travel within based on the following criteria:

• elevation (lower is better)
• steepness of the terrain (flatter is better)
• degree of human alteration of the landscape (less is better)

The input data for this analysis includes a DEM (Digital Elevation Model), and a dataset showing the degree of human modification to the landscape.

In general, weighted overlay analysis can be divided into three steps:

1. Normalization: The pixels in the input raster datasets are reclassified to a common scale of numeric values based on their suitability according to the analysis criteria.
2. Weighting: The normalized datasets are assigned a percent influence based on their importance to the final result by multiplying them by values ranging from 0.0 - 1.0. The sum of the values must equal 1.0.
3. Summation: The sum of the weighted datasets is calculated to produce a final analysis result.

We'll begin by connecting to the GIS and accessing the data for the analysis.

## Connect to the GIS

In [1]:
`# import GIS from the arcgis.gis modulefrom arcgis.gis import GIS# Connect to the GIS.try:       web_gis = GIS("https://dev004543.esri.com/arcgis", 'djohnsonRA')       print("Successfully connected to {0}".format(web_gis.properties.name))except:       print("")`
`Enter password:········ Successfully connected to ArcGIS Enterprise A`

## Search the GIS for the input data for the analysis

### Human Modified Index

In [2]:
`# Search for the Human Modified Index imagery layer item by titleitem_hmi = web_gis.content.search('title:Human Modified Index', 'Imagery Layer')[0]item_hmi`
Out[2]:
Human Modified Index
A measure of the degree of human modification, the index ranges from 0.0 for a virgin landscape condition to 1.0, for the most heavily modified areas.Imagery Layer by djohnsonRA

### Elevation

In [3]:
`# Search for the DEM imagery layer item by titleitem_dem = web_gis.content.search('title:USGS NED 30m', 'Imagery Layer')[0]item_dem`
Out[3]:
USGS NED 30m
The National Elevation Dataset (NED) is the primary elevation data product of the USGS. This version was resampled to 30m from source data at 1/3 arc-second resolution and projected to an Albers Equal Area coordinate system.Imagery Layer by djohnsonRA

### Study area boundary and extent

In [4]:
`# Search for the Ventura County feature layer item by titleitem_studyarea = web_gis.content.search('title:State of Washington, USA',                                         'Feature Layer')[0]item_studyarea`
Out[4]:
State of Washington, USA
State of Washington, USAFeature Layer Collection by djohnsonRA
In [5]:
`# Get a reference to the feature layer from the portal itemlyr_studyarea = item_studyarea.layers[0]lyr_studyarea`
Out[5]:
`<FeatureLayer url:"https://dev004543.esri.com/server/rest/services/Hosted/Washington/FeatureServer/1">`

#### Get the coordinate geometry of the study area

In [6]:
`# Query the study area layer to get the boundary featurequery_studyarea = lyr_studyarea.query(where='1=1')# Get the coordinate geometry of the study area.# The geometry will be used to extract the Elevation and Human Modified Index data.geom_studyarea = query_studyarea.features[0].geometry# Set the spatial reference of the geometry.geom_studyarea['spatialReference'] = query_studyarea.spatial_reference`

#### Get the extent of the study area

In [7]:
`# Import the geocode functionfrom arcgis.geocoding import geocode# Use the geocode function to get the location/address of the study areageocode_studyarea = geocode('State of Washington, USA',                            out_sr= query_studyarea.spatial_reference)`
In [8]:
`# Get the geographic extent of the study area# This extent will be used when displaying the Elevation, Human Modified Index, # and final result data.extent_studyarea = geocode_studyarea[0]['extent']extent_studyarea`
Out[8]:
`{'xmax': -1451059.3770040546,  'xmin': -2009182.5321227335,  'ymax': 1482366.818700374,  'ymin': 736262.260048952}`

# Display the analysis data

### Human Modified Index

In [9]:
`# Get a reference to the imagery layer from the portal itemlyr_hmi = item_hmi.layers[0]# Set the layer extent to geographic extent of study area and display the data.lyr_hmi.extent = extent_studyarealyr_hmi`
Out[9]:

### Elevation

In [10]:
`# Get a reference to the imagery layer from the portal itemlyr_dem = item_dem.layers[0]# Set the layer extent to the geographic extent of study area and display the data.lyr_dem.extent = extent_studyarealyr_dem`
Out[10]:

### Slope (derived from elevation via the Slope raster function)

In [11]:
`# Import the raster functions from the ArcGIS API for Python (new to version 1.2!)from arcgis.raster.functions import *`
In [12]:
`# Derive a slope layer from the DEM layer using the slope functionlyr_slope = slope(dem=lyr_dem,slope_type='DEGREE', z_factor=1)# Use the stretch function to enhance the display of the slope layer.lyr_slope_stretch = stretch(raster=lyr_slope, stretch_type='StdDev', dra='true')# Display the stretched slope layer within the extent of the study area.lyr_slope_stretch.extent= extent_studyarealyr_slope_stretch`
Out[12]:

# Extract the data within the study area geometry

Use the Clip raster function to extract the analysis data from within the study area geometry

#### Human Modified Index

In [13]:
`# Extract the Human Modified Index data from within the study area geometryhmi_clipped = clip(raster=lyr_hmi, geometry=geom_studyarea)hmi_clipped`
Out[13]:

#### Elevation#Elevation

In [14]:
`# Extract the Elevation data from within the study area geometryelev_clipped = clip(raster=lyr_dem, geometry=geom_studyarea)elev_clipped`
Out[14]:

#### Slope#Slope

In [15]:
`# Extract the Slope data from within the study area geometryslope_clipped = clip(raster=lyr_slope, geometry=geom_studyarea)# Apply the Stretch function to enhance the display of the slope_clipped layer.slope_clipped_stretch = stretch(raster=slope_clipped, stretch_type='StdDev',                                 dra='true')slope_clipped_stretch`
Out[15]:

# Perform the analysis

### Step 1: Normalization

Use the Remap function to normalize each set of input data to a common scale of 1 - 9, where 1 = least suitable and 9 = most suitable.

In [16]:
`# Create a colormap to display the analysis results with 9 colors ranging # from red to yellow to green.clrmap=  [[1, 230, 0, 0], [2, 242, 85, 0], [3, 250, 142, 0], [4, 255, 195, 0],          [5, 255, 255, 0], [6, 197, 219, 0], [7, 139, 181, 0], [8, 86, 148, 0],         [9, 38, 115, 0]]`
In [17]:
`# Normalize the elevation dataelev_normalized = remap(raster=elev_clipped,                        input_ranges=[0,490, 490,980, 980,1470, 1470,1960, 1960,2450,                                       2450,2940, 2940,3430, 3430,3700, 3920,4100],                        output_values=[9,8,7,6,5,4,3,2,1], astype='U8')# Display color-mapped image of the reclassified elevation datacolormap(elev_normalized, colormap=clrmap) `
Out[17]:
In [18]:
`# Normalize the slope dataslope_normalized = remap(raster=slope_clipped,                                                  input_ranges=[0,1, 1,2, 2,3, 3,5, 5,7, 7,9, 9,12, 12,15,                                       15,100],                        output_values=[9,8,7,6,5,4,3,2,1],  astype='U8')  # Display a color-mapped image of the reclassified slope datacolormap(slope_normalized, colormap=clrmap)`
Out[18]:
In [19]:
`# Normalize the Human Modified Index datahmi_normalized = remap(raster=hmi_clipped,                                        input_ranges=[0.0,0.1, 0.1,0.2, 0.2,0.3, 0.3,0.4, 0.4,0.5,                                    0.5,0.6, 0.6,0.7, 0.7,0.8, 0.8,1.1],                                        output_values=[9,8,7,6,5,4,3,2,1],  astype='U8')# Display a color-mapped image of the reclassified HMI datacolormap(hmi_normalized, colormap=clrmap)`
Out[19]:

### Step 2: Weighting

Use the overloaded multiplication operator * to assign a weight to each normalized dataset based on their relative importance to the final result.

In [20]:
`# Apply weights to the normalized data using the overloaded multiplication # operator "*".# - Human Modified Index: 60%# - Slope: 25%# - Elevation: 15%hmi_weighted = hmi_normalized * 0.6slope_weighted = slope_normalized * 0.25elev_weighted = elev_normalized * 0.15`

### Step 3: Summation

Add the weighted datasets together to produce a final analysis result.

In [21]:
`# Calculate the sum of the weighted datasets using the overloaded addition # operator "+". result_dynamic = colormap(hmi_weighted + slope_weighted + elev_weighted,                           colormap=clrmap, astype='U8')result_dynamic`
Out[21]:

### The same analysis can also be performed in a single operation

In [22]:
`result_dynamic_one_op = colormap(       raster=       (      # Human modified index layer              0.60 * remap(raster=clip(raster=lyr_hmi, geometry=geom_studyarea),         input_ranges=[0.0,0.1, 0.1,0.2, 0.2,0.3, 0.3,0.4, 0.4,0.5,                       0.5,0.6, 0.6,0.7, 0.7,0.8, 0.8,1.1],                              output_values=[9,8,7,6,5,4,3,2,1])              +                # Slope layer              0.25 * remap(raster=clip(raster=lyr_slope, geometry=geom_studyarea),         input_ranges=[0,1, 1,2, 2,3, 3,5, 5,7, 7,9, 9,12, 12,15,                       15,100],                              output_values=[9,8,7,6,5,4,3,2,1])                +              # Elevation layer              0.15 * remap(raster=clip(raster=lyr_dem, geometry=geom_studyarea),         input_ranges=[-90,250, 250,500, 500,750, 750,1000, 1000,1500,                       1500,2000, 2000,2500, 2500,3000, 3000,5000],                             output_values=[9,8,7,6,5,4,3,2,1])           ),       colormap=clrmap,  astype='U8')result_dynamic_one_op`
Out[22]:

## Generate a persistent analysis result via distributed server based raster processing.

Portal for ArcGIS has been enhanced with the ability to perform distributed server based processing on imagery and raster data. This technology enables you to boost the performance of raster processing by processing data in a distributed fashion, even at full resolution and full extent.

You can use the processing capabilities of ArcGIS Pro to define the processing to be applied to raster data and perform processing in a distributed fashion using their on premise portal. The results of this processing can be accessed in the form of a web imagery layer that is hosted in their ArcGIS Organization.

In [23]:
`# Does the GIS support raster analytics?import arcgisarcgis.raster.analytics.is_supported(web_gis)`
Out[23]:
`True`
In [24]:
`# The .save() function invokes generate_raster from the arcgis.raster.analytics# module to run the analysis on a GIS server at the source resolution of the # input datasets and store the result as a persistent web imagery layer in the GIS.result_persistent = result_dynamic.save("NaturalAndAccessible_WashingtonState")result_persistent`
Out[24]:
NaturalAndAccessible_WashingtonState
Analysis Image Service generated from GenerateRasterImagery Layer by djohnsonRA
In [25]:
`# Display the persistent resultlyr_result_persistent = result_persistent.layers[0]lyr_result_persistent.extent = extent_studyarealyr_result_persistent`
Out[25]:
Data Credits:
Human Modified Index:
A measure of the degree of human modification, the index ranges from 0.0 for a virgin landscape condition to 1.0 for the most heavily modified areas. The average value for the United States is 0.375. The data used to produce these values should be both more current and more detailed than the NLCD used for generating the cores. Emphasis was given to attempting to map in particular, energy related development. Theobald, DM (2013) A general model to quantify ecological integrity for landscape assessment and US Application. Landscape Ecol (2013) 28:1859-1874 doi: 10.1007/s10980-013-9941-6
USGS NED 30m:
Data available from the U.S. Geological Survey. See USGS Visual Identity System Guidance for further details. Questions concerning the use or redistribution of USGS data should be directed to: ask@usgs.gov or 1-888-ASK-USGS (1-888-275-8747). NASA Land Processes Distributed Active Archive Center (LP DAAC) Products Acknowledgement: These data are distributed by the Land Processes Distributed Active Archive Center (LP DAAC), located at USGS/EROS, Sioux Falls, SD.
State of Washington: Esri Data & Maps

# Realistic Waters Effects in JavaScript Web Apps

Posted by rcarmichael-esristaff Jul 3, 2017

At last year's Developer SummitJesse van den Kieboom demonstrated how realistic water effects can be applied to a JavaScript based web application (see slides, demo and source).  The Prototype Lab modified Jesse's code to work with coastal inundation areas hosted in an AGOL feature service.  This sample is based on version 4.3 of the ArcGIS API for JavaScript and three.js.