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-05-2017

import numpy as np
a = np.arange(4*6).reshape(4, 6)
b = np.zeros(4*6).reshape(4, 6)

for i in range (5):
    for j in range (7):
        v1 = a[i:j]
        v2 = a[i:j] + 2
        b[i:j] = v1 + v2

print("array a \n{}\narray b \n{}".format(a, b))‍‍‍‍‍‍‍‍‍‍‍

you seem to be mixing up commas and colons...

See this generic example... ie. it is a lesson, not a direct solution

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.000  4.000  6.000  8.000  10.000  12.000]
 [ 14.000  16.000  18.000  20.000  22.000  24.000]
 [ 26.000  28.000  30.000  32.000  34.000  36.000]
 [ 38.000  40.000  42.000  44.000  46.000  48.000]]‍‍‍‍‍‍‍‍‍‍

If you want array b to be integers, b = np.zeros(3*4, dtype='int').reshape(3,4)

XinXia1 · ‎01-06-2017

Thank you Dan! That was very instructive.

I now have a very small sample of the data:

array a 
[[ 2.50291276  2.50291276  2.50291276 ...,  2.20816231  2.20816231
   2.20816231]
 [ 2.50291276  2.50291276  2.50291276 ...,  2.20816231  2.20816231
   2.20816231]] 
 array b 
[[ 453.25619507  453.25619507  453.25619507 ...,  411.51553345
   411.51553345  411.51553345]
 [ 453.25619507  453.25619507  453.25619507 ...,  411.51553345
   411.51553345  411.51553345]]
array c 

[[ 455.75909424  455.75909424  455.75909424 ...,  413.72369385
   413.72369385  413.72369385]
 [ 455.75909424  455.75909424  455.75909424 ...,  413.72369385
   413.72369385  413.72369385]]
array c 

[[ 458.26202393  458.26202393  458.26202393 ...,  415.93185425
   415.93185425  415.93185425]
 [ 458.26202393  458.26202393  458.26202393 ...,  415.93185425
   415.93185425  415.93185425]]
array c 

[[ 460.7649231   460.7649231   460.7649231  ...,  418.14001465
   418.14001465  418.14001465]
 [ 460.7649231   460.7649231   460.7649231  ...,  418.14001465
   418.14001465  418.14001465]]‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

I'm not sure how to make the entire array print, but here's the 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)

a = numpy.zeros(4*6).reshape(4, 6)
b = numpy.zeros(4*6).reshape(4, 6)
c = 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 = v1
        b = v2

print("array a \n{} \n array b \n{}".format(a, b)) 

for x in range (1, 4):
    c = a * x
    c = c + b

    print "array c \n"
    print c
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

It definitely lets me loop through more than 2 iterations. Should I try saving it back to a raster to see what is causing me to get that 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.‍‍‍‍‍‍

DanPatterson_Retired · ‎01-06-2017

the line numbers dont match up now... where is 'lowerleft' defined?

just skip it, repost your full code once you have cleaned it up and use

arcpy.NumPyArrayToRaster(myArray, x_cell_size=1)

if that works, then the problem is with the lower left

ps get ready for python... print is not a function

ie

print("this is my array\n{}".format(myarray))‍

XinXia1 · ‎01-08-2017

Hi Dan,

I don't quite understand when you say print is not a function?

Do you mean never do this:

print myArray‍

Instead always do this:

print("this is my array\n{}".format(myarray))‍

If so what is wrong with doing this:

print myArray‍

This is my cleaned up code:

import arcpy
import numpy

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

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

a = numpy.zeros(4*6).reshape(4, 6)
b = numpy.zeros(4*6).reshape(4, 6)
c = 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 = v1
        b = v2 

print("array a \n{} \n array b \n{}".format(a, b)) 

for x in range (1, 4):
    c = a * x
    c = c + b

    print ("array c \n{}".format(c))

    newRas = arcpy.NumPyArrayToRaster(c, x_cell_size = 1)
    newRas.save("C:/PyTest/testResult.gdb/Year%d" %x)
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Why is does the resulting raster look like a thin line of pixels instead of a 4*6 shaped raster?

Does this bit of code read a 4*6 area of the original raster starting from the top left or just the first 24 cells?

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

        a = v1
        b = v2 ‍‍‍‍‍‍‍

Since I used this bit of code, does this mean my rasters have a different cell size than the original rasters?

newRas = arcpy.NumPyArrayToRaster(c, x_cell_size = 1)‍

DanPatterson_Retired · ‎01-08-2017

re print... the new syntax print( ) is for python 3.x and is supported by python 2.x... your code will break in the future if you don't plan ahead

a=v1

b=v2

if that is indeed what you are trying to do

re the cell size... get it to work first, then change the cell size to match

I am just trying to rule out potential sources of error in your original code

XinXia1 · ‎01-10-2017

Hi Dan,

I have changed my code to include:

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

This is the error I get:

Traceback (most recent call last):
  File "C:\PyTest\petCode.py", line 19, in <module>
    a = v1
ValueError: could not broadcast input array from shape (0,6318) into shape (0,6)
‍‍‍‍‍‍‍‍‍

Here is the code:

import arcpy
import numpy

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

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

a = numpy.zeros(4*6).reshape(4, 6)
b = numpy.zeros(4*6).reshape(4, 6)
c = 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{} \n array b \n{}".format(a, b)) 

for x in range (1, 4):
    c = a * x
    c = c + b

    print ("array c \n{}".format(c))

    newRas = arcpy.NumPyArrayToRaster(c, x_cell_size = 1)
    newRas.save("C:/PyTest/testResult.gdb/Year%d" %x)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

I thought I would test using a 4*6 clip of the the original (5400, 6318) raster starting from the top left corner in order to see if I can resolve that previous error.

DanPatterson_Retired · ‎01-10-2017

well... there have been 26 views to this thread... no one else has had any comments... and I still can't get a straight answer as to what you are doing or why you keep insisting on processing on a cell by cell basis... In this situation, it is time to bail and let the other viewers leap in with their interpretation of what is being done..

So in summary... you dont work with numpy arrays on a cell-by-cell basis... you use vectorized operations at worse or full array and/or scalar operations as normal. Once you or someone can clarify this situation I will defer to the viewing gallery to leap in with their insights

XinXia1 · ‎01-10-2017

Thank you for all the time you have put in to help me, Dan. Thank you!

When I did this:
arrayA = arrayA * x

arrayC = arrayB + arrayA

where x multiplies the entire arrayA instead of going through cell by cell I could not get x to loop beyond 2.

So I wanted to use a smaller section of the original array to post as a sample to test. But I didn't understand that trying to get a smaller section was still going through an array cell by cell.

I don't want to go through an array cell by cell. But nothing else was working. Why does the code only work when going cell by cell.

DanPatterson_Retired · ‎01-10-2017

this is from a related thread... here is how to work with rasters and conditional statement

https://community.esri.com/thread/188417-how-to-make-raster-conditionals-more-efficient