Intermediate raster leads to incorrect final output in arcpy map algebra

705
2
Jump to solution
08-16-2017 07:03 AM
JacobNederend
New Contributor II

I want to calculate an index and can't figure out why 2 seemingly identical scripts produce very different results. The input raster is a multiband GeoTiff (32-bit float) and I need to use the RGB bands as inputs. The RGB bands first need to be normalized to a scale of 0-1 (e.g. Red/Red_max), then the chromatic coordinates need to be determined (e.g. red/red+green+blue), and finally the index (ExG = 2*green-red-blue).

Since my data are already on the 0-1 scale, I can skip that step but still encounter problems. Can someone please explain why this:

red = Float(Raster(inraster+"/Band_1")/(Raster(inraster+"/Band_1")+Raster(inraster+"/Band_2")+Raster(inraster+"/Band_6")))
green = Float(Raster(inraster+"/Band_2")/(Raster(inraster+"/Band_1")+Raster(inraster+"/Band_2")+Raster(inraster+"/Band_6")))
blue = Float(Raster(inraster+"/Band_6")/(Raster(inraster+"/Band_1")+Raster(inraster+"/Band_2")+Raster(inraster+"/Band_6")))
ExG = 2*green - red - blue‍‍‍‍

produces the correct output. But the following does not:

red = Raster(inraster+"/Band_1")
green = Raster(inraster+"/Band_2")
blue = Raster(inraster+"/Band_6")
red_cc = Float(red/red+green+blue)
green_cc = Float(green/red+green+blue)
blue_cc = Float(blue/red+green+blue)
ExG = 2*green_cc - red_cc - blue_cc‍‍‍‍‍‍‍‍‍‍‍‍‍‍

I did the math for a single pixel by hand and the former is definitely correct. Eg.

The original inputs prior to finding the chromatic coordinates are: Red = 0.427023, Green = 0.405449, Blue =0.349515

The ExG calculates to 0.02907. 

When using the first method, I get 0.029070, but using the second I get 0.080466. 

Tags (2)
0 Kudos
1 Solution

Accepted Solutions
curtvprice
MVP Esteemed Contributor

The "/" operator in Python and arcpy map algebra does integer division, that is 5 / 3 = 1  (not 1.66) if both operands are integer rasters. You need to float one of the operands:

rgb = Float(red + green + blue)  # Float tool (not float)
red_cc = red / rgb
green_cc = green / rgb
blue_cc = blue / rgb‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Divide—Help | ArcGIS Desktop 

View solution in original post

0 Kudos
2 Replies
curtvprice
MVP Esteemed Contributor

The "/" operator in Python and arcpy map algebra does integer division, that is 5 / 3 = 1  (not 1.66) if both operands are integer rasters. You need to float one of the operands:

rgb = Float(red + green + blue)  # Float tool (not float)
red_cc = red / rgb
green_cc = green / rgb
blue_cc = blue / rgb‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Divide—Help | ArcGIS Desktop 

0 Kudos
JacobNederend
New Contributor II

Thanks! This worked so I am marking it as the answer. It did, however, introduce a new problem. The ExG map gets promoted to 64-bit double precision, which unfortunately doesn't work with the Binary Threshold function I need to apply from the Image Analysis Window. I posted about that previously, but never did manage to implement it myself:

https://community.esri.com/thread/197481-why-does-arcpy-produce-different-results-than-raster-calcul...  

0 Kudos