Looping through rasters only give correct result on first iteration

XinXia1 · ‎01-01-2017

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 · ‎01-03-2017

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

XinXia1 · ‎01-03-2017

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.

DanPatterson_Retired · ‎01-03-2017

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

XinXia1 · ‎01-03-2017

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.

DanPatterson_Retired · ‎01-03-2017

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

XinXia1 · ‎01-03-2017

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

DanPatterson_Retired · ‎01-03-2017

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

XinXia1 · ‎01-04-2017

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?

DanPatterson_Retired · ‎01-04-2017

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

XinXia1 · ‎01-05-2017

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