Skip to main content

Lab 3A - Spectral Indices

Overview

In this lab, we will work with the spectral characteristics in our data to visualize and extract insights that go beyond basic visual interpretation. We will work with the different spectral bands offered by Landsat 8 to find unique patterns that can help us solve problems and conduct analysis. By the end of this lab, you should be able to understand how to build and visualize existing indices, as well as construct your own, identify how different indices can help your use case, and understand the mechanism behind how they work.

Spectral Indices

Spectroscopy is the study of how radiation is absorbed, reflected and emitted by different materials. While this discipline has its origins in chemistry and physics, we can utilize the same techniques to identify different land cover types from satellite data. In the chart below, land cover types have unique spectral characteristics. Snow has a major peak at lower wavelengths and is near zero above 1.5 micrometer, whereas soil has very low reflectance at lower levels of wavelength but relatively strong and steady reflectance after ~0.75 micrometers. Spectral indices are built to leverage these unique characteristics and isolate specific types of land cover.

Land covers are separable at different wavelengths. Vegetation curves (green) have high reflectance in the NIR range, where radiant energy is scattered by cell walls (Bowker, 1985) and low reflectance in the red range, where radiant energy is absorbed by chlorophyll. We can leverage this information to build indices that help us differentiate vegetation from urban areas. In the next few sections, we will cover several of the most important indices in use.

Land Cover Reflectance

Important Indices

Normalized Difference Vegetation Index (NDVI)

The Normalized Difference Vegetation Index (NDVI) has a long history in remote sensing, and is one of the most widely used measures. The typical formulation is:

NDVI=(NIRred)/(NIR+red)\text{NDVI} = (\text{NIR} - \text{red}) / (\text{NIR} + \text{red})

Where NIR refers to the near infrared band and red refers to the red peak in the visible spectrum.

Because NDVI is a popular and well-known index, we can use the built-in functionality within Earth Engine normalizedDifference()to calculate NDVI. You can follow the steps below to build your own indices.

First, build a baseline true color image around our region of interest, Niger. We will work with the Landsat 8 Collection 1 Tier 1 TOA Reflectance data from 2015, sort by cloud cover and extract the first image.

#!pip install geemap
import ee, geemap, pprint
#ee.Authenticate()
def build_map(lat, lon, zoom, vizParams, image, name):
map = geemap.Map(center = [lat, lon], zoom = zoom)
map.addLayer(image, vizParams, name)
return map
ee.Initialize()

var point = ee.Geometry.Point([-80.42, 37.22]);
var landsat = ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA")
var image = ee.Image(landsat
.filterBounds(point)
.filterDate('2015-06-01', '2015-09-01')
.sort('CLOUD_COVER')
.first());
var trueColor = {bands: ['B4', 'B3', 'B2'],
min: 0, max: 0.3};
Map.centerObject(point, 12);
Map.addLayer(image, trueColor, 'image');

image

Now that we have the true color baseline image, we can build the NDVI index and visualize it. For visualization, we are creating a custom palette, where low values trend towards white and high values trend towards green.

var point = ee.Geometry.Point([-80.42, 37.22]);
var landsat = ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA")
var image = ee.Image(landsat
.filterBounds(point)
.filterDate('2015-06-01', '2015-09-01')
.sort('CLOUD_COVER')
.first());
var ndvi = image.normalizedDifference(['B5', 'B4']);
var vegPalette = ['white', 'green'];
Map.centerObject(point, 12);
Map.addLayer(ndvi, {min: -1, max: 1,
palette: vegPalette}, 'NDVI');

image

Question 1A: What are some of the sample pixel values of the NDVI in the below categories. Indicate which parts of the images you used and how you determined what each of their values were.

  1. Vegetation
  2. Urban features
  3. Bare earth
  4. Water

Enhanced Vegetation Index (EVI)

The Enhanced Vegetation Index (EVI) is designed to minimize saturation and background effects in NDVI (Huete, 2002).

EVI=2.5(NIRred)/(NIR+6red7.5blue+1)\text{EVI} = 2.5 * (\text{NIR} - \text{red}) / (\text{NIR} + 6 * \text{red} - 7.5 * \text{blue} + 1)

Since it is not a normalized difference index, we need to build a unique expression and then identify all of the different segments. Programmatically, bands are specifically referenced with the help of an object that is passed as the second argument to image.expression() (everything within the curly brackets).

