CODEWITHSUNDEEP
rasterterra
jmutua
jmutua

Reputation: 293

Calculating cumulative sum of rasters in a stack

I am calculating the cumulative sum of rasters in a stack using methods provided by both terra and raster packages

First, I try to do it step by step as follows:

library(terra)
library(raster)

r <- rast(ncols=9, nrows=9)
values(r) <- 1:ncell(r)
s <- c(r, r, r, r, r, r)
s <- s * 1:6
x <- s[[1]]+s[[2]]
x <- x+s[[3]]
x <- x+s[[4]]
x <- x+s[[5]]
rr1 <- x+s[[6]]

When using terra, the output seem not correct

rr2 <- rapp(s, s[[1]], nlyr(s), function(i) max(cumsum(i)))

However, when using raster, the values seem OK.

rs<- raster::stack(s)
rs <- calc(rs, fun = cumsum)
rr3 <- rs[[6]]

What could be the problem?

Upvotes: 0

Views: 954

Answers (1)

Robert Hijmans
Robert Hijmans

Reputation: 47146

Your example data

library(terra)
#terra 1.6.30    
r <- rast(ncols=9, nrows=9)
values(r) <- 1:ncell(r)
s <- c(r, r, r, r, r, r) * 1:6

With terra 1.6-30 you can use cumsum (there is a bug in the CRAN version that prevents this from working)

cumsum(s)
#class       : SpatRaster 
#dimensions  : 9, 9, 6  (nrow, ncol, nlyr)
#resolution  : 40, 20  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 
#source(s)   : memory
#names       : lyr.1, lyr.1, lyr.1, lyr.1, lyr.1, lyr.1 
#min values  :     1,     3,     6,    10,    15,    21 
#max values  :    81,   243,   486,   810,  1215,  1701

A work-around is

app(s, cumsum)
#class       : SpatRaster 
#dimensions  : 9, 9, 6  (nrow, ncol, nlyr)
#resolution  : 40, 20  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 
#source(s)   : memory
#names       : lyr.1, lyr.1, lyr.1, lyr.1, lyr.1, lyr.1 
#min values  :     1,     3,     6,    10,    15,    21 
#max values  :    81,   243,   486,   810,  1215,  1701

But in case you only want the last layer, you can of course do

sum(s)
#class       : SpatRaster 
#dimensions  : 9, 9, 1  (nrow, ncol, nlyr)
#resolution  : 40, 20  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 
#source(s)   : memory
#name        :  sum 
#min value   :   21 
#max value   : 1701

rapp does not apply here, but you could use your formula with app

app(s, \(i) max(cumsum(i)))

You can install terra version 1.6-30 with:

install.packages('terra', repos='https://rspatial.r-universe.dev')

Upvotes: 1

Related Questions