topic Re: Looping through rasters only give correct result on first iteration in Python Questions

Looping through rasters only give correct result on first iteration

XinXia1 — Sat, 11 Dec 2021 13:50:52 GMT

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)

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

Re: Error with Raster Calculations and the range() Function in Loop

DanPatterson_Retired — Sat, 11 Dec 2021 13:50:55 GMT

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.

Re: Error with Raster Calculations and the range() Function in Loop

XinXia1 — Mon, 02 Jan 2017 06:01:20 GMT

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?

Re: Error with Raster Calculations and the range() Function in Loop

DanPatterson_Retired — Mon, 02 Jan 2017 06:53:58 GMT

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

Re: Error with Raster Calculations and the range() Function in Loop

XinXia1 — Mon, 02 Jan 2017 18:35:54 GMT

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?

Re: Error with Raster Calculations and the range() Function in Loop

DanPatterson_Retired — Mon, 02 Jan 2017 22:09:45 GMT

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?

Re: Error with Raster Calculations and the range() Function in Loop

XinXia1 — Sat, 11 Dec 2021 13:50:58 GMT

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.

Re: Looping through rasters only give correct result on first iteration

DanPatterson_Retired — Sat, 11 Dec 2021 13:51:01 GMT

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?

Re: Looping through rasters only give correct result on first iteration

XinXia1 — Sat, 11 Dec 2021 13:51:04 GMT

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

Re: Looping through rasters only give correct result on first iteration

DanPatterson_Retired — Sat, 11 Dec 2021 13:51:07 GMT

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]])‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Re: Looping through rasters only give correct result on first iteration

DanPatterson_Retired — Sat, 11 Dec 2021 13:51:10 GMT

here is a simpler version explaining what happens to x, i and j

a = np.arange(4*6).reshape(4, 6)[:3,:3]
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 x in range(3):
    for i in range(2):
        for j in range(2):
            v1 = a[i,j]
            v1 = v1*x
            v2 = b[i,j]
            tot = v1 + v2
            print("x... {} v1*x +  v2 =>  {} + {} = {}".format(x, v1, v2, tot))
            result.append([x, i, j, v1, v2, tot])
result = np.array(result)
print("final result...\n x, i, j, v1, v2, tot\n{}".format(result))‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

with the result

final result...
 x, i, j, v1, v2, tot
[[ 0  0  0  0  2  2]
 [ 0  0  1  0  3  3]
 [ 0  1  0  0  8  8]
 [ 0  1  1  0  9  9]
 [ 1  0  0  0  2  2]
 [ 1  0  1  1  3  4]
 [ 1  1  0  6  8 14]
 [ 1  1  1  7  9 16]
 [ 2  0  0  0  2  2]
 [ 2  0  1  2  3  5]
 [ 2  1  0 12  8 20]
 [ 2  1  1 14  9 23]]‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Re: Looping through rasters only give correct result on first iteration

XinXia1 — Sat, 11 Dec 2021 13:51:12 GMT

Hi Dan,

I know going through an array cell by cell is very inefficient. So the first code I tried did not got through my array cell by cell. The first code I tried looked 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)


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

This code would also give me the incorrect values after the first iteration, and it would crash after only 2 iterations giving me an error 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.‍‍‍‍‍‍

I don't know how to fix this error. Are my arrays to big to fit into RAM? Do I have to clear or release the memory after each iteration? Do I need to go multiprocessing?

I will try to subsample my arrays and return with the results.

Re: Looping through rasters only give correct result on first iteration

DanPatterson_Retired — Sat, 11 Dec 2021 13:51:15 GMT

You don't go through arrays cell by cell.. If you want to 'do something' to an array, you apply the function or whatever all at once. ... last crack..

>>> a = np.arange(3*4).reshape(3,4)
>>> a
array([[ 0, 1, 2, 3],
       [ 4, 5, 6, 7],
       [ 8, 9, 10, 11]])
>>> a*3
array([[ 0, 3, 6, 9],
       [12, 15, 18, 21],
       [24, 27, 30, 33]])
>>> (a * 3)/a.max()
array([[ 0.000, 0.273, 0.545, 0.818],
       [ 1.091, 1.364, 1.636, 1.909],
       [ 2.182, 2.455, 2.727, 3.000]])‍‍‍‍‍‍‍‍‍‍‍‍‍
‍‍‍‍‍‍‍‍‍‍‍‍‍

Re: Looping through rasters only give correct result on first iteration

XinXia1 — Sat, 11 Dec 2021 13:51:18 GMT

