using function inside raster::calc with loop

Question

I have a raster stack with 12 raster layers (NDVI time series). My goal is to count PIXELWISE the number of managing events (grazing, mowing) and receive a raster layer that contains for every pixel the number of events.

For demonstration a random raster stack is produced:

r1 <- raster(nrows = 1, ncols = 1, res = 0.5, xmn = -1.5, xmx = 1.5, ymn = -1.5, ymx = 1.5, vals = 0.3)
rr <- lapply(1:12, function(i) setValues(r1,runif(ncell(r1))))
rr<-(stack(rr))

These events are indicated by:

1. a fall (>= 0.1) of the NDVI value and
2. a rise (< -0.05 & > - 0.32) of the NDVI value after the event

My idea is to loop over the number of rasterlayers - 2 (last two dates can't be used),to insert a number into a empty string IF a event occurs and return the length of the vector. Putting the function into raster::calc should do this for every pixel.

 count_events <- function(x){
 events <- c()
 for (i in dim(x)[3]-2) {
 if(isTRUE((x[[i]]-x[[i+1]]) >= 0.1 &
 (x[[i+1]]-x[[i+2]]) < - 0.05 &
 (x[[i+1]]-x[[i+2]]) > -0.32)){
 events <- c(events,1)
 }
 }
 return(length(events))
}
raster::calc(rr,count_events)
class : RasterLayer 
dimensions : 6, 6, 36 (nrow, ncol, ncell)
resolution : 0.5, 0.5 (x, y)
extent : -1.5, 1.5, -1.5, 1.5 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs 
source : memory
names : layer 
values : 0, 0 (min, max)

The resulting raster has always the min/ max values 0/0. I'm sure that there are these events, and the random rasterstacks gives the same result. I think there is mistake in the function using the brackets. It would be very kind if somebody could help me to find the error.

Answer 1

I would use lag function from dplyr package:

Your data:

library(raster)
r1 <- raster(nrows = 1, ncols = 1, res = 0.5, xmn = -1.5, xmx = 1.5, ymn = -1.5, ymx = 1.5, vals = 0.3)
rr <- lapply(1:12, function(i) setValues(r1,runif(ncell(r1))))
rr<-(stack(rr))

The custom function with the description you provided:

library(dplyr)
custom_fun <- function(x){
 x <- x[1:10] # only first 10 NDVI values
 diff_ <- x-lag(x)
 events_ <- sum(diff_ >= 0.1 | (diff_< -0.05 & diff_ > -0.32), na.rm = T)
}

In this case, always the first value is of diff_ is NA since it doesn't have any value to compare. Then the result:

plot(calc(rr,custom_fun))

enter image description here

• talocodat: Thanks so far! The lag function was what I was looking for! But how can I include the temporal dimension into the function, so that the conditions are met for two values next to each other?
• aldo_tapia: That sounds like a different question. Please ask a new question and explain your problem with examples