var point = ee.Geometry.Point([-80.42, 37.22]);
var landsat = ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA")
var image = ee.Image(landsat
.filterBounds(point)
.filterDate('2015-06-01', '2015-09-01')
.sort('CLOUD_COVER')
.first());
// Build the expression
var exp = '2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))';
var evi = image.expression( exp,
{'NIR': image.select('B5'),
'RED': image.select('B4'),
'BLUE': image.select('B2')
});
var vegPalette = ['white', 'green'];
Map.centerObject(point, 12);
Map.addLayer(evi,
{min: -1, max: 1, palette: vegPalette},
'EVI');

image

Question 1B: Compare EVI to NDVI across those same land use categories as in the previous question. What do you observe -- how are the images and values similar or different across the two indices?

Normalized Difference Water Index (NDWI)

The Normalized Difference Water Index (NDWI) was developed by Gao (1996) as an index to identify the water content within vegetation. SWIR stands for short-wave infrared, which is the Landsat band 6. This is not an exact implementation of NDWI, according to the OLI spectral response, since OLI does not have a band in the right position (1.26 𝛍m) - but for our purposes, this is an approximation that does an acceptable job of identifying water content.

NDWI=(NIRSWIR))/(NIR+SWIR)\text{NDWI} = (\text{NIR} - \text{SWIR})) / (\text{NIR} + \text{SWIR})

var point = ee.Geometry.Point([-80.42, 37.22]);
var landsat = ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA")
var image = ee.Image(landsat
.filterBounds(point)
.filterDate('2015-06-01', '2015-09-01')
.sort('CLOUD_COVER')
.first());
var ndwi = image.normalizedDifference(['B5', 'B6']);
var waterPalette = ['white', 'blue'];
Map.centerObject(point, 12);
Map.addLayer(ndwi,
{min: -1, max: 1, palette: waterPalette},
'NDWI');


image

Normalized Difference Water Body Index (NDWBI)

The fact that two different NDWI indices were independently invented in 1996 complicates things. While the NDWI looks at water content within vegetation, the NDWBI is built to identify bodies of water (rivers, lakes, oceans). To distinguish, define the Normalized Difference Water Body Index (NDWBI) as the index described in McFeeters (1996):

NDWBI=(greenNIR)/(green+NIR)\text{NDWBI} = (\text{green} - \text{NIR}) / (\text{green} + \text{NIR})

As previously, implement NDWBI with normalizedDifference() and display the result.

var point = ee.Geometry.Point([-80.42, 37.22]);
var landsat = ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA")
var image = ee.Image(landsat
.filterBounds(point)
.filterDate('2015-06-01', '2015-09-01')
.sort('CLOUD_COVER')
.first());
var waterPalette = ['white', 'blue'];
var ndwbi = image.normalizedDifference(['B3', 'B5']);
Map.addLayer(ndwbi,
{min: -1,
max: 0.5,
palette: waterPalette},
'NDWBI');


image

You can combine the code blocks to compare the actual values at different pixel locations. Use inspector to test out different land areas.

Question 2: Compare NDWI and NDWBI. What do you observe? What do each of the indices try to focus on?

Normalized Difference Bare Index (NDBI)

The Normalized Difference Bare Index (NDBI) was developed by Zha, 2003) to aid in the differentiation of urban areas by using a combination of the shortwave and near infrared.

NDBI=(SWIRNIR)/(SWIR+NIR)\text{NDBI} = (\text{SWIR} - \text{NIR}) / (\text{SWIR} + \text{NIR})

Note that NDBI is the negative of NDWI. Compute NDBI and display with a suitable palette. (Check this reference to demystify the palette reversal)

var point = ee.Geometry.Point([-80.42, 37.22]);
var landsat = ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA")
var image = ee.Image(landsat
.filterBounds(point)
.filterDate('2015-06-01', '2015-09-01')
.sort('CLOUD_COVER')
.first());
var waterPalette = ['white', 'blue'];
var barePalette = waterPalette.slice().reverse();
var ndbi = image.normalizedDifference(['B6', 'B5']);
Map.addLayer(ndbi, {min: -1, max: 0.5, palette: barePalette}, 'NDBI');


image

Burned Area Index (BAI)

The Burned Area Index (BAI) was developed by Chuvieco et al. (2002) to assist in the delineation of burn scars and assessment of burn severity. It is based on maximizing the spectral characteristics of charcoal reflectance. To examine burn indices, load an image from 2013 showing the Rim fire in the Sierra Nevadas. We'll start by creating a true image of the area to see how well this index highlights the presence of wildfire.

