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

Out:

```<AxesSubplot:title={'center':'Landsat 8 - Normalized Difference Vegetation Index (NDVI)'}>
```

## 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
)
np.unique(ndvi_landsat_class)
```

Out:

```masked_array(data=[1, 2, 3, 4, 5, --],
mask=[False, False, False, False, False,  True],
fill_value=999999)
```

## 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]

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

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

Gallery generated by Sphinx-Gallery