0
votes

I want to analyze a raster stack and determine the first layer in the stack that meets a condition for every pixel.

# Generate stack of random values
library(raster)

r1 <- r2 <- r3 <- r4 <-r5 <- r6 <- r7 <- raster(ncol=10, nrow=10)
r1[] <- rpois(ncell(r1), 1)
r2[] <- rpois(ncell(r2), 1)
r3[] <- rpois(ncell(r3), 1)
r4[] <- rpois(ncell(r4), 1)
r5[] <- rpois(ncell(r5), 1)
r6[] <- rpois(ncell(r6), 1)
r7[] <- rpois(ncell(r7), 1)

# stack them
s <- stack(r1,r2,r3,r4,r5,r6,r7)

I want to take this stack and determine the first layer to meet a condition like: "Value = 1 for 3 consecutive layers" or "2 of last 3 layers have value = 1", for example.

Can you think of any function to do this?

If I run rasterToPoints I can convert the stack into a dataframe with columns representing each layer value. This may make things easier.

values<-as.data.frame(rasterToPoints(s))
values<-values[,3:ncol(values)]  #eliminates two columns at the beginning
values$first<-NA

At this point I think I would want to loop through and fill in the values$first column with the first column meeting my condition.

Once I determine the correct layer for each pixel, I then want an output raster with that correct layer number as the value for each pixel.

I'd appreciate any help!

Somewhat similar question posted here, though that's counting the number of layers that meet a condition.

1

1 Answers

0
votes

You can write a function that can do this for a vector and then use it in raster::calc.

For example, to find the first layer that has a value of 3

f <- function(x) which(x == 3)[1] 

(and handles the case that none of the values is 3 well, by returning NA)

For example

library(raster)
s <- stack(system.file("external/rlogo.grd", package="raster"))
f <- function(x) which(x == 255)[1] 
x <- calc(s, f)
plot(x)