Looping through rasters only give correct result on first iteration

[XinXia1]

Hi,

I'm new to python. I need to generate a new raster C from two rasters A and B by first multiplying a constant value with all elements of raster A and then adding A and B. I need to generate raster C for multiple years so I have it looping over the years I need. However, line 41 below is not giving the correct value after multiplication except on the first iteration through the loop at line 27. Can anyone explain why that might be? Why does this code only work on the first iteration? 

EDIT: the code also works without using the Decimal module but still only gives the correct raster on the first iteration

import arcpy
import numpy
import decimal

#inRas = arcpy.Raster("C:/VegModel/CCSSRawData.gdb/vegnorm")
#lowerLeft = arcpy.Point(inRas.extent.XMin,inRas.extent.YMin)
#cellSize = inRas.meanCellWidth

inRasPETcgcm = arcpy.Raster("C:/VegModel/CCSSRawData.gdb/petCgcmRS100Clip")
inRasPETnorm = arcpy.Raster("C:/VegModel/CCSSRawData.gdb/petNormRS100Clip")

lowerLeft = arcpy.Point(inRasPETnorm.extent.XMin,inRasPETnorm.extent.YMin)
cellSize = inRasPETnorm.meanCellWidth

inArrayPETcgcm = arcpy.RasterToNumPyArray(inRasPETcgcm)
#inArrayPETcgcm = inArrayPETcgcm * 30
inArrayPETnorm = arcpy.RasterToNumPyArray(inRasPETnorm)

rasHeight, rasWidth = inArrayPETnorm.shape

#outArrayPET = inArrayPETnorm + inArrayPETcgcm

outArrayPET = inArrayPETnorm

decimal.getcontext().prec = 20

for x in range (0, 3):
    for i in range(0, 2):
        for j in range(0, 2):
            value1 = str(inArrayPETcgcm[i, j])
            value2 = str(inArrayPETnorm[i, j])
            loopValue = str(x)
            
            
            #print "value1 is: "
            #print value1
        
            #print "value2 is: "
            #print value2

            value1 = decimal.Decimal(value1) * decimal.Decimal(loopValue)

            #print "value1 after multiplication is: "
            #print value1
            value3 = value1 + decimal.Decimal(value2)

            #print "value3 is: "
            #print value3

            outArrayPET[i, j] = value3
    
    newRasPET = arcpy.NumPyArrayToRaster(outArrayPET, lowerLeft, cellSize)
    newRasPET.save("C:/VegModel/Revise.gdb/pet100ReviseYear%dB" %x)

‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[DanPatterson_Retired]

Why are you using the Decimal module?  You already have a numpy array, just continue doing your work with it then save it back out to arcmap if you really need to.  

ft = {'bool': lambda x: repr(x.astype('int32')),
      'float': '{: 0.10f}'.format}
np.set_printoptions(edgeitems=3, linewidth=100, precision=10,
                    suppress=True, threshold=100,
                    formatter=ft)
a = np.arange(9.).reshape(3,3)
b = a* 0.2
a
array([[ 0.0000000000,  1.0000000000,  2.0000000000],
       [ 3.0000000000,  4.0000000000,  5.0000000000],
       [ 6.0000000000,  7.0000000000,  8.0000000000]])
b
array([[ 0.0000000000,  0.2000000000,  0.4000000000],
       [ 0.6000000000,  0.8000000000,  1.0000000000],
       [ 1.2000000000,  1.4000000000,  1.6000000000]])
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

The dtype is float64 which given your equations, I am not sure Decimal is needed given that its use can be circumvented in many situations by scaling and/or use of integer math.

[XinXia1]

Hi Dan,

Thank you for the tip. I did try it also without using the Decimal module, which works fine. However, my problem is that I am only getting a correct raster result in the first loop through at line 27. All subsequent iterations gives me the wrong values in the rasters, which I have checked by using the identify tool and picking random values to do the math. I should also mention that my actual raster is something like 6000 by 3000. I'm pretty sure arcpy can handle the math of multiplication (line 41) and addition (line 45), but I have no idea why only the first iteration give a correct result. Is it a memory problem? Or do I need to clear the values after each iteration?

[DanPatterson_Retired]

well if you are trying to save each iteration to a new file, then your last 2 lines need to be indented more since they are outside of the loop it appears.

