Calculate time since event in a raster time series in R

Question

I have a time series of binomial data and I want to know the number of time steps since the last event (1). I am just not sure how to do this in R.
Here is an example data set. For each cell in the raster, I want to know the number of timesteps since a 1 was recorded.
r <- raster(nrow=100, ncol=100)
r[] = rbinom(10000,1,.1)

for (i in 1:9){
  tmp = raster(nrow=100, ncol=100)
  tmp[] = rbinom(10000,1,.1)
  r = stack(r,tmp)
}

r = brick(r)

So for each cell of the raster brick, you would have a series (e.g. c(0,0,0,0,1,0,0,0,1,0)). In this case the time to the event would by 4. I just can't figure out how to do this in a function that can be applied to the raster brick.

Jeffrey Evans · Accepted Answer

Here is an alternative using rle.
Write function to calculate run length to specified value. If the specified run value (y) is the first value it returns a 0 and if it is missing then a NA is returned.
time.to.event <- function(x, y=1, dir=c("LR", "RL")) {
    if(dir[1] == "RL") x <- rev(x)
      x.idx <- rle(x)$values
      e.idx <- rle(x)$lengths 
    if(x[1] == y) {
      e <- 0
    } else {
      e <- e.idx[which(x.idx == y)-1][1]
    }   
  return( e )
}

Test the function on your example, left-to-right (LR is default) and right-to-left (dir="RL").
time.to.event(c(0,0,0,0,1,0,0,0,1,0))
time.to.event(c(0,0,0,0,1,0,0,0,1,0),dir="RL")

Now, apply it to a raster stack example
library(raster)
r <- do.call(raster::stack, replicate(20,raster::raster(matrix(sample(
             c(0,1), 1000, replace=TRUE), 100, 100))))
             
( t2e <- calc(r, fun=time.to.event) )

Spacedman · Answer

You can use calc with a function that finds the first one. The only complication is coping with not having any ones in the stack at any point, which I workaround here by testing and returning -1.
> fone= calc(r, fun=function(x){W=which(x==1); if(length(W)==0){return(-1)};return(min(W))})

To check, the first pixel out is:
> fone[1,1]
[1] 7

which comes from:
> r[1,1]
     layer.1.1.1 layer.2.1.1 layer.1.2.1 layer.2.2.1 layer.1.1.2 layer.2.1.2
[1,]           0           0           0           0           0           0
     layer.1.2.2 layer.2.2.2 layer.1 layer.2
[1,]           1           0       0       0
>

which has a 1 in the 7th position.
Couple more tests:
> fone[50,50]
[1] -1

So there shouldn't be any 1s in the raster at (50,50):
> r[50,50]
     layer.1.1.1 layer.2.1.1 layer.1.2.1 layer.2.2.1 layer.1.1.2 layer.2.1.2
[1,]           0           0           0           0           0           0
     layer.1.2.2 layer.2.2.2 layer.1 layer.2
[1,]           0           0       0       0

and (20,21) has more than 1:
> r[20,21]
     layer.1.1.1 layer.2.1.1 layer.1.2.1 layer.2.2.1 layer.1.1.2 layer.2.1.2
[1,]           1           0           0           0           0           0
     layer.1.2.2 layer.2.2.2 layer.1 layer.2
[1,]           0           0       1       0

and we catch the first one:
> fone[20,21]
[1] 1

There's probably a neater expression for the function that does it without an "if" but my brain is not up to it right now. Will edit if I think of it.
Thought.
Stick a "1" on the end of the vector of values being tested as a "guard value", then you know that which will definitely find something, and if it returns the number of layers in r plus 1, its found the guard value:
firstone = function(x){min(which(c(x,1)==1))}
fone = calc(r, firstone)

and if you want to set the pixels with no ones to something:
fone[fone == nlayers(r)+1] = NA

or whatever.
You say that
c(0,0,0,0,1,0,0,0,1,0))

should return 4, but my code will return 5. Subtract one from the value in my function.

Calculate time since event in a raster time series in R

2 Answers

Add your own answers!

Ask a Question