Perform Analysis - Using out of the box tools

There are many out-of-the-box raster analysis tools available in ArcGIS API for Python. In this guide, we will demonstrate a few examples such as creating viewshed, interpolating points, converting raster to feature, and converting feature to raster.

Getting prepared

Let's first import the packages and make a connection to our GIS.

import arcgis
from arcgis.gis import GIS

gis = GIS(url='https://pythonapi.playground.esri.com/portal', username='arcgis_python', password='amazing_arcgis_123')

Creating viewshed

Viewshed indentifies the cells in an input raster that can be seen from one or more observation points or lines. Each cell in the output raster receives a value that indicates how many points or lines can be seen from each location. Viewshed is very useful when you want to know how visible objects might be, such as finding well-exposed places for communication towers.

Figure 1. An example of creating viewshed with a single observation point

In order to create viewshed, we need an input surface raster layer and an observer feature layer. In this example, we will use an elevation layer as the input raster and a few points representing schools as the observer layer.

items1 = gis.content.search("title: elevation, owner:portaladmin", item_type="Imagery Layer")
items1[0]
elevation
testImagery Layer by portaladmin
Last Modified: September 19, 2019
0 comments, 0 views
input_raster = items1[0].layers[0] # elevation/DEM raster layer
map1 = gis.map('Stowe, Vermont, USA')
map1
map1.add_layer(input_raster)
map1.legend = True
items2 = gis.content.search("title: schools, owner:portaladmin", item_type="Feature Layer")
items2[0]
schools
Schools in Stowe Vermont. This point feature class dataset is intended to be used as part of the Spatial Analyst tutorial.Feature Layer Collection by portaladmin
Last Modified: September 19, 2019
0 comments, 0 views
observer_layer = items2[0].layers[0] # school points as observer layer
map1.add_layer(observer_layer)

As you can see, we have added five schools on the map. Now we have the input layers in place, we are ready to create viewshed.

from arcgis.raster.analytics import create_viewshed

res = create_viewshed(input_elevation_surface = input_raster,
                      input_observer_features= observer_layer)
res
CreateViewshed_CGU80D
Analysis Image Service generated from CreateViewshedImagery Layer by portaladmin
Last Modified: September 24, 2019
0 comments, 0 views

Let's visualize the viewshed output.

map2 = gis.map('Stowe, Vermont, USA')
map2
map2.add_layer(res)

Let's enable the legend to make more sense of the viewshed result.

map2.legend = True

As we can see, there are five different colors, each representing a number. For example, red pixels means four schools can been seen from them.

Interpolating points

Interpolation predicts values for cells in a raster from a limited number of sample data points. It can be used to predict unknown values for any geographic point data, such as elevation, rainfall, chemical concentrations, and noise levels.

Figure 2. An example of interpolating rainfall point data

The interpolation tools are generally divided into deterministic and geostatistical methods. The deterministic interpolation methods assign values to locations based on the surrounding measured values and on specified mathematical formulas that determine the smoothness of the resulting surface. The deterministic methods include IDW (inverse distance weighting), Natural Neighbor, Trend, and Spline.

The geostatistical methods are based on statistical models that include autocorrelation (the statistical relationship among the measured points). Because of this, geostatistical techniques not only have the capability of producing a prediction surface but also provide some measure of the certainty or accuracy of the predictions. Ordinary Kriging and Empirical Bayesian Kriging are geostatistical methods of interpolation.

In this example, we use ozone concentration point samples (measurements) in California to produce a continuous surface (map) predicting the values of ozone concentration for every location throughtout the state.

items3 = gis.content.search("title: O3_Sep06_3pm, owner:portaladmin", item_type="Feature Layer")
items3[0]
O3_Sep06_3pm
This dataset is used as an example of environmental data for geostatistical spatial data analysis using Geostatistical Analyst extension.Feature Layer Collection by portaladmin
Last Modified: September 20, 2019
0 comments, 0 views
ozone_layer = items3[0].layers[0]

