python map algebra for NDVI

Discussion created by mwbeckumn on Jul 10, 2011
Latest reply on Aug 12, 2011 by wgivens

I'm attempting to create NDVI layers from a four band raster GRID where layer 1 is red and layer 4 is infrared.  NDVI needs to be calculated as usual, i.e. (red-infrared)/(red+infrared).  Here's my code...

import arcgisscripting,os

gp = arcgisscripting.create()
gp.OverwriteOutput = 1
gp.Workspace = "C:\\Projects\\GIS data\\OBIA nearshore\\temporary\\"

DOW = '01008900'
lake = os.path.join("\\".join(gp.Workspace.split("\\")[:-1]),"buffers all",DOW + "_all_buff.img")
lake_buff = os.path.join("\\".join(gp.Workspace.split("\\")[:-1]),"buffers all",DOW + "_100m.shp")

    print "\nExtracting area of " + DOW + " for processing"
    gp.CopyRaster_management (lake, "g1tmp.img")
    gp.RasterToOtherFormat_conversion("g1tmp.img", gp.Workspace, "GRID")
    gp.ExtractByMask_sa("g1tmp", lake_buff, "g2tmp")

    print "Processing"
    express = str("(" + gp.Workspace + "\\g2tmpc4 - " + gp.Workspace + "\\g2tmpc1) / ( " +
        gp.Workspace + "\\g2tmpc4 + " + gp.Workspace + "\\g2tmpc1)")

    out = str(DOW + "_ndvi.img")

    gp.SingleOutputMapAlgebra_sa(express, out)


    print gp.GetMessages()

print "\ndone"


The call to 'express' returns the following text string, which I thought was correct for map algebra formatting:

'(C:\\Projects\\GIS data\\OBIA nearshore\\temporary\\g2tmpc4 - C:\\Projects\\GIS data\\OBIA nearshore\\temporary\\g2tmpc1) / (C:\\Projects\\GIS data\\OBIA nearshore\\temporary\\g2tmpc4 + C:\\Projects\\GIS data\\OBIA nearshore\\temporary\\g2tmpc1)'

The code runs without errors but I am left with a binary raster layer where the cells take only values zero or negative one.  Any idea what I'm doing wrong?  I'm using Python 2.5.2 with ArcMap9.3.