After the “What is NDVI?” post I decided to create a Python Script for NDVI. This script works with a *.tif NIR band and a *.tif RED band and gives as output a new NDVI image in *.tif format. To run this script, it’s important the name of the RED band to be “red.tif“, the name of the NIR band “nir.tif“ and all of them to be in the same folder (Python Script, nir.tif and red.tif). Also you could change the names of the images from lines 21 and 26. GDAL and Python must be installed on your computer. This script is under the terms of the GNU GENERAL PUBLIC LICENSE.
# NDVI Python Script # # GNU GENERAL PUBLIC LICENSE # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or #(at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # Created by Alexandros Falagas. # from osgeo import gdal # this allows GDAL to throw Python Exceptions gdal.UseExceptions() from gdalconst import * import numpy as np import sys red=gdal.Open('red.tif', GA_ReadOnly) if red is None: print 'Could not open file' sys.exit(1) r=np.array(red.GetRasterBand(1).ReadAsArray(), dtype=float) nir=gdal.Open('nir.tif', GA_ReadOnly) if nir is None: print 'Could not open file' sys.exit(1) n=np.array(nir.GetRasterBand(1).ReadAsArray(), dtype=float) geotr=red.GetGeoTransform() proj=red.GetProjection() tableshape=r.shape np.seterr(divide='ignore', invalid='ignore') #Ignore the divided by zero or Nan appears ndvi=(n-r)/(n+r) # The NDVI formula driver=gdal.GetDriverByName('GTiff') dst_ds=driver.Create('ndvi.tif', tableshape, tableshape, 1, gdal.GDT_Float32) dst_ds.SetGeoTransform(geotr) dst_ds.SetProjection(proj) dst_ds.GetRasterBand(1).WriteArray(ndvi) dst_ds=None # save, close print 'The NDVI image is saved.'
Here is the result after using the script in a Landsat 8 image (Edited in QGIS):