Now we have the ozone point samples, we are ready to interpolate these points to create a continuous surface. By default, interpolate_points() uses Empirical Bayesian kriging under the hood.

from arcgis.raster.analytics import interpolate_points

res = interpolate_points(input_point_features = ozone_layer,
                         interpolate_field = "OZONE",
                         output_name = "interpolation_res1")
res
FunctionOutput(output_raster=<Item title:"interpolation_res1" type:Imagery Layer owner:portaladmin>, process_info=<IPython.core.display.HTML object>)

interpolate_points() can generate multiple raster outputs based on output_prediction_error parameter and hence it returns a named tuple that contains:

  • output_raster (the output_raster item description is updated with the process_info),
  • process_info (if run in a non-Jupyter environment, use process_info.data to get the HTML data) and
  • output_error_raster (if output_prediction_error is set to True).
res.output_raster # access the result item through the `output_raster` property
interpolation_res1
Analysis Image Service generated from InterpolatePointsImagery Layer by portaladmin
Last Modified: September 24, 2019
0 comments, 0 views
map3 = gis.map('California, USA')
map3
map3.add_layer(res.output_raster)
map3.add_layer(ozone_layer)
map3.legend = True

Converting feature to raster

You can convert any feature class (polygon, polyline, or point) to a raster with convert_feature_to_raster(). When you convert points, cells are usually given the value of the points found within each cell. Cells that do not contain a point are given the value of NoData. If more than one point is found in a cell, the cell is given the value of the first point it encounters when processing. Using a smaller cell size during conversion may alleviate this. More information can be found at here.

In this guide, we are going to take the same ozone data from the example above and convert it into a raster layer. value_field is the field that will be used to assign values to the output raster.

from arcgis.raster.analytics import convert_feature_to_raster

res = convert_feature_to_raster(input_feature = ozone_layer,
                                value_field = "OZONE",
                                output_cell_size = {"distance":3000,"units": "meters"})
res
ConvertFeatureToRaster_BT5DP0
Analysis Image Service generated from ConvertFeatureToRasterImagery Layer by portaladmin
Last Modified: September 20, 2019
0 comments, 0 views

As you can see, the returned item is an Imagery Layer. Now let's add it onto a map.

map4 = gis.map('California, USA')
map4
map4.add_layer(res)
raster_layer = res.layers[0]

As we can see, cells that do not contain a point are given the value of NoData, which are transparent.

Converting raster to feature

In the meanwhile, we can convert raster data to feature class as well. convert_raster_to_feature() supports all polygon, polyline, and point data. Let's take the result from the last step - raster_layer and convert it back into a feature class.

from arcgis.raster.analytics import convert_raster_to_feature

res = convert_raster_to_feature(input_raster = raster_layer,
                                output_type = "Point")
res
RasterToFeature_9BA39C
RasterToFeature_9BA39CFeature Layer Collection by portaladmin
Last Modified: September 20, 2019
0 comments, 0 views

Let's load the result into a spatially enabled dataframe.

import pandas as pd
from arcgis.features import GeoAccessor, GeoSeriesAccessor

sdf = pd.DataFrame.spatial.from_layer(res.layers[0])
len(sdf)
167

We can also dive into the attributes of these 167 points.

sdf.head()
SHAPEgrid_codeobjectidpointid
0{"x": -220664.8320000004, "y": 415180.40939999...0.043011
1{"x": -352664.8320000004, "y": 313180.40940000...0.021022
2{"x": -202664.8320000004, "y": 283180.40939999...0.039433
3{"x": -133664.83199999854, "y": 280180.4093999...0.050044
4{"x": -193664.8320000004, "y": 274180.40939999...0.045055

As can be seen, we have got 167 points with different geographic coordinates. Notice that we have a grid_code column that represents the cell value when the data was in raster form.

In this guide, we have demonstrated how to create viewshed, interpolate points, convert raster to feature, and convert feature to raster, but these are just a few examples. There are many other out-of-the-box tools available and can be found at API reference.

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