0

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?

enter image description here

3
  • 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
  • 2
    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 whichand 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!
    – Charlotte
    Dec 12, 2020 at 19:13

1 Answer 1

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
1
  • 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.
    – Charlotte
    Dec 12, 2020 at 19:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.