And again, dump the Decimal module... it is giving you false precision.  Rasters may be 32 bit floats (ie. in the case of esri grids)..., numpy arrays 64 bit at least then back to a raster... 

Finally, if you are working with arrays, you don't work on a cell-by-cell basis as in my example...I haven't examined your code closely, to see what its intent was...just got caught up on the false precision

[XinXia1]

Hi Dan,

Thank you, I understand what you are staying about not needing the Decimal module. However, I have run it without the using Decimal and I am still only getting a correct raster on the first iteration of the loop at line 27. I want to save each iteration to a new file after the outer most for loop at line 27. If I indent it more doesn't it mean I'll be saving each iteration to a new file on the inner nested for loops? 

[DanPatterson_Retired]

And did you try it?  As I said, I am not clear on the intent, therefore not clear on the desired output... all I do know is do you want the file saved after the j loop is completed or the i loop?

[XinXia1]

Hi Dan,

I want the file saved after the i loop is completed or after each iteration of the x loop. The intent is that I have raster A and B. I need to multiple 1, 2, 3, ... 48, 49, 50 to raster A and add the result to raster B. I have outlined what I have tried below.

The initial code I tried is this:

import arcpy
import numpy

inRasPETcgcm = arcpy.Raster("C:/VegModel/CCSSRawData.gdb/petCgcmRS100Clip")
inRasPETnorm = arcpy.Raster("C:/VegModel/CCSSRawData.gdb/petNormRS100Clip")

lowerLeft = arcpy.Point(inRasPETnorm.extent.XMin,inRasPETnorm.extent.YMin)
cellSize = inRasPETnorm.meanCellWidth

inArrayPETcgcm = arcpy.RasterToNumPyArray(inRasPETcgcm)

inArrayPETnorm = arcpy.RasterToNumPyArray(inRasPETnorm)


for x in range (1, 51):
    inArrayPETcgcm = inArrayPETcgcm * x

    outArrayPET = inArrayPETnorm + inArrayPETcgcm
    
    newRasPET = arcpy.NumPyArrayToRaster(outArrayPET, lowerLeft, cellSize)
    newRasPET.save("C:/VegModel/Revise.gdb/pet100ReviseYear%dB" %x)

‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

But that code crashes after 2 iterations and give an error message like this:

Traceback (most recent call last):
  File "C:\VegModel\PETyearlyChange.py", line 20, in <module>
    newRasPET = arcpy.NumPyArrayToRaster(outArrayPET, lowerLeft, cellSize)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\__init__.py", line 2292, in NumPyArrayToRaster
    return _NumPyArrayToRaster(*args, **kwargs)
TypeError: Cannot create raster for numpy array.‍‍‍‍‍‍

Plus, of the 2 rasters that got saved before crashing only the raster in the first iteration is correct. So then I tried code that looks like this:

import arcpy
import numpy

inRasPETcgcm = arcpy.Raster("C:/VegModel/CCSSRawData.gdb/petCgcmRS100Clip")
inRasPETnorm = arcpy.Raster("C:/VegModel/CCSSRawData.gdb/petNormRS100Clip")

lowerLeft = arcpy.Point(inRasPETnorm.extent.XMin,inRasPETnorm.extent.YMin)
cellSize = inRasPETnorm.meanCellWidth

inArrayPETcgcm = arcpy.RasterToNumPyArray(inRasPETcgcm)

inArrayPETnorm = arcpy.RasterToNumPyArray(inRasPETnorm)

rasHeight, rasWidth = inArrayPETnorm.shape

outArrayPET = inArrayPETnorm

for x in range (1, 51):
    for i in range(0, 2):
        for j in range(0, 2):
            value1 = inArrayPETcgcm[i, j]
            value2 = inArrayPETnorm[i, j]

            value1 = value1 * x

            value3 = value1 + value2

            outArrayPET[i, j] = value3
    
    newRasPET = arcpy.NumPyArrayToRaster(outArrayPET, lowerLeft, cellSize)
    newRasPET.save("C:/VegModel/Revise.gdb/pet100ReviseYear%dB" %x)

‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

This code is able to make it to 50 without crashing. It also doesn't use Decimal. But again, only the raster from the first iteration is correct.

