CODEWITHSUNDEEP
rrastercumulative-summulticoreterra
tabumis
tabumis

Reputation: 106

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

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.

Upvotes: 0

Views: 399

Answers (3)

pdbentley
pdbentley

Reputation: 440

Here is an alternative approach inspired by here Cumulative sum that resets when 0 is encountered ; this approach will reset the counter to zero

library(terra)

v <-  c(1,1,1,1,0,0,0,0,1,1,1,0,0,0,1,1)

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

ff = function(x)
{
  cs = cumsum(x)
  cs - cummax((x == 0) * cs)
}

r2 <- app(r, ff)
plot(r2)

Upvotes: 1

Robert Hijmans
Robert Hijmans

Reputation: 47146

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
        i <- i[1]
        x[i] <- x[i] - v[i]
    }
    v
}

cumsumreset(v)
#[1]  1  3  0 35 45 55

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,    35,    45,    55 
#max values  :     1,     3,     0,    35,    45,    55

Upvotes: 1

tabumis
tabumis

Reputation: 106

UPDATED. Apparently the function I posted before was incorrect. Below is a solution which benefits from parallelization of the operation via cores argument in the terra::app function. Though I assume this is not the most optimal suggestion.

library(terra)
v <-  c(1, 2, -40, 35, 10, 10)
r <- rast(ncol=1, nrow=1, nlyr=6, vals=v)
r

cumsumreset <- function(x) {
  for (i in 2:length(x)) {
    if (x[i-1] < 0) x[i-1] <- 0
    x[i] <- x[i] + x[i-1]
    if (x[i] < 0) x[i] <- 0
  }
  x
}


a <- app(r, cumsumreset, cores=5)

a

P.S. My appologies and thanks to @RobertHijmans

Upvotes: 1

Related Questions