PCI Geomatics Help Center

How can we help you today?

Analysis Ready Data (ARD) - Vegetation Monitoring Workflow

PCI Geomatics -

New Analysis Ready Data (ARD) tools were added to Geomatica 2018 and are available in Geomatica Banff. With these tools you can create industry-leading workflows to produce data that is corrected geometrically and radiometrically at scale. You can use the tools to generate standardized illumination conditions and surface-reflectance products to perform meaningful multi-temporal analysis.

 Raw sensor imagery is affected by various factors including sensor calibration and degradation, illumination geometry (different solar angle, Earth-sun distance), atmosphere (absorption, scattering) and Intensity variations due to solar geometry and slope/aspect. Each of these effects are removed in the steps of the ARD workflow. 

This tutorial outlines an Analysis Ready Data workflow that can be used to produce outputs for vegetation monitoring. The data used is Landsat-8 L1TP. Landsat-8 data from the USGS EarthExplorer website - https://earthexplorer.usgs.gov/ . The imagery, DEM and script for this tutorial are available from here: https://pcigeomatics.sharefile.com/d-sf73df0252dd421d9

Each of the following processing steps will be completed in the ARD workflow:

  1. Per-pixel solar and sensor calculations - solviewzaz
  2. TOA Reflectance Calibration - dn2reflectance
  3. MRA Fusion - mrafusion
  4. Up sample per-pixel angle file - reproj
  5. Spectral Pre-Classification - speclass
  6. Topographic Normalization - toposolnorm
  7. Generate outputs for vegetation analysis
    • Tasseled Cap Transformation – tassel
    • Vegetation Indices - vegindex

In this tutorial, these algorithms were run in python, however you can also run each of these algorithms individually in Focus > Algorithm Librarian or in EASI. 

Note that more information on creating python scripts using Geomatica libraries is available from the Developer Zone. If you are just starting out in Python the Getting Started page is quite helpful. There are also a number of Geomatica Developer tutorials.

All of Geomatica algorithm parameters are outlined in the Geomatica Help. You can check the help for detailed information on each parameter.

mceclip0.png mceclip1.png
Raw Landsat-8 MS image

Raw Landsat-8 PAN image

 

Set up the python script

You will first need to import the Geomatica library, Geomatica Python API and Python OS module

Additionally a helper function is defined at the beginning of the script – get_geo_info. This function reads in a file and outputs a three-element tuple with the map units string, upper left coordinate and lower right coordinate of the image. This function will be run later in the script to get the geographic information required to run reprojection (REPROJ).

from pci import algo
import os
from pci.api import datasource as ds
from pci.api import cts

# Function to obtain resolution and image coordinates for use when upsampling the solviewzaz multispectral output
def get_geo_info(input_file):

## get_geo_info function reads in a file and outputs a three-element tuple with Map Units String, Upper Left Coordinate and
## Lower Right Coordinate
##
## :param input_file: path to input raster file
## :return: six-element tuple with Map Units String, Upper Left Coordinate, Lower Right Coordinate, Width, Height,
## and resolution

# Open the dataset
with ds.open_dataset(input_file, ds.eAM_READ) as dataset:

map_units = cts.crs_to_mapunits(dataset.crs) # Get map units string

# Get dataset height and width in pixels
height = dataset.height

width = dataset.width

# Get upper left and lower right pixels
upper_left = dataset.geocoding.raster_to_map(0, 0)

lower_right = dataset.geocoding.raster_to_map(width, height)
resolution = [dataset.geocoding.resolution[0], dataset.geocoding.resolution[1]]

return map_units, upper_left, lower_right, width, height, resolution

We will now set the input and output variables that point to the raw data and DEM. The call to os.path.exists ensures that the output directory exists and if it does not, creates that directory.

# Define Landsat imagery inputs
ms_input = r"E:\raw\LC08_L1TP_045028_20180928_20181009_01_T1_MTL.txt-MS"
pan_input = r"E:\raw\LC08_L1TP_045028_20180928_20181009_01_T1_MTL.txt-PAN"
dem_input = r"E:\DEM.pix"

# Define output folder for results
output_dir = r"E:\ARD_Output"

# Make sure that output directory exists
if not os.path.exists(output_dir):

os.makedirs(output_dir)

Before starting to run the ARD algorithms we will also import the raw images into .pix format.

