Identifying natural and accessible areas using weighted overlay analysis

Many people vacation to scenic areas free from everyday noise and congestion. However, many scenic areas can be difficult to reach and challenging to navigate once there. Many places also vary in their degree of remoteness from everyday human activity that travelers like to escape. This sample identifies areas in the State of Washington that are more "natural" and easy to get to and visit based on the following criteria:

  • elevation (lower elevations are easier to travel)
  • steepness of the terrain (lower slopes are easier to travel)
  • degree of human alteration of the landscape (less altered landscapes are more natural)

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

Weighted overlay analysis

The weighted overlay is a standard GIS analysis technique often used for solving multicriteria problems such as generating surfaces representing site-suitability and travel-cost. Weighted overlay is used when a number of factors of variying importance should be considered to arrive at a final decision. This sample shows how raster anlaysis and raster arithmetic can be used to perform such analysis to solve spatial problems. The graphic below explains the logic behind weighted overlay, refer to this help for a detailed review of weighted overlay analysis

In the illustration, the two input rasters have been reclassified to a common measurement scale of 1 to 3. Each raster is assigned a percentage influence. The cell values are multiplied by their percentage influence, and the results are added together to create the output raster. For example, consider the upper left cell. The values for the two inputs become (2 * 0.75) = 1.5 and (3 * 0.25) = 0.75. The sum of 1.5 and 0.75 is 2.25. Because the output raster from Weighted Overlay is integer, the final value is rounded to 2.

Connect to the GIS

In [1]:
# import GIS from the arcgis.gis module
from arcgis.gis import GIS
In [2]:
# Connect to the GIS
gis = GIS('https://python.playground.esri.com/portal', 'arcgis_python', 'amazing_arcgis_123')
print("Successfully connected to {0}".format(gis.properties.name))
Successfully connected to Python Playground

Access the data for analysis

In [3]:
# Search for the Human Modified Index imagery layer item by title
hmna_item = gis.content.search('title:Human Modification for North America', 'Imagery Layer')[0]
hmna_item
Out[3]:
Human Modification for North America
Human Modification for North America transboundary study area.Imagery Layer by arcgis_python
Last Modified: February 21, 2018
0 comments, 0 views
In [4]:
# Search for the DEM imagery layer item by title
elev_item = gis.content.search('title:National Elevation Dataset (NED)', 'Imagery Layer')[0]
elev_item
Out[4]:
National Elevation Dataset (NED)
National Elevation Dataset (NED) for CONUS resampled and reprojected from source data in NAD 83 at 1 arc-second resolution to an Albers Equal Area Conic projection at 270m resolution.Imagery Layer by arcgis_python
Last Modified: February 21, 2018
0 comments, 0 views

Get the study area geometry using data from the Living Atlas

In [5]:
# Access the USA States item from the Living Atlas using the item id value
states_item = gis.content.search("id:99fd67933e754a1181cc755146be21ca", 'Feature Layer', outside_org=True)[0]
states_item
Out[5]:
USA States (Generalized)
This layer presents the 50 states and the District of Columbia of the United States. This layer is available at no cost to users with an ArcGIS Online subscription.Feature Layer Collection by esri_livingatlas
Last Modified: October 06, 2016
0 comments, 0 views
In [6]:
# Access the USA states feature layer
states_lyr = states_item.layers[0]

We choose the State of Washington as our study area for this example.

In [7]:
# Access the feature for the State of Washington
study_area_query = states_lyr.query("STATE_NAME='Washington'", return_geometry=True)
In [8]:
# Get the geometry of the State of Washington feature.
# We will use this geometry to extract the input data for the study area.
study_area_geom= study_area_query.features[0].geometry
study_area_geom['spatialReference'] = study_area_query.spatial_reference

Get the coordinates of the study area extent using geocoding

In [9]:
# Import the geocode function
from arcgis.geocoding import geocode
# Use geocoding to get the location of the study area in the spatial reference of the input data for the analysis.
study_area_gcd = geocode(address='State of Washington, USA', out_sr=hmna_item.layers[0].extent['spatialReference'])
# Get the geographic extent of the study area.
# This extent will be used for displaying the input data and output results.
study_area_extent = study_area_gcd[0]['extent']
study_area_extent
Out[9]:
{'xmax': -1451059.3770040546,
 'xmin': -2009182.5321227335,
 'ymax': 1482366.818700374,
 'ymin': 736262.260048952}

Preview the analysis data

Human Modification for North America

In [10]:
# Get a reference to the imagery layer from the portal item
hmna_lyr = hmna_item.layers[0]
# Set the layer extent to geographic extent of study area and display the data.
hmna_lyr.extent = study_area_extent
hmna_lyr
Out[10]:

Elevation

In [11]:
# Get a reference to the imagery layer from the portal item
elev_lyr = elev_item.layers[0]
# Set the layer extent to the geographic extent of study area and display the data.
elev_lyr.extent = study_area_extent
elev_lyr
Out[11]:

Slope (derived from elevation via the Slope raster function)

In [13]:
# Import the raster functions from the API
from arcgis.raster.functions import *
In [14]:
# Derive a slope layer from the DEM layer using the slope function
slope_lyr = slope(dem=elev_lyr, slope_type='DEGREE', z_factor=1)
slope_lyr.extent = study_area_extent
# Use the stretch function to enhance the display of the slope layer.
stretch(raster=slope_lyr, stretch_type='StdDev', dra='true')
Out[14]: