plot_calculate_classify_ndvi.py

"""
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, vmin=-1, vmax=1)

###############################################################################
# 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()