Spectral Indices and Transformation

A Spectral Index is an indicator of a spectral characterstics that can be obtained by combining the reflectance values of two or more bands. On this section, we will take a look at four indices and examine different attributes of a land cover. Later, we will also see how to transform pixel values and detect patterns from a satellite image.

The topics covered on this chapter are:

1. Spectral Idices
    - NDVI, NDBI, NDMI, SAVI
2. Tassled Cap Transformation
    - Brigthness, Greeness, Wetness

Let's start by loading the libraries.

In [1]:
# For reading files
from osgeo import gdal

# For visualization
%matplotlib inline
import matplotlib.pyplot as plt
from rasterio import plot
from rasterio.plot import show

# For Summary and Analysis
import numpy as np

The indicies we will calculate would require a multi-spectral image with at least the RGB,NIR and SWIR bands. Next, we'll simply Load the satellite image we've been using so far and convert raster values to arrays as Float data type. Finally, we'll access bands based on their position and and save them as array.

In [2]:
# Load data
image_data = gdal.Open(r"data/2001_dc_sub.tif")

# Read as Array and Convert to float
image_data_arr = image_data.ReadAsArray()
image_data_f32 = image_data_arr.astype(np.float32)

# Access bands
BLUE = image_data_f32[0]
RED =  image_data_f32[1]
NIR =  image_data_f32[3]
SWIR = image_data_f32[4]

Definitions

NDVI (Normalized DIfference Vegetation Index) is the the most commonly used vegetation index in Remte Sensing. In simple terms, NBDI is a standardized indicator of the greenness of a land.

NDMI (Normalized Difference Moisture Index ) is sensitive to the moisture levels in vegetation. It is used to monitor droughts as well as monitory fuel levels in fire-prone areas. It uses NIR and SWIR bands to create a ratio designed to mitigate illumination and atmospheric effect

NDBI (Normalized Difference Builtup Index) highlights urban areas with higher reflectance in the shortwave-infrared (SWIR) region, compared to the Near Infra-red (NIR) region. Applications include watershed runoff predictions and land-use planning

SAVI (Soil Adjusted Vegetation Index) was developed as a modification of the Normalized Difference Vegetation Index to correct for the influence of soil brightness when vegetative cover is low.

For all the four indices, we will calculate values using the equations for Landsat 4 TM bands. Let's put each equation inside a function.

In [3]:
def NDVI(RED, NIR):
    ndvi = (NIR - RED) / (NIR + RED)
    return ndvi

def NDBI(SWIR, RED):
    ndbi = (SWIR - RED) / (SWIR + RED)
    return ndbi

def NDMI(NIR, SWIR):
    ndmi = (NIR - SWIR) / (NIR + SWIR)
    return ndmi

def SAVI(RED,NIR, L=0.5):
    savi = ((NIR-RED) / (NIR + RED + L)) * (1+L)
    return savi

Next, we will simply specify the bands we accessed earlier.

In [4]:
ndvi_vals = NDVI(RED, NIR)
ndbi_vals = NDBI(SWIR,RED)
ndmi_vals = NDMI(NIR,SWIR)
savi_vals = SAVI(RED, NIR)

Great! Values are calculated. It's time visualize and export.

In [5]:
fig = plt.subplots(2,2,figsize=(20,16))

plt.subplot(221), plt.imshow(ndvi_vals, cmap='RdYlGn'),plt.colorbar(),plt.title("NDVI")
plt.subplot(222), plt.imshow(ndbi_vals, cmap='Greys'),plt.colorbar(),plt.title("NDBI")
plt.subplot(223), plt.imshow(ndmi_vals, cmap='PuBuGn'),plt.colorbar(),plt.title("NDMI")
plt.subplot(224), plt.imshow(savi_vals, cmap='Oranges'),plt.colorbar(),plt.title("SAVI")

plt.show()

The steps of exporting (writing) a processed raster file in python-gdal environment is similar to the one we followed when exporting the PCAs.

Tassled Cap Transformation

The Tasseled-Cap Transformation (TCT) is process of converting the original bands of an image into a new set of bands using weighted sums measures. In other words, the surface reflectance of a satellite image will be converted to new components that can help us better understand the brightness, greeness and wetness of a land.

The function bellow is created by K.Arthur Endsely and we will use later for transforming pixels of each bands with TCT cofficients.

In [6]:
def tasseled_cap_tm(rast, reflectance=True):
    '''
    Applies the Tasseled Cap transformation for TM data. Assumes that the TM
    data are TM reflectance data (i.e., Landsat Surface Reflectance). The
    coefficients for reflectance factor data are taken from Crist (1985) in
    Remote Sensing of Environment 17:302.
    '''
    if reflectance:
        # Reflectance factor coefficients for TM bands 1-5 and 7; they are
        #   entered here in tabular form so they are already transposed with
        #   respect to the form suggested by Kauth and Thomas (1976)
        r = np.array([
            ( 0.2043, 0.4158, 0.5524, 0.5741, 0.3124, 0.2303),
            (-0.1603,-0.2819,-0.4934, 0.7940,-0.0002,-0.1446),
            ( 0.0315, 0.2021, 0.3102, 0.1594,-0.6806,-0.6109),
            (-0.2117,-0.0284, 0.1302,-0.1007, 0.6529,-0.7078),
            (-0.8669,-0.1835, 0.3856, 0.0408,-0.1132, 0.2272),
            ( 0.3677,-0.8200, 0.4354, 0.0518,-0.0066,-0.0104)
        ], dtype=np.float32)

    else:
        raise NotImplemented('Only support for Landsat TM reflectance has been implemented')

    shp = rast.shape

    # Can accept either a gdal.Dataset or numpy.array instance
    if not isinstance(rast, np.ndarray):
        x = rast.ReadAsArray().reshape(shp[0], shp[1]*shp[2])

    else:
        x = rast.reshape(shp[0], shp[1]*shp[2])

    return np.dot(r, x).reshape(shp)

Let's go ahead and pass the function to the image in array form and to visualize the transformed bands, we will subset each components

In [7]:
tassled_cap = tasseled_cap_tm(image_data_arr)
ts_1 = tassled_cap[0] # Brightness
ts_2 = tassled_cap[1] # Greeness
ts_3 = tassled_cap[2] # Wetness

The TCT feature space below shows us the relationship between Brightness and Greeness.

In [8]:
fig, ax = plt.subplots(figsize=(10,8))
plt.scatter(ts_1, ts_2, alpha = 0.3)
plt.title("Landsat-4 Tassled Cap Components")
plt.xlabel('Tassled Cap Brightness(TC 1)')
plt.ylabel('Tassled Cap Greeness (TC 2)')
plt.show()

Let's also plot the three transformed components.

In [9]:
fig,ax = plt.subplots(1,3,figsize=(20,6))

plt.subplot(131), plt.imshow(ts_1, cmap='Reds'), plt.colorbar(),plt.title("BRIGHTNESS")
plt.subplot(132), plt.imshow(ts_2, cmap='Greens'),plt.colorbar(),plt.title("GREENESS")
plt.subplot(133), plt.imshow(ts_3, cmap='Blues'),plt.colorbar(),plt.title("WETNESS")

plt.show()

The TCT components resembles what we generally know about the area. The Brightness is relatively higher in urban or builtup areas and the Greeness is similar to NDVI values generated ealier. The Wetness component also clearly identified the water bodies.

Summary

On this document, we covered spectral indices and Tassled Cap Transformation. These indicators serve as instruments to evaluate conditions mostly in relateion to vegetation cover. The next and final part will cover image classification.