NDVI calculation from Landsat8 images with Python 3 and Rasterio - Tutorial

The NDVI is a vegetation index widely used for environmental impact assessment, agricultural evaluation, and land use change metrics. The procedure for the calculation of the NDVI is simple and straighforward in softwares of Geographical Information Systems (GIS) as QGIS. However the efficiency is relevant only we deal with one image, but what would happend if we have to analyze a series of images, or if we have large images with limited computational resources, or if we have to apply some custom filters and preprocessing to the dataset, then we have to look for another methods to calculate the NDVI on a faster and more efficient way according to the computational resources available.

Satellite images are georasters, these images are a regular array of columns and rows (a matrix per band) with a georeferenciation. Python is a programming and data analysis language very versatile for the matrix algebra with the Numpy library, however there was no efective and simple way to process a georaster until the development of the Rasterio package.

Rasterio is a library to open, write, explore and analyze georasters in Python. The library uses GeoTIFF images along with other formats and is capable to work with satellite images, digital elevation models, and drone generated imagery.

This tutorial show the complete procedure to analyse the NDVI from a Landsat 8 image with Python 3 and Rasterio. The scripting and representation was performed on a interactive enviroment called Jupyter Notebook, finally the result georaster was opened in QGIS and compared with some background images.

Tutorial

Python code

#import required libraries
import rasterio
from rasterio import plot
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

import os
os.listdir('../Landsat8/')

#import bands as separate 1 band raster
band4 = rasterio.open('../Landsat8/LC08_L1TP_042035_20180603_20180615_01_T1_B4_clip.tif') #red
band5 = rasterio.open('../Landsat8/LC08_L1TP_042035_20180603_20180615_01_T1_B5_clip.tif') #nir

#number of raster rows
band4.height

#number of raster columns
band4.width

#plot band 
plot.show(band4)

#type of raster byte
band4.dtypes[0]

#raster sytem of reference
band4.crs

#raster transform parameters
band4.transform

#raster values as matrix array
band4.read(1)

#multiple band representation
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))
plot.show(band4, ax=ax1, cmap='Blues') #red
plot.show(band5, ax=ax2, cmap='Blues') #nir
fig.tight_layout()

#generate nir and red objects as arrays in float64 format
red = band4.read(1).astype('float64')
nir = band5.read(1).astype('float64')

nir

#ndvi calculation, empty cells or nodata cells are reported as 0
ndvi=np.where(
    (nir+red)==0., 
    0, 
    (nir-red)/(nir+red))
ndvi[:5,:5]

#export ndvi image
ndviImage = rasterio.open('../Output/ndviImage.tiff','w',driver='Gtiff',
                          width=band4.width, 
                          height = band4.height, 
                          count=1, crs=band4.crs, 
                          transform=band4.transform, 
                          dtype='float64')
ndviImage.write(ndvi,1)
ndviImage.close()

#plot ndvi
ndvi = rasterio.open('../Output/ndviImage.tiff')
fig = plt.figure(figsize=(18,12))
plot.show(ndvi)

Input data

You can download the input data and scripts used for the tutorial on this link.