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.
#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
#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)
You can download the input data and scripts used for the tutorial on this link.