python map algebra for NDVI

Question

Hi,

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.CheckOutExtension("Spatial")
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")

try:
    print "
Extracting 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)

except:

print gp.GetMessages()

print "
done"

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.

Thanks,

Marcus

WadeGivens · Answer

How would this work if you had a multi-band tiff?

Thanks,
Wade

MarcusBeck · Answer

Ah yes that was the problem.  I added a command (as a list comprehension) to change the integer grids to float and *greatly* simplified the map algebra expression.  Everything works as it should.  Thanks Dan!
import arcgisscripting,os

gp = arcgisscripting.create()
gp.OverwriteOutput = 1
gp.CheckOutExtension("Spatial")
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")

try:
    print "
Extracting 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 "Converting to float"
    [gp.float_sa("g2tmp" + x, "g2tmp" + x + "_flt") for x in ("c1","c4")]

print "Processing"
    express = "((g2tmpc4_flt - g2tmpc1_flt) / (g2tmpc4_flt + g2tmpc1_flt))"
    out = DOW + "_ndvi.img"
    gp.singleoutputmapalgebra_sa(express, out)

except:

print gp.GetMessages()

print "
done"

-Marcus

DanPatterson_Retired · Answer

a pretty ugly expression...but I suspect that you need to "float" one of the terms since integer division is occurring and everything is getting truncated to 0 or 1.  ie in Python
>>> print 1/3
0
>>> print 1.0/3
0.333333333333
>>>
surprised by the first result?