Using Python to calculate NDVI with multiband imagery


The following is a helpful section of python code to assist with calculating NDVI (Normalised Difference Vegetation Index) using the arcpy site package.  The Image Analysis Window allows you to see this visually, but if you want the values in the raster itself, then you will need to use this code.

Thanks go to the Esri team for explaining this new notation in the Spatial Analyst module of arcpy (arcpy.sa).  The arcpy.sa module calculates rasters more effectively than the Raster Calculator at 9.x, as it uses the number crunching power of Python.  Also, dot notation for multiple bands no longer works in the Raster Calculator at 10, so this code is a optimised replacement for this loss in functionality.

In this example, Band 3 is Red, Band 5 is Near IR.

To run this script, open IDLE or another IDE (like PythonWin, PyChecker, or WingIDE), and create a new script.  Copy and paste the below script and paste into that new script.  Then you need to read the green comments, and replace the CAPITALISED sections.  It will create a few outputs, like NDVI and Red band rasters.  These can be deleted once you are happy with the result.

____________________copy below into IDE__________________

#import modules arcpy and string

import arcpy, string

#import the environment settings and spatial analyst extension

from arcpy import env
from arcpy.sa import *

#Check out the Spatial Analyst extension (must be available to check out!)

arcpy.CheckOutExtension(“spatial“)

#Add your workspace here

env.workspace = r”C:YOURWORKSPACE

#set the input raster here (will need to be inside the workspace shown above)

input = r”YOURRASTERFILE.tif

#Create raster directory and specify inputs. It references the input you specified above
result = “FINALNDVIOUTPUTNAME.tif
NIR = input + “\Band_5
Red = input + “\Band_3

#Define outputs – Note:  These will need to be deleted if you need to run this script again
NIR_out = “NIR.tif
Red_out = “Red.tif

#Copy raster and map algebra – makes new rasters of Red and NearIR
arcpy.CopyRaster_management(NIR, NIR_out)
print “Copied NIR band as raster
arcpy.CopyRaster_management(Red, Red_out)
print “Copied Red band as raster

#Create Numerator and Denominator rasters as variables and NDVI output (note that arcpy.sa.Float returns a floating point raster)

Num = arcpy.sa.Float(Raster(NIR_out) – Raster(Red_out))
Denom = arcpy.sa.Float(Raster(NIR_out) + Raster(Red_out))
NIR_eq = arcpy.sa.Divide(Num, Denom)
print “Dividing

#Saving output to result output you specified above
NIR_eq.save(result)
print “Successful

____________________________________________

If you don’t think you are using Spatial Analyst to it’s full extent, perhaps have a look at the Spatial Analyst training course?

Happy NDVI-ing.

Kath S

