AnsweredAssumed Answered

Looping through rasters only give correct result on first iteration

Question asked by clair0girl on Jan 1, 2017
Latest reply on Jan 10, 2017 by Dan_Patterson


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)"C:/VegModel/Revise.gdb/pet100ReviseYear%dB" %x)