Soil-Adjasted Vegetation Index (SAVI) Python Script

Empirically derived NDVI products have been shown to be unstable, varying with soil color, soil moisture, and saturation effects from high density vegetation. SAVI is a transformation technique that minimizes soil brightness influences from spectral vegetation indices involving red and near-infrared wavelengths.

The index is given as:

  {\displaystyle {\mbox{SAVI}}={\frac {({\mbox{1}}+{\mbox{L}})({\mbox{NIR}}-{\mbox{Red}})}{({\mbox{NIR}}+{\mbox{Red}}+{\mbox{L}})}}}

where L is a canopy background adjustment factor. An L value of 0.5 in reflectance space was found to minimize soil brightness variations and eliminate the need for additional calibration for different soils. The transformation was found to nearly eliminate soil-induced variations in vegetation indices.

Click to the video bellow to see how this script works on a Linux Ubuntu 16.04 System.

Python, GDAL, Numpy and Time must be installed on your computer. The script bellow is under the terms of the 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 2 of the License, or     *
#* (at your option) any later version.                                   *
#*************************************************************************
#copyright: (C) 2016 by Falagas Alexandros

#SAVI VEGETATION INDEX

def SAVI (n, r, L):
np.seterr(divide='ignore', invalid='ignore') #Ignore the divided by zero or Nan appears
savi=((n-r)/(n+r+L))*(1+L)
return savi
import os
import sys
import gdal
# this allows GDAL to throw Python Exceptions
gdal.UseExceptions()
from gdalconst import *
import numpy as np
import time
print '***************SAVI Index calculation ***************'
cwd=os.getcwd()
print 'Your Current Working Directory:', cwd
for file in os.listdir(cwd): # list all files from the directory.
print file
form=str(raw_input('Press the image format (for ex. .TIF, .JPG etc):'))
for file in os.listdir(cwd):
if file.endswith(form):
print file # list all files ends with the specific image format.
redim=str(raw_input('Press the code name of the RED image:'))
nirim=str(raw_input('Press the code name of the NIR image:'))
saviim='SAVI'+form
L=float(raw_input('Press L soil factor (0-Very high vegetation regions, 1-No vegetation regions):'))
if L>1 and L<0:
print 'L factor out of limits.'
sys.exit(1)
startTime = time.time()
red=gdal.Open(redim, 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(nirim, 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
savi=SAVI (n, r, L)
driver=gdal.GetDriverByName('GTiff')
dst_ds=driver.Create(saviim, tableshape[1], tableshape[0], 1, gdal.GDT_Float32)
dst_ds.SetGeoTransform(geotr)
dst_ds.SetProjection(proj)
dst_ds.GetRasterBand(1).WriteArray(savi)
dst_ds=None # save, close
print 'The SAVI image '+saviim+' is saved.'
elapsedTime = time.time() - startTime
print 'The SAVI image is saved. The algorithm is completed after {} sec.'.format(round(elapsedTime, 2))
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