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

[JacobNederend]

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)
    outcive.save(CIVEFolder+"//"+"CIVE_"+raster[-6:-4]+".tif")
    print outcive
print "Done..."

Binary thresholding of CIVE from ArcMap

Binary thresholding of CIVE from ArcPy

Example pixel showing that ArcMap and ArcPy CIVEs are identical

[JacobNederend]

Dan,

Thanks so much for this. I have not yet had time to run it (my Python skills have not yet extended to defining functions) but I trust that it is correct so I am marking it as the answer. 

[curtvprice]

Hi Jacob,

I have a guess that Raster Calculator may be different then arcpy because you may be running Raster Calculate from ArcMap  (in x32) and running arcpy in x64 Python you may get different results.

If the bit depth of your arcpy map algebra output is not what you need, you may want to try using the Copy Raster function instead of running .save(), which allows you to specify the bit depth and NoData value. Or maybe try running the script using arcpy x32. (C:\Python27\ArcGIS10.5\python.exe).

[Luke_Pinner]

You probably just need to add a "Calculate Statistics" step in your script.