# Create live links to Landsat MTL file for MS and PAN images
ms_link = os.path.join(output_dir, "L8_MS.pix")

pan_link = os.path.join(output_dir, "L8_PAN.pix")
algo.link(fili=ms_input, filo=ms_link)
algo.link(fili=pan_input, filo=pan_link)

Per-pixel solar and sensor calculations: SOLVIEWZAZ

The first algorithm you need to run is SOLVIEWZAZ which derives per-pixel solar and sensor view zenith and azimuth angles. The output is required for reflectance calibration and topographic normalization.

For Landsat-8 the per-pixel solar and sensor viewing angle channels are calculated using Landsat Angles Creation Tools (l8_angles) and the solar Illumination and Sensor Viewing Angle Coefficient File provided with Landsat-8

# Derive per-pixel solar and sensor view zenith and azimuth angles
sol_ms = os.path.join(output_dir, "solviewzaz_MS.pix")

sol_pan = os.path.join(output_dir, "solviewzaz_PAN.pix")
algo.solviewzaz(fili=ms_link, filo=sol_ms, chntype="32R")
algo.solviewzaz(fili=pan_link, filo=sol_pan, chntype="32R")
mceclip2.png mceclip3.png
MS SOLVIEWZAZ result

PAN SOLVIEWZAZ result

TOA Reflectance Calibration: DN2Reflectance

DN2REFLECTANCE converts from radiometrically calibrated input data to absolute reflectance pixel values. This conversion requires that the solar effects are removed, and therefore these input data types require the knowledge of solar zenith angles. We will use the solviewzaz file as an input to DN2REFLECTANCE.

The Landsat-8 data contains relative reflectance data, so the angle file from SOLVIEWZAZ is required as an input. 

# Convert to absolute reflectance pixel values with radiometric calibration information
dn2ref_ms = os.path.join(output_dir, "dn2reflectance_MS.pix")

dn2ref_pan = os.path.join(output_dir, "dn2reflectance_PAN.pix")
algo.dn2reflectance(fili=ms_link, fileang=sol_ms, filo=dn2ref_ms)
algo.dn2reflectance(fili=pan_link, fileang=sol_pan, filo=dn2ref_pan)

MRA Pan-Sharpening: MRAFUSION

MRAFUSION applies an image-fusion algorithm to increase the resolution of multispectral (color) imagery by using a high-resolution, panchromatic (black and white) image. Known as multi-resolution analysis (MRA), this image-fusion technique applies wavelet decomposition. MRAFUSION produces similar-to-superior sharpening results in comparison to other image fusion techniques (ex. PANSHARP) while preserving the spectral characteristics of the original images.

# Create high-resolution color images by fusing black-and-white panchromatic (PAN) and multispectral(MS) color images.
mrafus = os.path.join(output_dir, "mrafusion.pix")

algo.mrafusion(fili=dn2ref_ms, dbic=[1, 2, 3, 4, 5, 6, 7, 8], dbic_ref=[2, 3, 4, 5], fili_pan=dn2ref_pan, dbic_pan=[1],
filo=mrafus)
mceclip2.png  mceclip3.png

Original MS image

MRAFUSION result

Up-sample per-pixel angle file: REPROJ

Since the MS image has now been pansharpened to a higher resolution we also need to up sample the multispectral SOLVIEWZAZ layer to the same pansharpened resolution. Both the pansharpened image and up-sampled SOLVIEWZAZ layers will be required in the topographic normalization (TOPOSOLNORM) step.

# Up sample solviewzaz MS to match the PAN resolution and extents
sol_ms_reproj = os.path.join(output_dir, "solviewzaz_MS_reproj.pix")

geo_info = get_geo_info(pan_link)
algo.reproj(fili=sol_ms, dbic=[1, -32], dbsl=[1], filo=sol_ms_reproj, dbsz=[geo_info[3], geo_info[4]], pxsz=geo_info[5],
maxbnds="NO", mapunits=geo_info[0], ulx=str(geo_info[1][0]), uly=str(geo_info[1][1]),
lrx=str(geo_info[2][0]), lry=str(geo_info[2][1]), resample="NEAR")

Spectral Pre-Classification: SPECLASS

