Calculate and Classify Normalized Difference Results with EarthPy

Learn how to calculate and classify normalized difference indices in Python using EarthPy. This example shows how to calculate and classify the normalized difference vegetation index (NDVI) using Landsat 8 data.

Calculating Normalized Difference in Python Using EarthPy

Note

The examples below will show you how to use the normalized_diff() function to calculate the normalized difference vegetation index (NDVI), a commonly used remotely sensed index for quantifying vegetation health.

The example below walks you through a typical workflow for calculating the normalized difference vegetation index (NDVI) using Landsat 8 data with EarthPy. NDVI provides a measure of healthy vegetation and ranges in value from -1 to 1. Values closer to 1 represent healthy, green vegetation. NDVI can be calculated from Landsat 8 data using band 4 (red) and band 5 (near-infrared).

First, you will create a stack of bands using Landsat 8 data and then calculate NDVI using the normalized_diff() function. Then, you will plot the NDVI results using a colorbar legend with continuous values. Last, you will classify the NDVI results using threshold values and plot the classified data with a categorical legend.

Import Packages

To begin, import the needed packages. You will use a combination of several EarthPy modules including spatial and plot.

import os
from glob import glob
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
import earthpy as et
import earthpy.spatial as es
import earthpy.plot as ep

# Get data and set your home working directory
data = et.data.get_data("vignette-landsat")

Import Example Data

To get started, make sure your directory is set. Then, create a stack from all of the Landsat .tif files (one per band). The nodata value for Landsat 8 is -9999 so you can use the nodata= parameter when you call the stack() function.

os.chdir(os.path.join(et.io.HOME, "earth-analytics"))

# Stack the Landsat 8 bands
# This creates a numpy array with each "layer" representing a single band
landsat_path = glob(
    "data/vignette-landsat/LC08_L1TP_034032_20160621_20170221_01_T1_sr_band*_crop.tif"
)
landsat_path.sort()
arr_st, meta = es.stack(landsat_path, nodata=-9999)

Calculate Normalized Difference Vegetation Index (NDVI)

You can calculate NDVI for your dataset using the normalized_diff function from the earthpy.spatial module. Math will be calculated (b1-b2) / (b1 + b2).

# Landsat 8 red band is band 4 at [3]
# Landsat 8 near-infrared band is band 5 at [4]
ndvi = es.normalized_diff(arr_st[4], arr_st[3])

Plot NDVI With Colorbar Legend of Continuous Values

You can plot NDVI with a colorbar legend of continuous values using the plot_bands function from the earthpy.plot module.

titles = ["Landsat 8 - Normalized Difference Vegetation Index (NDVI)"]

# Turn off bytescale scaling due to float values for NDVI
ep.plot_bands(
    ndvi, cmap="RdYlGn", cols=1, title=titles, scale=False, vmin=-1, vmax=1
)
../_images/sphx_glr_plot_calculate_classify_ndvi_001.png

Classify NDVI

Next, you can categorize (or classify) the NDVI results into useful classes. Values under 0 will be classified together as no vegetation. Additional classes will be created for bare area and low, moderate, and high vegetation areas.

# Create classes and apply to NDVI results
ndvi_class_bins = [-np.inf, 0, 0.1, 0.25, 0.4, np.inf]
ndvi_landsat_class = np.digitize(ndvi, ndvi_class_bins)

# Apply the nodata mask to the newly classified NDVI data
ndvi_landsat_class = np.ma.masked_where(np.ma.getmask(ndvi), ndvi_landsat_class)
np.unique(ndvi_landsat_class)

Plot Classified NDVI With Categorical Legend - EarthPy Draw_Legend()

You can plot the classified NDVI with a categorical legend using the draw_legend() function from the earthpy.plot module.

# Define color map
nbr_colors = ["gray", "y", "yellowgreen", "g", "darkgreen"]
nbr_cmap = ListedColormap(nbr_colors)

# Define class names
ndvi_cat_names = [
    "No Vegetation",
    "Bare Area",
    "Low Vegetation",
    "Moderate Vegetation",
    "High Vegetation",
]

# Get list of classes
classes = np.unique(ndvi_landsat_class)
classes = classes.tolist()
# The mask returns a value of none in the classes. remove that
classes = classes[0:5]

# Plot your data
fig, ax = plt.subplots(figsize=(12, 12))
im = ax.imshow(ndvi_landsat_class, cmap=nbr_cmap)

ep.draw_legend(im_ax=im, classes=classes, titles=ndvi_cat_names)
ax.set_title(
    "Landsat 8 - Normalized Difference Vegetation Index (NDVI) Classes",
    fontsize=14,
)
ax.set_axis_off()

# Auto adjust subplot to fit figure size
plt.tight_layout()
../_images/sphx_glr_plot_calculate_classify_ndvi_002.png

Total running time of the script: ( 0 minutes 7.533 seconds)

Gallery generated by Sphinx-Gallery