Raster analysis - advanced concepts
Raster functions are lightweight and process only the pixels visible on your screen, in memory, without creating intermediate files. They are powerful because you can chain them together and apply them on huge rasters and mosaics on the fly. This guide will introduce you to the different capabilities of Imagery Layers and how to apply raster functions. As a sample, it uses Landsat imagery data.
Note: This guide requires bokeh python plotting library. Install it in your conda environment using the command below. This command should be typed in the terminal, not in the notebook
conda install bokeh
Table of contents
- Access Landsat imagery
- Interactive raster processing in Jupyter notebook
- Setting an area of interest
- Exporting images from Imagery Layer
- Exporting images from an imagery layer to which a raster function has been applied
- Vegetation index
- Custom bands
- Image attributes
- Spectral profile from the sampled values at a location
- Clipping to an area of interest
- Select images by where clause, geometry and the time range
- Resolving overlapping pixels in selected rasters
- Change detection
- Persisting your analysis for visualization or analysis
Access Landsat imagery¶
ArcGIS Online provides multispectral landsat imagery layer, that we'll be using for this tutorial. Let's connect to ArcGIS Online and query for the "Landsat Multispectral" imagery layer:
from arcgis.gis import GIS
gis = GIS("your_enterprise_profile")
landsat_item = gis.content.search(query='Landsat GLS Multispectral AND owner:esri_livingatlas',
item_type='Imagery Layer', outside_org=True)[0]
landsat_item
View Landsat imagery layer item description¶
from IPython.display import HTML
HTML(landsat_item.description)
Access the layers available with the Landsat Imagery Layer item¶
This imagery layer item contains the imagery layer that we'll be using for this tutorial. Let's save a reference to the layer in the landsat
variable. Querying the variable will in the Jupyter notebook will quickly render it as an image
landsat = landsat_item.layers[0]
landsat
Explore different wavelength bands¶
import pandas as pd
pd.DataFrame(landsat.key_properties()['BandProperties'])
Visualize the layer in the map widget¶
Let's visualize the layer by creating a map widget around our area of interest and adding the Landsat Imagery layer to it:
m = gis.map('Redlands, CA')
m
m.add_layer(landsat)
Apply built-in raster functions¶
The multispectral imagery layer can be rendered using several different raster functions (also known as raster function templates). Each raster function template incorporates different satellite bands to highlight different land cover features. Obtain list of predefined raster function templates defined by the service backing the imagery layer:
for rasterfunc in landsat.properties.rasterFunctionInfos:
print(rasterfunc.name)
Let's apply the 'Color Infrared with DRA' raster function to visualize the color infrared view. This can be done using the apply
function in the arcgis.raster.functions
module. This function applies a server defined raster function template, given it's name, to the Imagery layer.
from arcgis.raster.functions import apply
color_infrared = apply(landsat, 'Color Infrared with DRA')
m = gis.map('Redlands, CA')
m.add_layer(color_infrared)
m
A true-color image uses red, green and blue satellite bands to create an image that looks like a photograph. The color infrared view, on the other hand, uses the near infrared, red and green satellite bands to create an image. As a result, vegetation appears bright red in the image above.
Interactive raster processing in Jupyter Notebook¶
Imagery layers can be queried in the Jupyter Notebok. This displays a raster representation of the layer as an image:
color_infrared
Setting an area of interest¶
When working with Imagery layers, you are typically working in an area of interest. You can set the extent of the Imagery layer to that area of interest, and query it to visualize the imagery layer with that extent within the notebook:
from arcgis.geocoding import geocode
area = geocode('Redlands, CA', out_sr=landsat.properties.spatialReference)[0]
color_infrared.extent = area['extent']
color_infrared
Since we will be using the landsat layer further down in this notebook, let's set it's extent to our area of extent as well:
landsat.extent = area['extent']
landsat
Exporting Images from Imagery Layer¶
In addition to visualizing the imagery layers in the Jupyter Notebook, or using a map widge, they can be exported using the export_image
method provided by ImageryLayers:
from IPython.display import Image
img = landsat.export_image(bbox=area['extent'], size=[1200,450], f='image')
Image(img)
The exported image can be saved to a file by specifying a folder and file name:
savedimg = landsat.export_image(bbox=area['extent'], size=[1200,450], f='image', save_folder='.', save_file='img.jpg')
savedimg
from IPython.display import Image
Image(filename=savedimg, width=1200, height=450)
Exporting images from an imagery layer to which a raster function has been applied¶
The above techniques can be used to visualize imagery layers to which raster functions (or a chain of raster functions) have been applied:
color_infrared = apply(landsat, 'Color Infrared with DRA')
img = color_infrared.export_image(bbox=area['extent'], size=[1200, 450], f='image')
Image(img)
Vegetation Index¶
An index combines two or more wavelengths to indicate the relative abundance of different land cover features, like vegetation or moisture.
Let's use the 'NDVI Colorized' raster function to visualize healthy, green vegetation:
ndvi_colorized = apply(landsat, 'NDVI Colorized')
ndvi_colorized
NDVI is an index that indicates the presence of healthy, green vegetation (seen in green above).
Custom Bands¶
You can also create your own indexes and band combinations, as well as specify stretch and gamma values to adjust the image contrast.
The code below first extracts the [3 (Red), 2 (Green), 1 (Blue)] bands using the extract_bands
function and passes it's output to the stretch
function to enhance the image:
from arcgis.raster.functions import stretch, extract_band
naturalcolor = stretch(extract_band(landsat, [3,2,1]),
stretch_type='percentclip', min_percent=0.1, max_percent=0.1, gamma=[1, 1, 1], dra=True)
naturalcolor
This is a true-color image (similar to a photograph), created using the red, green, and blue satellite bands. Notice how much easier it was to find healthy vegetation using the NDVI vegetation index, compared to the true-color image here.
Image Attributes¶
The get_samples()
method returns pixel values of the source data (i.e before any rendering rules or raster functions are applied) for a given geometry as well as image attributes, like Scene ID, Acquisition Date, and Cloud Cover.
import arcgis
g = arcgis.geometry.Geometry(area['extent'])
samples = landsat.get_samples(g, sample_count=50,
out_fields='AcquisitionDate,OBJECTID,GroupName,Category,SunAzimuth,SunElevation,CloudCover')
samples[0]
import datetime
value = samples[0]['attributes']['AcquisitionDate']
datetime.datetime.fromtimestamp(value /1000).strftime("Acquisition Date: %d %b, %Y")
Here's the same information visualized as a table using Pandas:
pd.DataFrame(samples[0]['attributes'], index=[0])
Spectral profile from the sampled values at a location¶
The get_samples
method can also estimate which type of land cover most closely matches any point on the map by comparing their spectral profiles. The example below adds a 'click handler' to the map widget. When the map is clicked, a spectral profile for the clicked location is displayed using the bokeh
charting library.
m = gis.map('Redlands, CA')
m
m.add_layer(landsat)
import bokeh
from bokeh.models import Range1d
from bokeh.plotting import figure, show, output_notebook
from IPython.display import clear_output
output_notebook()
def get_samples(mw, g):
clear_output()
m.draw(g)
samples = landsat.get_samples(g, pixel_size=30)
values = samples[0]['value']
vals = [float(int(s)/100000) for s in values.split(' ')]
x = ['1','2', '3', '4', '5', '6', '7', '8']
y = vals
p = figure(title="Spectral Profile", x_axis_label='Spectral Bands', y_axis_label='Data Values', width=600, height=300)
p.line(x, y, legend="Selected Point", line_color="red", line_width=2)
p.circle(x, y, line_color="red", fill_color="white", size=8)
p.y_range=Range1d(0, 1.0)
show(p)
print('Click anywhere on the map to plot the spectral profile for that location.')
m.on_click(get_samples)
Clipping to an area of interest¶
Imagery layers can be clipped to an area of interest using the clip
raster function. The clipping geometry can be obtained from a feature layer or spatial dataframe as well. Here's a simple example of clipping the landsat layer:
from arcgis.geometry import Geometry, buffer
The unit
parameter indicates the units for calculating each buffer distance. If not specified, the units are derived from bufferSR
. If bufferSR
is not specified, the units are derived from in_sr
. In the example below, we are calling the buffer
method with unit
set to be a constant string "9001"
:
poly = buffer(geometries=[Geometry(area['location'])],
in_sr=102100, distances=6000,
unit="9001")[0]
Equivalently we can use a constant Enum
, e.g. LengthUnits.METER
, or the integer 9001
to specify the same requirement.
from arcgis.geometry import LengthUnits
geom_buffered_b = buffer(geometries = [polyline1],
in_sr=3857,
distances=50,
unit=LengthUnits.METER,
out_sr=None,
buffer_sr=None,
union_results=None,
geodesic=None,
gis=gis)
from arcgis.raster.functions import clip
redclip = clip(landsat, poly)
The clipped layer can also be added to a map widget:
m = gis.map('Redlands, CA')
m
m.add_layer(redclip)
Select images by where clause, geometry and time range¶
You can select images by attributes using a where clause as well as using spatial and temporal filters, using the filter_by
method.
The code snippet below limits the images available to those that have less than 10% cloud cover and which intersect with our area of interest:
from arcgis.geometry.filters import intersects
selected = landsat.filter_by(where="(Category = 1) AND (CloudCover <=0.10) AND (WRS_Row = 36)",
geometry=intersects(area['extent']))
You can query the filtered rasters as a FeatureSet:
fs = selected.query(out_fields="AcquisitionDate, GroupName, Best, CloudCover, WRS_Row, Month, Name",
return_geometry=True,
return_distinct_values=False,
order_by_fields="AcquisitionDate")
Attributes of the selected rasters can be queried using a Pandas dataframe using the 'df' property of the FeatureSet.
df = fs.sdf
df.head()
Looking at the shape of the dataframe we see that 46 scenes match the specified criteria:
df.shape
The footprints of the rasters matching the criteria can be drawn using the map widget:
df['Time'] = pd.to_datetime(df['AcquisitionDate'], unit='ms')
df['Time'].head(10)
Resolving overlapping pixels in selected rasters¶
When a setof rasters are selected by filtering an Imagery layer, they may have overlapping pixels.
The Imagery layer has methods like first(), last(), min(), max(), mean(), blend() and sum() to resolve overlap pixel values from first or last raster, use the min, max or mean of the pixel values, or blend them:
Here's the same information visualized using the map widget. This shows the selected rasters covering only our area of interest:
m = gis.map('Redlands, CA', 7)
display(m)
m.add_layer(selected.last())
m = gis.map('Redlands, CA', 7)
display(m)
m.add_layer(selected.first())
Change Detection¶
It's quite common to compare images of the same area from two different times. The example below selects the old and new images using the filter_by()
method ans specifying an object id. The object id can be obtained from the query()
method described above.
old = landsat.filter_by('OBJECTID=1139')
new = landsat.filter_by('OBJECTID=463490')
Difference Image¶
Difference Image mode illustrates all the changes in NDVI (vegeration index) between the two dates:
increases are shown in green, and decreases are shown in magenta.
from arcgis.raster.functions import *
diff = stretch(composite_band([ndvi(old, '5 4'),
ndvi(new, '5 4'),
ndvi(old, '5 4')]),
stretch_type='stddev', num_stddev=3, min=0, max=255, dra=True, astype='u8')
diff
Applying a threshold mask¶
The difference can also be computed using map algebra, as shown below:
ndvi_diff = ndvi(new, '5 4') - ndvi(old, '5 4')
ndvi_diff
However, in the image above is hard to visualize in which areas the vegetation index changed by a specified threshold. The example below renders the areas where the change is above the specified threshold and renders it using a green color:
threshold_val = 0.1
masked = colormap(remap(ndvi_diff,
input_ranges=[threshold_val, 1],
output_values=[1],
no_data_ranges=[-1, threshold_val], astype='u8'),
colormap=[[1, 124, 252, 0]], astype='u8')
Image(masked.export_image(bbox=area['extent'], size=[1200,450], f='image'))
The difference image and threshold mask can be visualized using the map widget:
m = gis.map('Redlands, CA')
m
m.add_layer(diff)
m.add_layer(masked)
Persisting your analysis for visualizaion or analysis¶
The save()
method on ImageryLayer class persists this imagery layer to the GIS as an Imagery Layer item. If for_viz
parameter is True, a new Item is created that uses the applied raster functions for visualization at display resolution using on-the-fly image processing. If for_viz
is False, distributed raster analysis is used for generating a new raster information product by applying raster functions at source resolution across the extent of the output imagery layer.
In the example below, the threshold mask is being saved as item for visualization:
lyr = masked.save('Test_viz_layer3', for_viz=True)
lyr