SPECLASS creates a spectral classification from optical data by converting it to TOA reflectance and extracting per-pixel features and indices. The output is a spectral pre-classification containing relevant metadata for radiometric normalization.  SPECLASS can be used for the generation of a general quality layer to guide further radiometric and geometric processing. The resulting cloud and water mask can be used for targeted ground control point and tie point collection leading to better image triangulation.

The output from this call to SPECLASS will be used as the classification layer in the topographic normalization step.

# Create a spectral classification from optical data by converting it to TOA reflectance and extracting per-pixel features and indices.
spec_fus = os.path.join(output_dir, "speclass.pix")

algo.speclass(fili=mrafus, filo=spec_fus)
mceclip7.png mceclip5.png
Legend

SPECLASS result

 

Topographic Normalization: TOPOSOLNORM

TOPOSOLNORM corrects pre-classified radiometric values based on variable incidence angles. This implementation of topographic normalization accounts for surfaces with different BRDF properties but is less sensitive to lower quality DEMs while accounting for various levels of skylight that is not eliminated during atmospheric correction. You can see in the results that areas of varied elevation such as mountains appear quite different. In the screenshot below the shadow that appears in the original imagery is significantly reduced by toposolnorm. 

In this tutorial TOPOSOLNORM will be run using the fused imagery from MRAFUISION, user provided DEM, classification file from SPECLASS and the up sampled per-pixel angles file from SOLVIEWZAZ.

# Correct pre-classified radiometric values based on variable incidence angles
toposol = os.path.join(output_dir, "toposolnorm.pix")

algo.toposolnorm(fili=mrafus, filedem=dem_input, filecls=spec_fus, fileang=sol_ms_reproj, filo=toposol)
mceclip8.png  mceclip9.png
MRAFUSION result TOPOSOLNORM result

OUTPUTS

With the final TOPOSOLNORM image, the data is ready for different analysis options. Brightness, greenness, and wetness rasters can be obtained, and a selection of vegetation indices can also be calculated for the imagery.

Output Option A: TASSEL

TASSEL calculates the tasseled-cap transformation for the sensor and bands specified. To calculate the tasseled-cap transformation, the input file must contain the specific bands necessary. The number of bands in the file must range from three through six (Blue/Green/Red/NIR/SWIR1.6/SWIR2.2). You must also specify the necessary bands (channels) in the correct order.

The output will always be three channels: brightness, greenness, and wetness. The tasseled-cap transformation for brightness is written to the first output channel, greenness to the second channel, and wetness to the third channel of the output file.

# Ouput A: Calculate the tasseled-cap transformation for the sensor and bands specified
tassel_topo = os.path.join(output_dir, "topo_tassel.pix")

algo.tassel(fili=toposol, filo=tassel_topo)
mceclip13.png mceclip10.png

TASSEL Greeness (brighter pixels indicate higher greenness)

TASSEL Wetness (lighter areas indicate higher moisture)

Output Option B: VEGINDEX

VEGINDEX calculates one or more vegetation indices using the reflectance data (top of atmosphere (TOA) or absolute reflectance) in the input file, channels, or both. If the input data is DN values, they will be converted internally to TOA reflectance values for calculation. The results are written to a new file, or to the file and channels you specify.

To avoid scaling of the results produced by VEGINDEX, make sure they are written to 32-bit real channels. When file and band metadata information is available, VEGINDEX uses this information to select appropriate data for the input index or indices.

A full list of possible vegetation indices can be found here.

# Output B: Compute selected vegetation indices
vegindex_topo = os.path.join(output_dir, "topo_vegindex.pix")

algo.vegindex(fili=toposol, filo=vegindex_topo, index="NDVI")

 mceclip12.png

VEGINDEX NDVI result (lighter pixels indicate a greater amount of vegetation detected)

The final vegetation monitoring results are now generated. You could also run pixel-based or object based classification to further extract vegetation information. Additionally, you could run this same workflow on the same area of interest over various dates to track the change in vegetation over time.

Full Script

from pci import algo
import os
from pci.api import datasource as ds
from pci.api import cts

# Function to obtain resolution and image coordinates for use when upsampling the solviewzaz multispectral output
def get_geo_info(input_file):

## get_geo_info function reads in a file and outputs a three-element tuple with Map Units String, Upper Left Coordinate and
## Lower Right Coordinate
##
## :param input_file: path to input raster file
## :return: six-element tuple with Map Units String, Upper Left Coordinate, Lower Right Coordinate, Width, Height,
## and resolution

