NDVI calculation in Python

1   Introduction

Implement the normalized difference vegetation index (NDVI).

2   Python code

from osgeo import gdal
import matplotlib.pyplot as plt

ds = gdal.Open("p:/courses/python/m_3008101_ne_17_060_20191110.tif")

red_band = ds.GetRasterBand(1)
nir_band = ds.GetRasterBand(4)

red_arr = red_band.ReadAsArray().astype("float")
nir_arr = nir_band.ReadAsArray().astype("float")

# NDVI calculation is very easy because numpy does the heavy lifting for us!
ndvi_arr = (nir_arr - red_arr) / (nir_arr + red_arr)

driver = gdal.GetDriverByName("GTiff")
out_ds = driver.Create("p:/courses/python/ndvi.tif", ndvi_arr.shape[1], ndvi_arr.shape[0], 1, gdal.GDT_Float32)
out_ds.SetProjection(ds.GetProjection())
out_ds.SetGeoTransform(ds.GetGeoTransform())

out_band = out_ds.GetRasterBand(1)
out_band.WriteArray(ndvi_arr)
out_band.FlushCache()
out_band.ComputeStatistics(False)

plt.imshow(ndvi_arr)
plt.show()