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(profile="your_organization_profile")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: Stowe_DEM AND owner:api_data_owner", item_type="Imagery Layer")
items1[0]
Imagery Layer by api_data_owner
input_raster = items1[0].layers[0] # elevation/DEM raster layermap1 = gis.map('Stowe, Vermont, USA')
map1
map1.content.add(input_raster)map1.legend.enabled = Trueitems2 = gis.content.search("title: schools, owner:portaladmin", item_type="Feature Layer")
items2[0]
Feature Layer Collection by portaladmin
observer_layer = items2[0].layers[0] # school points as observer layer
map1.content.add(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
Imagery Layer by portaladmin
Let's visualize the viewshed output.
map2 = gis.map('Stowe, Vermont, USA')
map2
map2.content.add(res)Let's enable the legend to make more sense of the viewshed result.
map2.legend.enabled = TrueAs 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]
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")
resFunctionOutput(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` propertymap3 = gis.map('California, USA')
map3
map3.content.add(res.output_raster)
map3.content.add(ozone_layer)map3.legend.enabled = TrueConverting 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
Imagery Layer by portaladmin
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.content.add(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
Feature Layer Collection by portaladmin
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()| SHAPE | grid_code | objectid | pointid | |
|---|---|---|---|---|
| 0 | {"x": -220664.8320000004, "y": 415180.40939999... | 0.0430 | 1 | 1 |
| 1 | {"x": -352664.8320000004, "y": 313180.40940000... | 0.0210 | 2 | 2 |
| 2 | {"x": -202664.8320000004, "y": 283180.40939999... | 0.0394 | 3 | 3 |
| 3 | {"x": -133664.83199999854, "y": 280180.4093999... | 0.0500 | 4 | 4 |
| 4 | {"x": -193664.8320000004, "y": 274180.40939999... | 0.0450 | 5 | 5 |
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.