[DanPatterson_Retired]

a simpler version of what you have so far..

rows = 4
cols = 6
N = 5    # your 51
a = np.arange(rows*cols).reshape(rows, cols)
b = np.arange(rows*cols).reshape(rows, cols) + 2
c = np.zeros_like(b)
print("Input arrays... a and b\n{}\n...\n{}".format(a, b))
for x in range (1, N):
    for i in range(0, 2):
        for j in range(0, 2):
            value1 = a[i, j]
            value2 = b[i, j]
            value1 = value1 * x
            value3 = value1 + value2
            c[i][j] = value3     # note change
    print("c at {}, {}\n{}".format(i, j, c))‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

and the results so far

Input arrays... a and b
[[ 0  1  2  3  4  5]
 [ 6  7  8  9 10 11]
 [12 13 14 15 16 17]
 [18 19 20 21 22 23]]
...
[[ 2  3  4  5  6  7]
 [ 8  9 10 11 12 13]
 [14 15 16 17 18 19]
 [20 21 22 23 24 25]]
c at 1, 1
[[ 2  4  0  0  0  0]
 [14 16  0  0  0  0]
 [ 0  0  0  0  0  0]
 [ 0  0  0  0  0  0]]
c at 1, 1
[[ 2  5  0  0  0  0]
 [20 23  0  0  0  0]
 [ 0  0  0  0  0  0]
 [ 0  0  0  0  0  0]]
c at 1, 1
[[ 2  6  0  0  0  0]
 [26 30  0  0  0  0]
 [ 0  0  0  0  0  0]
 [ 0  0  0  0  0  0]]
c at 1, 1
[[ 2  7  0  0  0  0]
 [32 37  0  0  0  0]
 [ 0  0  0  0  0  0]
 [ 0  0  0  0  0  0]]
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

which isn't what you want.  So if the output array c, which is your outArrayPet which is a clone of inarrayPetNorm.

Could you provide some insight on...

• the array dimensions  (ie a.ndim)
• What is each array's shape (ie a.shape )
• are you trying to perform some sort of running calculation on the array?  

[XinXia1]

I'm not sure what you mean when you say running calculation on the array. I'm trying to do exactly what your code is showing.

Each array contains floats and have the following:

array dimensions (ndim) = 2

array shape (shape) = (5400, 6318)

I actually don't want my array C (outArrayPET) to be a clone of inarrayPETnorm I didn't know how to initialize it. I'll try doing what you have with:

c = np.zeros_like(b)

Also could you explain what the + 2 is for in line 5 of your code:

b = np.arange(rows*cols).reshape(rows, cols) + 2

I ran this code again with print statements: 

import arcpy
import numpy

inRasPETcgcm = arcpy.Raster("C:/VegModel/CCSSRawData.gdb/petCgcmRS100Clip")
inRasPETnorm = arcpy.Raster("C:/VegModel/CCSSRawData.gdb/petNormRS100Clip")

lowerLeft = arcpy.Point(inRasPETnorm.extent.XMin,inRasPETnorm.extent.YMin)
cellSize = inRasPETnorm.meanCellWidth

inArrayPETcgcm = arcpy.RasterToNumPyArray(inRasPETcgcm)

inArrayPETnorm = arcpy.RasterToNumPyArray(inRasPETnorm)

outArrayPET = inArrayPETnorm

for x in range (1, 3):
    for i in range(0, 2):
        for j in range(0, 2):
            value1 = inArrayPETcgcm[i, j]
            value2 = inArrayPETnorm[i, j]
            
            print "value1 before calculation is: "
            print value1
        
            print "value2 is before calculation: "
            print value2

            value1 = value1 * x

            print "value1 after multiplication is: "
            print value1
            value3 = value1 + value2

            print "value3 after addition is: "
            print value3

            outArrayPET[i, j] = value3
    
    newRasPET = arcpy.NumPyArrayToRaster(outArrayPET, lowerLeft, cellSize)
    newRasPET.save("C:/VegModel/Revise.gdb/pet100ReviseYear%dB" %x)


‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

value2 has somehow changed in the second iteration (the bottom 4 groups) of the x loop. The output gave me the following:

value1 before calculation is: 
2.50291
value2 is before calculation: 
453.256
value1 after multiplication is: 
2.50291275978
value3 after addition is: 
455.759107828


value1 before calculation is: 
2.50291
value2 is before calculation: 
453.256
value1 after multiplication is: 
2.50291275978
value3 after addition is: 
455.759107828


value1 before calculation is: 
2.50291
value2 is before calculation: 
453.256
value1 after multiplication is: 
2.50291275978
value3 after addition is: 
455.759107828


value1 before calculation is: 
2.50291
value2 is before calculation: 
453.256
value1 after multiplication is: 
2.50291275978
value3 after addition is: 
455.759107828


value1 before calculation is: 
2.50291
value2 is before calculation: 
455.759
value1 after multiplication is: 
5.00582551956
value3 after addition is: 
460.764919758


value1 before calculation is: 
2.50291
value2 is before calculation: 
455.759
value1 after multiplication is: 
5.00582551956
value3 after addition is: 
460.764919758


value1 before calculation is: 
2.50291
value2 is before calculation: 
455.759
value1 after multiplication is: 
5.00582551956
value3 after addition is: 
460.764919758


value1 before calculation is: 
2.50291
value2 is before calculation: 
455.759
value1 after multiplication is: 
5.00582551956
value3 after addition is: 
460.764919758

‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

I have no clue where the value2 of 455.759 is coming from in the second iteration of the x loop

[DanPatterson_Retired]

I am still not sure why you are cycling through one cell in an array.

Here is a simpler demo... 

import numpy as np

a = np.arange(4*6).reshape(4, 6)
result = []
b = a + 2     # just add 2 to array a to make b different than a
print("array a\n{}\narray b\n{}".format(a, b))
for i in range(3):
    for j in range(3):
        v1 = a[i,j]
        v2 = b[i,j]
        tot = v1 + v2
        print("v1 + v2 =>  {} + {} = {}".format(v1, v2, tot))
        result.append([i, j, v1, v2, tot])
result = np.array(result)
print("final result...\n i, j, v1, v2, tot\n{}".format(result))‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

and here is the result

array a
[[ 0  1  2  3  4  5]
 [ 6  7  8  9 10 11]
 [12 13 14 15 16 17]
 [18 19 20 21 22 23]]
array b
[[ 2  3  4  5  6  7]
 [ 8  9 10 11 12 13]
 [14 15 16 17 18 19]
 [20 21 22 23 24 25]]
v1 + v2 =>  0 + 2 = 2
v1 + v2 =>  1 + 3 = 4
v1 + v2 =>  2 + 4 = 6
v1 + v2 =>  6 + 8 = 14
v1 + v2 =>  7 + 9 = 16
v1 + v2 =>  8 + 10 = 18
v1 + v2 =>  12 + 14 = 26
v1 + v2 =>  13 + 15 = 28
v1 + v2 =>  14 + 16 = 30
final result...
 i, j, v1, v2, tot
[[ 0  0  0  2  2]
 [ 0  1  1  3  4]
 [ 0  2  2  4  6]
 [ 1  0  6  8 14]
 [ 1  1  7  9 16]
 [ 1  2  8 10 18]
 [ 2  0 12 14 26]
 [ 2  1 13 15 28]
 [ 2  2 14 16 30]]
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

I suspect this isn't what you want because of your slicing the arrays by

So again, since you only have a 2D array, why are you slicing out single cells.

EDIT

Why don't you slice out your array and show some real numbers with a result, say for the first 5 values for X.

You can subsample your arrays, using slicing as follows

>>> a
array([[ 0,  1,  2,  3,  4,  5],
       [ 6,  7,  8,  9, 10, 11],
       [12, 13, 14, 15, 16, 17],
       [18, 19, 20, 21, 22, 23]])
>>> a1 = a[:3,:3]
>>> b
array([[ 2,  3,  4,  5,  6,  7],
       [ 8,  9, 10, 11, 12, 13],
       [14, 15, 16, 17, 18, 19],
       [20, 21, 22, 23, 24, 25]])
>>> b1 = b[:3,:3]
>>> 
>>> 
>>> a1
array([[ 0,  1,  2],
       [ 6,  7,  8],
       [12, 13, 14]])
>>> b1
array([[ 2,  3,  4],
       [ 8,  9, 10],
       [14, 15, 16]])‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