AnsweredAssumed Answered

Why does Arcpy produce different results than raster calculator/map algebra?

Question asked by JakeNederend on Jun 29, 2017
Latest reply on Aug 20, 2017 by lpinner

I need to do binary thresholding on several rasters of the Color Index of Vegetation Extraction (CIVE). First, I calculate the index in Arcpy, then I load each up raster individually in ArcMap to run binary thresholding from the Image Analysis Window. Unfortunately, I get an entirely black image. 


So, I then selected one date and calculated the CIVE manually in the Raster Calculator (this takes a long time since each of the bands in the CIVE formula need to be modified before calculating the index). The binary thresholding image from this raster is perfect. Vegetation (maize) pixels are clearly identified. 


As far as I can tell, the CIVE from ArcPy and the CIVE from ArcMap are identical (see attached images). I cannot figure out what is going on here, and really can't afford the time to do all this manually because I have the same problem with the Excess Red Index (ExG) and the Excess Green Minus Excess Red Index (ExGR). Oddly, I DON'T have this problem with the Excess Green Index even though it is run in the same script...


Here is a code snippet for the CIVE calculation. Note that red=Band_1, Green=Band_2, and Blue=Band_6 because I merged the orthomosaics from two filter-modified digital cameras.


print "Calculating and saving CIVE..."
for raster in rasters:
    red = Float(Raster(raster+"/Band_1"))/(Raster(raster+"/Band_1")+Raster(raster+"/Band_2")+Raster(raster+"/Band_6"))
    green = Float(Raster(raster+"/Band_2"))/(Raster(raster+"/Band_1")+Raster(raster+"/Band_2")+Raster(raster+"/Band_6"))
    blue = Float(Raster(raster+"/Band_6"))/(Raster(raster+"/Band_1")+Raster(raster+"/Band_2")+Raster(raster+"/Band_6"))
    outcive = Float((0.441*red)-(0.881*green)+(0.385*blue)+18.78745)"//"+"CIVE_"+raster[-6:-4]+".tif")
    print outcive
print "Done..."

Binary thresholding of CIVE from ArcMap

CIVE from ArcMap

Binary thresholding of CIVE from ArcPy

CIVE from ArcPy

Example pixel showing that ArcMap and ArcPy CIVEs are identical

An example pixel showing that the values are identical for both CIVE calculation methods