Hi Dan,

Thank you for all your help. I see this code as NOT go through the array cell by cell. It is doing exactly what you are showing above. I AM multiplying the entire array inArrayPETcgcm by the constant x. If inArrayPETcgcm was array A like you have above. Then array A would have:

array dimensions (ndim) = 2

array shape (shape) = (5400, 6318)

Loop x just goes from 1 to 50 as the constants that get multiplied to the entire array.

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

How are you reading this code as going through the array cell by cell? I'm honestly so confused? If I am wrong please let me know.

Re: Looping through rasters only give correct result on first iteration

DanPatterson_Retired — Tue, 03 Jan 2017 23:50:34 GMT

the whole slicing with for i and for j that is in all your code until this last one. that is the point I was trying to make

Re: Looping through rasters only give correct result on first iteration

XinXia1 — Sat, 11 Dec 2021 13:51:20 GMT

But the last version of my code crashes, gives me this error, and is still only correct on the first iteration:

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

Re: Looping through rasters only give correct result on first iteration

DanPatterson_Retired — Sat, 11 Dec 2021 13:51:23 GMT

unless it is a structural problem with the array or the saving location... try a shorter file name

x = 1
>>> fn = "c:/somefolder/somesub/somegdb.gdb/X{}db".format(x)
>>> fn
'c:/somefolder/somesub/somegdb.gdb/X1db'‍‍‍‍

particularly if it is an esri grid... since filenames > about 10 (13?) are not permitted

Re: Looping through rasters only give correct result on first iteration

XinXia1 — Sat, 11 Dec 2021 13:51:26 GMT

Hi Dan,

Thank you so much for all your help. I changed to shorter file names, but I am getting the same error.

Error:

Traceback (most recent call last):
  File "C:\PyTest\petCode.py", line 19, in <module>
    newRasPET = arcpy.NumPyArrayToRaster(outArrayPE, 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.‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

The code now makes it 2 iterations with correct raster values thank you for showing me the zeros_like() to properly initialize an array.

Latest code:

import arcpy
import numpy

inRasCgcm = arcpy.Raster("C:/PyTest/testResult.gdb/petCRsClip")
inRasNorm = arcpy.Raster("C:/PyTest/testResult.gdb/petNRsClip")

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

inArrayCgcm = arcpy.RasterToNumPyArray(inRasCgcm)
inArrayNorm = arcpy.RasterToNumPyArray(inRasNorm)

outArrayPE = numpy.zeros_like(inArrayNorm)

for x in range (1, 4):
    inArrayCgcm = inArrayCgcm * x
    outArrayPE = inArrayNorm + inArrayCgcm

    newRasPET = arcpy.NumPyArrayToRaster(outArrayPE, lowerLeft, cellSize)
    newRasPET.save("C:/PyTest/testResult.gdb/petR{}Year".format(x))‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

The structure of my array and saving location should be okay if I can make it to 2 iterations correctly shouldn't it?

Re: Looping through rasters only give correct result on first iteration

DanPatterson_Retired — Wed, 04 Jan 2017 22:02:20 GMT

if it makes through 2, then something is wrong with the data... as suggested, slice it to a smaller area, post a sample, or just use the subset for testing

Re: Looping through rasters only give correct result on first iteration

XinXia1 — Sat, 11 Dec 2021 13:51:29 GMT

Hi Dan

I'm not good with python yet. I was trying this code to post a sample, but it gave me this error. I think it's something simple that I'm not understanding with python.

Error:

Traceback (most recent call last):
  File "C:\PyTest\petCode.py", line 23, in <module>
    a = v1
IndexError: index 6 is out of bounds for axis 0 with size 6‍‍‍‍

Code:

import arcpy
import numpy

inRasCgcm = arcpy.Raster("C:/PyTest/testResult.gdb/petCRsClip")
inRasNorm = arcpy.Raster("C:/PyTest/testResult.gdb/petNRsClip")

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

inArrayCgcm = arcpy.RasterToNumPyArray(inRasCgcm)
inArrayNorm = arcpy.RasterToNumPyArray(inRasNorm)

#outArrayPE = numpy.zeros_like(inArrayNorm)

a = numpy.zeros(4*6).reshape(4, 6)
b = numpy.zeros(4*6).reshape(4, 6)

for i in range (5):
    for j in range (7):
        v1 = inArrayCgcm[i, j]
        v2 = inArrayNorm[i, j]

        a[i][j] = v1
        b[i][j] = v2

print "array a \n"
print a

print "array b \n"
print b‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