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

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:

In [1]:

from arcgis.gis import GIS

In [2]:

gis = GIS('https://python.playground.esri.com/portal', 'arcgis_python', 'amazing_arcgis_123')

In [3]:

landsat_item = gis.content.search('Multispectral Landsat', 'Imagery Layer', outside_org=True)[0]

In [4]:

landsat_item

Out[4]:

Multispectral Landsat
Landsat 8 OLI, 30m Multispectral 8 band scenes with visual renderings and indices. Updated daily. Based on the Landsat on AWS collections.

Imagery Layer by esri_livingatlas
Last Modified: October 06, 2016
0 comments, 0 views

View Landsat imagery layer item description¶

In [5]:

from IPython.display import HTML
HTML(landsat_item.description)

Out[5]:

This world-wide dynamic image service provides access to Landsat 8 Operational Land Imager (OLI) scenes at 30m. Landsat 8 collects new multispectral scenes for each location on Earth every 16 days, assuming limited cloud coverage. Newest and near cloud-free scenes are displayed by default on top. Attributes include acquisition date and estimated cloud cover percentage. Most scenes collected since 1st January 2015 are included as well as approximately 5 scenes for each path/row from 2013 and 2014. The service also includes scenes from the Global Land Survey* (circa 2010, 2005, 2000, 1990, 1975).

The service provides access to all multispectral OLI bands as follows:

Band Num.	Description	Wavelength (micrometers)
1	Coastal aerosol	0.43 - 0.45
2	Blue	0.45 - 0.51
3	Green	0.53 - 0.59
4	Red	0.64 - 0.67
5	Near Infrared (NIR)	0.85 - 0.88
6	SWIR 1	1.57 - 1.65
7	SWIR 2	2.11 - 2.29
8	Cirrus (Note in OLI this is band 9)	1.36 - 1.38

The service includes various selectable functions that render different band combinations and indices. Each scene is processed on-the-fly from the original Landsat (L1T) scenes. Top of Atmosphere reflectance (scaled between 5000 to 55000) is computed and then the appropriate renderer is applied. The default display is Agriculture (bands 6,5,2) with Dynamic Range Adjustment (DRA), which is a band combination that displays vegetation as green. The DRA version of each band combination enables visualization of the full dynamic range of the images.

The service is time-enabled so applications can restrict the displayed scenes to a specific date range.

By default the layer is set to display and query only the best 3 (most recent cloud free) scene. The filter can be removed to gain access to all the scenes, but note that identify may take time to access all pixels in all scenes. By setting the filter to Best < QQQQ one can control to see the best N scenes. Where QQQQ=N*1million.

Additional filters can be set to restrict and order the scenes based on attributes.

Overviews exist with a pixel size of 300m and are updated weekly based on the best and latest imagery. As you zoom out beyond approximately 1:1M scale the overviews are displayed. To work with an individual scene(s) at all scales, use the lock raster functionality. ‘Lock Raster’ should only be used on the service for short periods of time, as the service and associated record Object IDs may change daily. Alternatively add a query filter to restrict the display to a specified scene or collection of scenes.

Similar services, including Panchromatic and Pansharpened Images, are also accessible.

This ArcGIS Server dynamic image service can be used in Web Maps and ArcGIS Desktop, Web and Mobile applications using the REST based Image services API.

Users can export images, but the exported area is limited to maximum of 2,000 columns x 2,000 rows per request.

Data Source: The imagery in these services is sourced from the U.S. Geological Survey (USGS) and the National Aeronautics and Space Administration (NASA). The data for these services reside on the Landsat Public Datasets hosted on the Amazon Web Service cloud. Users can access full scenes from https://github.com/landsat-pds/landsat_ingestor/wiki/Accessing-Landsat-on-AWS, or alternatively access http://landsatlook.usgs.gov to review and download full scenes from the complete USGS archive.

For more information on Landsat 8 images, see http://landsat.usgs.gov/landsat8.php.

*The Global Land Survey includes images from Landsat 1 through Landsat 7. Band numbers and band combinations differ from those of Landsat 8, but have been mapped to the most appropriate band as in the above table. For more information about the Global Land Survey, visit http://landsat.usgs.gov/science_GLS.php.

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

In [6]:

landsat = landsat_item.layers[0]
landsat

Out[6]:

Explore different wavelength bands¶

In [33]:

import pandas as pd

In [8]:

pd.DataFrame(landsat.key_properties()['BandProperties'])

Out[8]:

	BandIndex	BandName	DatasetTag	WavelengthMax	WavelengthMin
0	0	CoastalAerosol	MS	450	430
1	1	Blue	MS	510	450
2	2	Green	MS	590	530
3	3	Red	MS	670	640
4	4	NearInfrared	MS	880	850
5	5	ShortWaveInfrared_1	MS	1650	1570
6	6	ShortWaveInfrared_2	MS	2290	2110
7	7	Cirrus	MS	1380	1360

Let's visualize the layer by creating a map widget around our area of interest and adding the Landsat Imagery layer to it:

In [9]:

m = gis.map('Redlands, CA')
m

In [10]:

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:

In [11]:

for rasterfunc in landsat.properties.rasterFunctionInfos:
    print(rasterfunc.name)

Agriculture with DRA
Bathymetric with DRA
Color Infrared with DRA
Geology with DRA
Natural Color with DRA
Short-wave Infrared with DRA
Agriculture
Bathymetric
Color Infrared
Geology
Natural Color
Short-wave Infrared
NDVI Colorized
Normalized Difference Moisture Index Colorized
NDVI Raw
NBR Raw
None

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.

In [12]:

from arcgis.raster.functions import apply

In [13]:

color_infrared = apply(landsat, 'Color Infrared with DRA')

In [14]:

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:

In [16]:

color_infrared

Out[16]: