topic Re: Why does Arcpy produce different results than raster calculator/map algebra? in Python Questions

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

JacobNederend — Sat, 11 Dec 2021 07:39:01 GMT

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

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

DanPatterson_Retired — Thu, 29 Jun 2017 13:06:31 GMT

could be a symbology issue, what is the nodata value for both and how do their symbologies compare? Also check the raster type (since the bit type may have been promoted in one)

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

JacobNederend — Sat, 01 Jul 2017 02:11:40 GMT

Thanks for the advice, Dan. I think you are right about the bit-depth because the NoData values and symbologies are identical. The Excess Green Index I produced is 32-bit floating point, whereas the CIVE and all other indices are 64-bit double-precision.

Could it be possible that the decimal coefficients in the CIVE index are messing up the bit-depth? If so, how could I fix the outputs to be 32-bit floating point? Some formulas are below for comparison.

CIVE = 0.441r - 0.881g +0.385b +18.78745

ExG = 2g-r-b

ExR = 1.4r-b

ExGR = Exg-ExR

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

JacobNederend — Sat, 01 Jul 2017 02:33:57 GMT

Also, the NoData value is -1.797693e+308.

I also found this discussion from 2015 which you were also involved in. This is almost certainly my problem, but the OP never updated with a solution.

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

DanPatterson_Retired — Sat, 01 Jul 2017 02:57:43 GMT

I suspect you are going to have to deal with the Nodata (Environment Settings) in your script and 'Map Down' I think, but it has been some time since I normally work with numpy arrays where you have direct explicit control over the data type.

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

JacobNederend — Sat, 11 Dec 2021 07:39:03 GMT

This didn't work, unfortunately. Numpy would be better but am not sure how I would go about it. Here's what I am thinking:

# assume that I already set up a list called composites, which are 6-band GeoTiffs for which I need to use bands 1,2, and 6
# Also assume that I have lists for all the relevant geotransform properties in order to convert output back to raster 
#N ot sure what to set NoData to in RasterToNumpyArray
nparrays = []
for i, raster in enumerate(composites):
    nparrays.append(arcpy.RasterToNumPyArray(raster, nodata_to_value=raster.noDataValue).astype(float))
# Apply the CIVE formula of 0.441r - 0.881g +0.385b +18.78745 
CIVE_arrays = []
for date, numpyarrays in enumerate(nparrays):
    CIVE = np.zeros_like (numpyarrays, np.float32)
    # **not sure how to actually run the formula from here** because I don't know how to select dimensions 0,1, and 5
    for i in range (0,5) #
        CIVE[i,:,:] = 0.441*[i=0]-0.881[i=2]+0.385[i=6]+18.78745   # obviously the i= are invalid syntax

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

DanPatterson_Retired — Sat, 11 Dec 2021 07:39:06 GMT

