NDVI Python Script

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[1], tableshape[0], 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):

NDVI

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s