R terra:: cumulative sum reset function across layers of a stack-CodePudding

I have a stack with large number of layers, and I want to run a pixel-wise cumulative sum function across the layers with reset when sum becomes <0. Ive compiled the following code:

# given that "A" is an input stack with lengthy number of layers

B<-A 

for(i in 2:nlyr(B)){
  B[[i]]<-ifel((A[[i]] B[[i-1]])>0, A[i]] B[[i-1]], 0 )
  
}

, which works well but I have a large stack of rasters (in terms of nlyrs). Im wondering if there is a faster way to implement this task.

CodePudding user response：

Well, looks like the following code does this job with a much faster speed:

B <- terra::app(A, function(x) {
  cumsum_x <- cumsum(x)
  cumsum_x[cumsum_x < 0] <- 0
  return(cumsum_x)
})

CodePudding user response：

Example data

library(terra)
v <-  c(1, 2, -40, 35, 10, 10)
cumsum(v)
#[1]   1   3 -37  -2   8  18

r <- rast(ncol=1, nrow=1, nlyr=6, vals=v)

If you want to set the negative values to zero, you can do that like this

x <- cumsum(r)
x <- ifel(x < 0, 0, x)
x
#class       : SpatRaster 
#dimensions  : 1, 1, 6  (nrow, ncol, nlyr)
#resolution  : 360, 180  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 
#source(s)   : memory
#names       : lyr.1, lyr.2, lyr.3, lyr.4, lyr.5, lyr.6 
#min values  :     1,     3,     0,     0,     8,    18 
#max values  :     1,     3,     0,     0,     8,    18

But that does not reset the cumulation at zero. If that is what you want you need something more complex. Perhaps like this:

cumsumreset <- function(x) {
    if (x[1] < 0) x[1] <- 0
    while(TRUE) {
        v <- cumsum(x)
        i <- which(v < 0)
        if (length(i) == 0) break
        x[i] <- x[i] - v[i]
    }
    v
}

a <- app(r, cumsumreset)

a
#class       : SpatRaster 
#dimensions  : 1, 1, 6  (nrow, ncol, nlyr)
#resolution  : 360, 180  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 
#source(s)   : memory
#names       : lyr.1, lyr.2, lyr.3, lyr.4, lyr.5, lyr.6 
#min values  :     1,     3,     0,    37,    47,    57 
#max values  :     1,     3,     0,    37,    47,    57