14 thoughts on “Using Python to calculate NDVI with multiband imagery

  1. Bayaraa

    Hi Kath

    I am using python2.4. but when i run my code, there is error about importing arcpy. do you have some advice. thanks.

    Reply
  2. ovidiu

    Hi,
    What kind of raster should be put here ? Do we need something else beside the 2 bands ?

    #set the input raster here (will need to be inside the workspace shown above)
    input = r”YOURRASTERFILE.tif“

    Reply
  3. ksund Post author

    The only imagery you need is a multi-band image of any format (tif, gif, jpeg, GRID, geodatabase raster format). If you preview the raster in ArcGIS, there should be multiple bands when you look at the Source tab of the raster. If there are multiple bands built into the raster, then you can use the above code. If the bands are in separate rasters, then you only need to specify the two rasters as the NIR_in and RED_in. For example:

    NIR_in = ‘C:\data\yourNIFimage.tif’
    RED_in = ‘C:\data\yourREDimage.tif’

    Then the rest of the script will work as you specified.

    Hope that helps.

    Reply
  4. Luke T

    Thanks for this Kath. I added arcpy.ListRasters to my script to loop through all my TIF tiles to make the NDVI tiles, lists are great.

    Reply
      1. Luke T

        Hi Jon,
        Sorry I just received your message. I think this is probably the script you’re after. I don’t know how wordpress will format it but all you have to know about the indentation is everything below the for loop is indented once.

        #script copied from here https://esriaustraliatechblog.wordpress.com/2012/03/14/using-python-to-calculate-ndvi-with-multiband-imagery/
        #modified by Luke T***** on 2/8/2012
        #WARNING: this script deletes input rasters

        #import modules arcpy and string
        import arcpy
        import string

        #import the environment settings and spatial analyst extension
        from arcpy import env
        from arcpy.sa import *

        #check out the spatial analyst extension (must be available to check out)
        arcpy.CheckOutExtension(“spatial”)

        #add your workspace here
        env.workspace = “C:/GIS/scratch_aerial”

        #get a list of TIFFs in the workspace(seeing as their TIFFs the workspace is a folder)
        rasterlist = arcpy.ListRasters(“*”, “TIF”)
        for inputraster in rasterlist:

        #Specify output directory and specify inputs. It references the input you specified in the for loop
        result = “C:/GIS/scratch_aerial/outputs/” + “out” + inputraster
        NIR = inputraster + “/Band_4”
        Red = inputraster + “/Band_1”

        #Define outputs – Note: Tese will need to be deleted if you need to run this script again
        NIR_out = “NIR.tif”
        Red_out = “Red.tif”

        #Copy raster and map algebra – makes new raster of Red and NearIR
        arcpy.CopyRaster_management(NIR, NIR_out)
        print “Copied NIR band as raster”
        arcpy.CopyRaster_management (Red, Red_out)
        print “Copied Red band as raster”

        #Create Numerator and Denominator rasters as variables and NDVI output
        #note that arcpy.sa.Fload returns a floating point raster)
        Num = arcpy.sa.Float(Raster(NIR_out) – Raster(Red_out))
        Denom = arcpy.sa.Float(Raster(NIR_out) + Raster(Red_out))
        NIR_eq = arcpy.sa.Divide(Num, Denom)
        print “Dividing”

        #Saving output to result output you specified above
        NIR_eq.save(result)
        print “Successful”

        #Delete the intermediate TIFFs
        arcpy.Delete_management(“NIR.tif”, “TIFF”)
        print arcpy.GetMessages()
        arcpy.Delete_management(“Red.tif”, “TIFF”)
        print arcpy.GetMessages()
        #Delete the input TIFFs since the space may be required depending on the number and size of tiles
        arcpy.Delete_management(inputraster)
        print arcpy.GetMessages()

  5. Cahya

    i have problem in saving the image…

    this is the error message:

    Traceback (most recent call last):
    File “d:\***Spasial\NDVI\ndvi.tbx#NDVI.py”, line 56, in
    RuntimeError: ERROR 010240: Could not save raster dataset to C:\Users\***\Documents\ArcGIS\Default.gdb\float_ras with output format FGDBR.

    Failed to execute (NDVI).

    any conclusion?

    i have succed saving the image once…but now failed

    THX before

    Reply
    1. ksund Post author

      Hi Cahya,

      You may need to try savind the output to a different geodatabase. If your default geodatabase is stored in your windows profile, it might be limiting the ability to save there (if you have a max size on your Windows profile). Try building a c:/gistemp folder, and creating a gistemp.gdb there. Then set the env.workspace = “C://gistemp/gistemp.gdb”. See if you can save it to the that geodatabase. If that doesn’t work, then you may need to look at permissions, or increasing the cell size. Also worth setting your env.scratchworkspace to the same temp geodatabase (again, it could be set to your windows profile.

      Hope that helps.

      Kath S

      Reply
  6. DM

    Awesome! Yeah! Yahoo!
    Thanks so much! This was a crux I had in my workflow that really needed to be solved. (vs using Arcmap, NDVI, then export – which usually crashed).

    Reply

Got something to say?

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