var point = ee.Geometry.Point(-120.083, 37.850);
var landsat = ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA")
var image = ee.Image(landsat
.filterBounds(point)
.filterDate('2013-08-17', '2013-09-27')
.sort('CLOUD_COVER')
.first());
var trueColor = {bands: ['B4', 'B3', 'B2'],
min: 0, max: 0.3};
Map.centerObject(point, 12);
Map.addLayer(image, trueColor, 'image');

image

Closely examine the true color display of this image. Can you spot where the fire occurred? If difficult, let's look at the burn index.

var point = ee.Geometry.Point(-120.083, 37.850);
var landsat = ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA")
var image = ee.Image(landsat
.filterBounds(point)
.filterDate('2013-08-17', '2013-09-27')
.sort('CLOUD_COVER')
.first());
var trueColor = {bands: ['B4', 'B3', 'B2'],
min: 0, max: 0.3};
Map.centerObject(point, 12);
Map.addLayer(image, trueColor, 'image');

// Build Burn Index expression
var exp = '1.0 / ((0.1 - RED)**2 + (0.06 - NIR)**2)';
var bai = image.expression(exp,
{'NIR': image.select('B5'),
'RED': image.select('B4') });

var burnPalette = ['green', 'blue', 'yellow', 'red'];
Map.addLayer(bai, {min: 0, max: 400, palette: burnPalette}, 'BAI');

image

The charcoal burn area is now very evident. Being that Landsat has historical data and a wide array of sensors, this can be a powerful way to understand natural phenomena.

Question 3: Compare NDBI and the BAI displayed results -- what do you observe?

Normalized Burn Ratio Thermal (NBRT)

The Normalized Burn Ratio Thermal (NBRT) was developed based on the idea that burned land has low NIR reflectance (less vegetation), high SWIR reflectance (think ash), and high brightness temperature (Holden et al. 2005). Unlike the other indices, a lower NBRT means a higher likelihood of recent burn (for visualization, reverse the scale). This index can be used to diagnose the severity of wildfires (see van Wagtendonk et al. 2004).

var point = ee.Geometry.Point(-120.083, 37.850);
var landsat = ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA")
var image = ee.Image(landsat
.filterBounds(point)
.filterDate('2013-08-17', '2013-09-27')
.sort('CLOUD_COVER')
.first());
var trueColor = {bands: ['B4', 'B3', 'B2'],
min: 0, max: 0.3};
Map.centerObject(point, 12);
Map.addLayer(image, trueColor, 'image');
// Build Burn Index expression
var exp = '(NIR - 0.0001 * SWIR * Temp) / (NIR + 0.0001 * SWIR * Temp)'
var nbrt = image.expression(exp,
{'NIR': image.select('B5'),
'SWIR': image.select('B7'),
'Temp': image.select('B11')
});
var burnPalette = ['green', 'blue', 'yellow', 'red'];
Map.addLayer(nbrt, {min: 1, max: 0.9, palette: burnPalette}, 'NBRT');

image

Normalized Difference Snow Index (NDSI)

The Normalized Difference Snow Index (NDSI) was designed to estimate the amount of a pixel covered in snow (Riggs et al. 1994).

NDSI=(greenSWIR)/(green+SWIR)\text{NDSI} = (\text{green} - \text{SWIR}) /(\text{green} + \text{SWIR})

Let's look at Aspen, Colorado and use Landsat 8 data in the winter. You can use the layer manager to turn on and off the snow layer to compare results with the true color image. How does it compare? Reference the spectral reflectance chart at the beginning of the lab and look at the profile for snow. You will see that it has a distinct profile. In the image below, it does a very good job of matching the true color image - valleys and roads that do not have snow on them are accurately shown.

var point = ee.Geometry.Point([-106.81, 39.19]);
var landsat = ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA")
var image = ee.Image(landsat
.filterBounds(point)
.filterDate('2013-11-17', '2014-03-27')
.sort('CLOUD_COVER')
.first());
var trueColor = {bands: ['B4', 'B3', 'B2'],
min: 0, max: 0.3};
Map.centerObject(point, 12);
Map.addLayer(image, trueColor, 'image');
var ndsi = image.normalizedDifference(['B3', 'B6']);
var snowPalette = ['red', 'green', 'blue', 'white'];
Map.addLayer(ndsi,
{min: -0.5, max: 0.5, palette: snowPalette},
'NDSI');

image