An example... I used an array with a shape of (5, 5, 6)... 5 rows, 5 columns and 6 bands... This can just as easily be done using reshaped array with the bands first, then rows and cols (ie (6, 5, 5)

a = np.zeros(shape=(5,5,6), dtype='float32')
a[..., 0] = 1
a[..., 1] = 2
a[..., 5] = 6
r = a[..., 0]      # = 1 from 'a' definition
g = a[..., 1]
b = a[..., 5]
args = [r, g, b]
print("Band 1\n{}\nBand 2\n{}\nBand 6\n{}".format(*args))
result = 0.441*r - 0.881*g +0.385*b + 18.78745
print("Result...\n{}".format(result))‍‍‍‍‍‍‍‍‍‍‍

which would look like this as results

Band 1
[[ 1.0  1.0  1.0  1.0  1.0]
 [ 1.0  1.0  1.0  1.0  1.0]
 [ 1.0  1.0  1.0  1.0  1.0]
 [ 1.0  1.0  1.0  1.0  1.0]
 [ 1.0  1.0  1.0  1.0  1.0]]
Band 2
[[ 2.0  2.0  2.0  2.0  2.0]
 [ 2.0  2.0  2.0  2.0  2.0]
 [ 2.0  2.0  2.0  2.0  2.0]
 [ 2.0  2.0  2.0  2.0  2.0]
 [ 2.0  2.0  2.0  2.0  2.0]]
Band 6
[[ 6.0  6.0  6.0  6.0  6.0]
 [ 6.0  6.0  6.0  6.0  6.0]
 [ 6.0  6.0  6.0  6.0  6.0]
 [ 6.0  6.0  6.0  6.0  6.0]
 [ 6.0  6.0  6.0  6.0  6.0]]
Result...
[[ 19.8  19.8  19.8  19.8  19.8]
 [ 19.8  19.8  19.8  19.8  19.8]
 [ 19.8  19.8  19.8  19.8  19.8]
 [ 19.8  19.8  19.8  19.8  19.8]
 [ 19.8  19.8  19.8  19.8  19.8]]
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

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

DanPatterson_Retired — Sat, 11 Dec 2021 07:39:09 GMT

deja vu.... copyraster or numpy... who would have thought ... memory is slipping ... enjoy the example

EDIT

I should add if you want to change the dtype back to something before saving, then you can do the conversion by making an array of the desired dtype and shape and copying the values over into it. Here is an example using the cive changing out to an int32 array rather than the float64 array you would get

cive_2 = np.zeros(cive.shape, dtype='int32')
cive_2[:] = cive

cive result...
[[ 19.4  19.4  19.4  19.4  19.4]
 [ 19.4  19.4  19.4  19.4  19.4]
 [ 19.4  19.4  19.4  19.4  19.4]
 [ 19.4  19.4  19.4  19.4  19.4]
 [ 19.4  19.4  19.4  19.4  19.4]]

cive_2 result...
[[19 19 19 19 19]
 [19 19 19 19 19]
 [19 19 19 19 19]
 [19 19 19 19 19]
 [19 19 19 19 19]]

cive.dtype, cive_2.dtype
(dtype('float64'), dtype('int32'))‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

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

JacobNederend — Sat, 11 Dec 2021 07:39:12 GMT

I made a mistake in the formula and forgot to carry over the normalization (e.g. r/r+g+b) from my original spatial analyst script. Altering your example with the division yields a divide by zero error:

numpyarrays = arcpy.RasterToNumPyArray(infile, nodata_to_value=infile.noDataValue).astype(float)

# Apply the CIVE formula of 0.441r - 0.881g +0.385b +18.78745
r = numpyarrays[0,:,:]
g = numpyarrays[1,:,:]
b = numpyarrays[5,:,:]
args = [r, g, b]
print("Band 1\n{}\nBand 2\n{}\nBand 6\n{}".format(*args))
CIVE = 0.441*(r/r+g+b) - 0.881*(g/r+g+b) +0.385*(b/r+g+b) + 18.78745
print("Result...\n{}".format(CIVE))
#----------------------------------------------------------------------------------------------------------------------------------
C:/Users/Default.Default-PC/.PyCharmEdu3.5/config/scratches/scratch_5.py:48: RuntimeWarning: invalid value encountered in divide
  CIVE = 0.441*(r/r+g+b) - 0.881*(g/r+g+b) +0.385*(b/r+g+b) + 18.78745
C:/Users/Default.Default-PC/.PyCharmEdu3.5/config/scratches/scratch_5.py:48: RuntimeWarning: divide by zero encountered in divide
  CIVE = 0.441*(r/r+g+b) - 0.881*(g/r+g+b) +0.385*(b/r+g+b) + 18.78745
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

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

DanPatterson_Retired — Sat, 11 Dec 2021 07:39:15 GMT

Numpy done, png in... tif out, The middle line is green in the lower figure believe it or not

def num_109():
    """raster band math... just a quick demo
    : note: Bands are number from 0 to n, 
    : references
    : http://conteudo.icmc.usp.br/pessoas/moacir/papers/Ponti_GRSL2013.pdf
    : http://www.sciencedirect.com/science/article/pii/S2214317315000347
    :
    :a = np.zeros(shape=(5, 5, 6), dtype='int32')
    a[..., 0] = vals[0]  # 1
    a[..., 1] = vals[1]  # 2
    a[..., 5] = vals[2]  # 6
    """
    import arcpy
    f = r'C:\Book_Materials\images\rgb.png'  # an image
    arcpy.overwriteOutput = True
    bands = arcpy.RasterToNumPyArray(f)
    a = bands[:3, :]
    r = bands[0, :]
    g = bands[1, :]
    b = bands[2, :]
    k = bands[3, :]
    # normalize
    denom = (r + g + b)
    denom = denom.astype('float')
    r1 = r/(denom)
    g1 = g/(denom)
    b1 = b/(denom)
    args = [r, g, b]
    print("Band 1\n{}\nBand 2\n{}\nBand 6\n{}".format(*args))
    cive = (0.441 * r) - (0.881 * g) + (0.385 * b) + 18.78745
    cive_2 = np.zeros(cive.shape, dtype='int32')
    cive_2 = np.zeros(cive.shape, dtype='int32')
    #exg = 2*g1 - r1 - b1
    print("CIVE result...\n{}\n ExG result...\n{}".format(cive, cive_2))
    out = arcpy.NumPyArrayToRaster(cive_2)
    out.save(f.replace('.png', '1.tif'))
    return a, r, g, b, k

# __main__ .... code section
if __name__ == "__main__":
    a, r, g, b, k = num_109()  # raster band math demo‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

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

JacobNederend — Tue, 04 Jul 2017 12:46:15 GMT

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.

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

curtvprice — Sat, 19 Aug 2017 03:13:29 GMT

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).

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

Luke_Pinner — Sun, 20 Aug 2017 23:23:17 GMT

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