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?
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')