Run function over different raster layers depending on the grid cell values

Question

I have a Rasterstack X of 92 layers (all days of June, July, August). Which maps the difference between the maximum temperature on that day and a threshold which I have defined.

class      : RasterStack 
dimensions : 201, 464, 93264, 92  (nrow, ncol, ncell, nlayers)
resolution : 0.25, 0.25  (x, y)
extent     : -40.5, 75.5, 25.25, 75.5  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 

To give you an idea of what these layers look like I have plotted the first layer at the bottom of this post: the difference T_max - T_threshold for June 1st. What I want to do now is to take the sum over the layers, but only where X > 0 for at least 3 consecutive days. The result would be just 1 RasterLayer.

I have created the following function to do this. The script does exactly what I want it to when I run it on a vector of numbers.

total <- stack()
sum <- stack()
consec <- stack() 

test <- function(X) {
  for (day in 1:nlayers(X)) {
    if (X[[day]] > 0) {
      consec = consec + 1 
      sum = sum + X[[day]]
    } 
    else if (consec >= 3) { total = total + sum 
    consec = 0
    sum = 0}
    else {consec = 0
    sum = 0 }
  }
  return(total)
}
calc(x = X, fun = test)

This gives me the following error:

Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) : 
  cannot use this function

I think I understand why it gives me this error: calc() usually runs a function on all layers, but I want it to run on different layers for different pixel values.

I have been breaking my head over this for a long time. Does someone know how I can best solve this?

Answer 1

Here is a rewrite of your function to make it useable by calc

test <- function(X) {
    total <- consec <- sum <- 0
    for (day in 1:length(X)) {
        if (X[day] > 0) {
            consec = consec + 1 
            sum = sum + X[day]
        } else if (consec >= 3) { 
            total = total + sum 
            consec = 0
            sum = 0
        } else {
            consec = 0
            sum = 0 
        }
    }
    # note this line that you need if the last day is above 0
    if (consec >= 3) { 
        total = total + sum 
    }
    total
}

test the function

d <- c(10, -1, 5, 5, 0, 3, 3, 3, 3, 0)
test(d)
#[1] 12

Now we can apply it to raster data. First create some example raster data:

library(raster)
r <- brick(ncol=10, nrow=10, nl=100)
set.seed(1)
values(r) <- sample(-5:5, 100*100, replace=TRUE)

Solution

x <- calc(r, test) 

This is a more concise, and probably much faster version of your function

f <- function(x) {
    sum(x[with(rle(x > 0), rep(lengths * values >= 3, lengths))])
}
f(d)
#[1] 12

y <- calc(r, f)

Compare

all(values(x) == values(y))
#[1] TRUE

See related approaches here

If had first understood your objective to sum the values starting on the third day that the threshold had been passed. For reference, it may be helpful to someone (including my future self), I also add the function I wrote for that here

ff <- function(x) {
    y <- x > 0
    z <- ave(y, cumsum(!y), FUN = cumsum)
    sum(x[z>2])
}
ff(d)
#[1] 6