Cumulative sum through a list of rasters

[CamiSunley]

Hi everyone, 

I have a list of rasters named: "ih1", "ih2", "ih3", ..., "ih28" that represent the interpolated rainfall height in the corresponding time step (hour). I want to calculate the cumulative rainfall height in each time step which is the sum of the rainfall in the same step and the cumulative rain of the previous one, e.g. CumRain(i) = ih(i) + CumRain(i-1). I am trying to code this in arcpy yet the "i" thing does not work for rasters.

If you have any idea how this could be solved, let me know.

Thanks a lot !

[DanPatterson]

Convert the raster s to numpy arrays

RasterToNumPyArray—ArcGIS Pro | Documentation

a0 
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11],
       [12, 13, 14, 15]])

a1
array([[16, 17, 18, 19],
       [20, 21, 22, 23],
       [24, 25, 26, 27],
       [28, 29, 30, 31]])

a2 
array([[32, 33, 34, 35],
       [36, 37, 38, 39],
       [40, 41, 42, 43],
       [44, 45, 46, 47]])

np.cumsum((a0, a1, a2), axis=0)

array([[[ 0,  1,  2,  3],
        [ 4,  5,  6,  7],
        [ 8,  9, 10, 11],
        [12, 13, 14, 15]],

       [[16, 18, 20, 22],
        [24, 26, 28, 30],
        [32, 34, 36, 38],
        [40, 42, 44, 46]],

       [[48, 51, 54, 57],
        [60, 63, 66, 69],
        [72, 75, 78, 81],
        [84, 87, 90, 93]]], dtype=int32)
​# ---- or  np.nancumsum(a, axis=0) if you have nodata cells

Now if you want an intermediate time step, just slice it

csum = np.nancumsum(a, axis=0)
csum[1]
array([[16, 18, 20, 22],
       [24, 26, 28, 30],
       [32, 34, 36, 38],
       [40, 42, 44, 46]], dtype=int32)

Then if you need raster s back

NumPyArrayToRaster—ArcGIS Pro | Documentation

[CamiSunley]

Thank you, I try the code below:

from arcpy import env

from arcpy import env

from arcpy.sa import *

import os

import numpy

arcpy.env.workspace = "C:/MSc_thesis/modf"

idw = arcpy.ListRasters("ih*", "GRID")

for i in idw:  

     inras = Raster(i)

     lowerLeft = arcpy.Point(inras.extent.XMin,inras.extent.YMin)

     cellSize = inras.meanCellWidth

     arr = arcpy.RasterToNumPyArray(inras,lowerLeft,cellSize, nodata_to_value=0)

     csum = numpy.cumsum(arr)

     newras = arcpy.NumPyArrayToRaster(csum,lowerLeft,cellSize,value_to_nodata=0)

but it returns the following error:

ValueError: Argument in_array: A two or three dimensional NumPy array is required.

[DanPatterson]

you aren't accumulating the rasters, hence you get one value 

arrays_in = []
for i in idw:  
     inras = Raster(i)
     lowerLeft = arcpy.Point(inras.extent.XMin,inras.extent.YMin)
     cellSize = inras.meanCellWidth
     arr = arcpy.RasterToNumPyArray(inras,lowerLeft,cellSize, nodata_to_value=0)
     arrays_in.append(arr)

main array = np.array(arrays_in)
c_sum = np.nancumsum(main_array, axis=0)
for i in c_sum:
    newras = arcpy.NumPyArrayToRaster(*** you need to provide an output name for "i"

.You have to get all the arrays together, then if you need individual output rasters, then you need to provide a name for each output raster.

Also, don't forget the "axis" for the cumulative sum.

[CamiSunley]

numpy.nancumsum returns error in my arcmap version so I replaced it with numpy.cumsum. Also, I might have an error in the code because now an ERROR 999998: Unexpected Error is appeared. I am not really familiar with all this stuff, sorry for the repetitive replies.

arrays_in = []
for i in idw:
    inras = Raster(i)
    lowerLeft = arcpy.Point(inras.extent.XMin,inras.extent.YMin)
    cellSize = inras.meanCellWidth
    arr = arcpy.RasterToNumPyArray(inras,lowerLeft,cellSize, nodata_to_value=0)
    arrays_in.append(arr)
main_arr = numpy.array(arrays_in)
c_sum = numpy.cumsum(main_arr, axis = 0)
for j in c_sum:
      newras = arcpy.NumPyArrayToRaster(j,lowerLeft,cellSize, value_to_nodata=0 )
      newras.save(os.path.join(arcpy.env.workspace,"c" + str(j) + ".tif"))

[DanPatterson]

what is one of the paths? sounds like it is bailing there.  workspace would have to be a folder since you can't save tiffs in a gdb

[NawajishNoman]

Hi Cami,

If your hourly rasters are named sequentially as ih1.tif, ih2.tif, ih3.tif, etc. and you want to create cumulative rainfall rasters CumRain1.tif, CumRain2.tif, CumRain3.tif, etc. where conceptually:

CumRain1.tif = ih1.tif 

CumRain2.tif = CumRain1.tif + ih2.tif

CumRain3.tif = CumRain2.tif + ih3.tif      

and so on...

You can do it using a few simple steps:

1. Copy the first hourly rainfall raster to the first cumulative rainfall raster

2. loop through from 2nd hourly rainfall raster to the last hourly rainfall raster

3. for each hour create a cumulative rainfall raster by adding the hourly rainfall raster to the previous cumulative rainfall raster.

Please see the code snippet below. 

Please note you can also use the Cell Statistics tool to add many rasters using the SUM option.

Hope it helps. Please let me know if you have any questions.

Thanks

Noman

#Code Snippet

import os
import sys
import arcpy
from arcpy.sa import *

#Input, Scratch, and Output Workspace
current_path = os.path.dirname(__file__)
data_folder = os.path.join(current_path, 'data')
out_folder = os.path.join(current_path, 'output')
scratch_folder = os.path.join(current_path, 'scratch')

#Geoprocessing environmnet
arcpy.env.scratchWorkspace = scratch_folder
arcpy.env.workspace = out_folder
arcpy.env.overwriteOutput = True

#Checkout Spatial Analyst extenson
arcpy.CheckOutExtension("spatial")

#Copy 1st hourly rainfall to 1st cumulative rainfall.
firstHourlyRainfallRasterName = "ih"+str(1)+".tif"
firstHourlyRainfallRaster = os.path.join(data_folder, firstHourlyRainfallRasterName)
firstCumulativeRainfallRasterName = "CumRain"+str(1)+".tif"
firstCumulativeRainfallRaster = os.path.join(out_folder, firstCumulativeRainfallRasterName)

arcpy.management.CopyRaster(firstHourlyRainfallRaster, firstCumulativeRainfallRaster)

#Loop through 2nd to last (10th) hourly rainfall.
for i in range(2, 11):

#Specifying input and output raster

PreviousCumulativeRainfallRasterName = "CumRain"+str(i-1)+".tif"
PreviousCumulativeRainfallRaster = os.path.join(out_folder, PreviousCumulativeRainfallRasterName)
hourlyRainfallRasterName = "ih"+str(i)+".tif"
hourlyRainfallRaster = os.path.join(data_folder, hourlyRainfallRasterName)
cumulativeRainfallRasterName = "CumRain"+str(i)+".tif"
cumulativeRainfallRaster = os.path.join(out_folder, cumulativeRainfallRasterName)

#adding hourly rainfall raster with previous cummulative rainfall raster 
CumRain = Raster(PreviousCumulativeRainfallRaster) + Raster(hourlyRainfallRaster)
CumRain.save(cumulativeRainfallRaster)