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

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? • Can you show the function "does exactly what I want it to when I run it on a vector of numbers." with an example? Given a working function, this would be easy. Dec 11, 2020 at 17:32
• You are writing your function like you are operating on a raster stack. It should be written like you are operating on a vector. The `calc` and `overlay` function`, internally return a vector of values associated with each cell. For example, you are using `nlayers(X)` and a double bracket for indexing. It would be `1:length(X)` and `X[day]`. I would create a dummy vector eg., `y<-runif(92,0,10)` to test your function. I would also suggest looking at the `which`and `rle` functions to optimize your code. Dec 11, 2020 at 17:32
• @JeffreyEvans, thank you a lot for this. My original function had single brackets and length() but I thought I had to adapt it to work on the RasterStack. Very good tip! Dec 12, 2020 at 19:13

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)
# 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)
# 12

y <- calc(r, f)
``````

Compare

``````all(values(x) == values(y))
# 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)
# 6
``````
• Thank you so much for this! Works perfectly! I had tried rle() before but I thought it only worked on values which are the same, not all values above a certain threshold. Dec 12, 2020 at 19:16