python map algebra for NDVI

[MarcusBeck]

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 "\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)

except:

    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.

Thanks,

Marcus

[DanPatterson_Retired]

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?

[MarcusBeck]

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 "\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 "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 "\ndone"

-Marcus

[WadeGivens]

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

Thanks,
Wade