# Open the dataset
with ds.open_dataset(input_file, ds.eAM_READ) as dataset:

map_units = cts.crs_to_mapunits(dataset.crs) # Get map units string

# Get dataset height and width in pixels
height = dataset.height

width = dataset.width

# Get upper left and lower right pixels
upper_left = dataset.geocoding.raster_to_map(0, 0)

lower_right = dataset.geocoding.raster_to_map(width, height)
resolution = [dataset.geocoding.resolution[0], dataset.geocoding.resolution[1]]

return map_units, upper_left, lower_right, width, height, resolution


# Define Landsat imagery inputs
ms_input = r"E:\raw\LC08_L1TP_045028_20180928_20181009_01_T1_MTL.txt-MS"
pan_input = r"E:\raw\LC08_L1TP_045028_20180928_20181009_01_T1_MTL.txt-PAN"
dem_input = r"E:\DEM.pix"

# Define output folder for results
output_dir = r"E:\ARD_Output"

# Make sure that output directory exists
if not os.path.exists(output_dir):

os.makedirs(output_dir)

# Create live links to Landsat MTL file for MS and PAN images
ms_link = os.path.join(output_dir, "L8_MS.pix")

pan_link = os.path.join(output_dir, "L8_PAN.pix")
algo.link(fili=ms_input, filo=ms_link)
algo.link(fili=pan_input, filo=pan_link)

# Derive per-pixel solar and sensor view zenith and azimuth angles
sol_ms = os.path.join(output_dir, "solviewzaz_MS.pix")

sol_pan = os.path.join(output_dir, "solviewzaz_PAN.pix")
algo.solviewzaz(fili=ms_link, filo=sol_ms, chntype="32R")
algo.solviewzaz(fili=pan_link, filo=sol_pan, chntype="32R")

# Convert to absolute reflectance pixel values with radiometric calibration information
dn2ref_ms = os.path.join(output_dir, "dn2reflectance_MS.pix")

dn2ref_pan = os.path.join(output_dir, "dn2reflectance_PAN.pix")
algo.dn2reflectance(fili=ms_link, fileang=sol_ms, filo=dn2ref_ms)
algo.dn2reflectance(fili=pan_link, fileang=sol_pan, filo=dn2ref_pan)

# Create high-resolution color images by fusing black-and-white panchromatic (PAN) and multispectral(MS) color images.
mrafus = os.path.join(output_dir, "mrafusion.pix")

algo.mrafusion(fili=dn2ref_ms, dbic=[1, 2, 3, 4, 5, 6, 7, 8], dbic_ref=[2, 3, 4, 5], fili_pan=dn2ref_pan, dbic_pan=[1],
filo=mrafus)


# Up sample solviewzaz MS to match the PAN resolution and extents
sol_ms_reproj = os.path.join(output_dir, "solviewzaz_MS_reproj.pix")

geo_info = get_geo_info(pan_link)
algo.reproj(fili=sol_ms, dbic=[1, -32], dbsl=[1], filo=sol_ms_reproj, dbsz=[geo_info[3], geo_info[4]], pxsz=geo_info[5],
maxbnds="NO", mapunits=geo_info[0], ulx=str(geo_info[1][0]), uly=str(geo_info[1][1]),
lrx=str(geo_info[2][0]), lry=str(geo_info[2][1]), resample="NEAR")


# Create a spectral classification from optical data by converting it to TOA reflectance and extracting per-pixel features and indices.
spec_fus = os.path.join(output_dir, "speclass.pix")

algo.speclass(fili=mrafus, filo=spec_fus)

# Correct pre-classified radiometric values based on variable incidence angles
toposol = os.path.join(output_dir, "toposolnorm.pix")

algo.toposolnorm(fili=mrafus, filedem=dem_input, filecls=spec_fus, fileang=sol_ms_reproj, filo=toposol)

# Ouput A: Calculate the tasseled-cap transformation for the sensor and bands specified
tassel_topo = os.path.join(output_dir, "topo_tassel.pix")

algo.tassel(fili=toposol, filo=tassel_topo)

# Output B: Compute selected vegetation indices
vegindex_topo = os.path.join(output_dir, "topo_vegindex.pix")

algo.vegindex(fili=toposol, filo=vegindex_topo, index="NDVI")

 

Have more questions? Submit a request

Comments

Powered by Zendesk