0

I am working with a rasterbrick "a" with thousands of layers, closer description is not necessary for my problem. I am using following function to create a rasterlayer of the total amount of runs of at least 5 days with values greater than 1 (one layer in brick is one day):

 indices<-rep(1:69,each=90)
 ff<-function(x,na.rm=TRUE){
 y<-x > 1
 n<- ave(y,cumsum(y == 0), FUN = cumsum)
 sum(n==5)
 }
 Y<-stackApply(a,indices,fun=ff)

This works great, I tested that. In a similar manner, I wrote new function:

 fff<-function(x,na.rm = TRUE){
 y <- x > 1
 n <- ave(y, cumsum(y == 0), FUN = cumsum)
 mean(n[n >= 5])
 }
 X<-stackApply(a,indices,fun=fff)

Using this function, I wanted to create a rasterlayer of average lengths of those runs greater than 5 days. It seems reasonable and fine, but it does not work correctly. For example, when there is a run of 6 days (satisfying my criterion of value>1), it counts two runs, one of 5 and another one of six, and thus the average is 5,5 instead of 6. I am not sure how to adjust my function fff. If there is a way to do it, it would be great, otherwise I would be greatful if anyone shares another way how to calculate means of those runs. Thanks!

  • Hi Adam Welcome to SO. It would be great if you provide a sample of your `a` and wha is your desired expected output. Consider reviewing https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example – Chriss Paul Mar 22 '21 at 17:34
  • My brick: class : RasterBrick dimensions : 201, 241, 48441, 6210 (nrow, ncol, ncell, nlayers) resolution : 0.25, 0.25 (x, y) extent : -15.125, 45.125, 24.875, 75.125 (xmin, xmax, ymin, ymax) Values range from some -2 to 2. My output should be a Rasterlayer. Values at each pixel represent the average lenght of runs of at least 5 days with values greater than 1. Layers in brick are days, so if value at one pixel exceeds 1 for at least 5 layers in a row, it is detected as a run. My first function ff sucessfully calculates layer of total amount of those runs. – Adam Dragula Mar 22 '21 at 17:58
  • This is not what comments are for. Please edit your question instead. Do not describe the RasterBrick but generate it with code so that we can run your code – Robert Hijmans Mar 22 '21 at 18:40
  • @RobertHijmans, You answered my question called finding out lengths of runs...and suggested this function and it really looked fine, but after test, I found out that it doesnt work. – Adam Dragula Mar 22 '21 at 18:54
  • Please edit your question and show that it does not work. In your question, you can add a hyperlink to a previous question. This website, and the raster documentation, has 1000s of examples of how to write a minimal, self-contained, reproducible example. – Robert Hijmans Mar 22 '21 at 20:36

1 Answers1

0

In the future, please include a minimal, reproducible, self-contained example. Do not describe the behavior of your code, but show it. Also, be very clear aobut the question. As-is it is hard to see that your question is not about raster data at all, as you are looking for a function that works on any numeric vector (that you may then apply to raster data).

You are looking for function that finds local maxima larger than 5, in the cumulated sum of neighbors that are > 1; and then average these local maxima.

You have this function

ff<-function(x,na.rm=TRUE){
   y<-x > 1
   n <- ave(y,cumsum(y == 0), FUN = cumsum)
   sum(n==5)
}

Example data

x <- c(-1:10, -1:3, -1:6)
x
#  [1] -1  0  1  2  3  4  5  6  7  8  9 10 -1  0  1  2  3 -1  0  1  2  3  4  5  6

ff(x)
# [1] 2

(two local maxima that are larger than 5)

To write the function you want we can start with what we have

y <-x > 1
n <- ave(y,cumsum(y == 0), FUN = cumsum)
n        
# [1] 0 0 0 1 2 3 4 5 6 7 8 9 0 0 0 1 2 0 0 0 1 2 3 4 5

In this case, you need to find the numbers 9 and 5. You can start with

n[n<5] <- 0
n
# [1] 0 0 0 0 0 0 0 5 6 7 8 9 0 0 0 0 0 0 0 0 0 0 0 0 5

And now we can use diff to find the local maxima. These are the values for which the difference with the previous value is negative. Note the zero added to n to consider the last element of the vector.

i <- which(diff(c(n, 0)) < 0)
i
# [1] 12 25
n[i]
# [1] 9 5

Such that we can put the above together in a function like this

f <- function(x) {
   y <- x > 1
   n <- ave(y,cumsum(y == 0), FUN = cumsum)
   n[n<5] <- 0
   i <- which(diff(c(n, 0)) < 0)
   mean(n[i])
}

f(x)
# [1] 7

If you have NAs you may do

f <- function(x) {
   y <- x > 1
   y[is.na(y)] <- FALSE
   n <- ave(y,cumsum(y == 0), FUN = cumsum)
   n[n<5] <- 0
   i <- which(diff(c(n, 0)) < 0)
   mean(n[i])
}
Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63