0

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.

1 Answer 1

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

2
  • 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?
    – talocodat
    May 10 at 8:47
  • That sounds like a different question. Please ask a new question and explain your problem with examples
    – aldo_tapia
    May 10 at 12:44

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.