1

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

Here's a reproducible code example generating rasters and making a stack.

# Generate stack of random values
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 analyze this stack to determine which layer is the first 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 matrix with columns representing each layer, this may make things easier.

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

A somewhat similar question has been previously posted as Counting layer number of raster stack with value under some conditions and write it into output layer?, though that's counting the number of layers that meet a condition.

  • Are all your conditions functions of a single pixel-stack? Write a function that returns the index of the element that satisfies a condition given a vector of pixels and then use calc to compute that over your stack and then find the minimum of that? – Spacedman Jul 11 at 18:59
  • Yes, everything references a single pixel at a time. For now I've: Run rasterToPoints to create a data frame of raster values Looped through that with a while loop to find the first entry that satisfies a condition: test<-as.data.frame(rasterToPoints(s)) for (j in 1:nrow(test)){ i=7 while ((test[j,i] + test[j,i-1] + test[j,i-2] + test[j,i-3] + test[j,i-4]) < 2) { i <- i+1 } if (i==(ncol(test)-1)) i<-0 test$first[j]<-i } Any help vectorizing that or running that on the raster stack itself would be helpful though. – Scott Z Jul 11 at 21:11
1

If you can write a function that works on a single vector of values then you can use calc to run that over a raster. For example here's a function that returns the position of the first run of three ones in a vector, or NA if no run of three ones exists:

threeones <- function(x){
    r = rle(x) # use rle to count runs of things
    got31s = r$lengths >=3 & r$values == 1 # any runs of three 1s?
    n = NA
    if(any(got31s)){
        pos = min(which(got31s)) # get the first run of three
        n = sum(r$lengths[1:(pos)]) - r$lengths[pos]+1 # get pos of first 1
    }
    return(n)
}

You can then test this with some test values:

> tv = list(c(1,1,1,0,2,3), c(0,1,2,1,2,3), c(0,1,1,1,0,0), c(0,0,0,0,1,1), c(0,1,1,1,1,2,1,1,1))
> lapply(tv,threeones)
[[1]]
[1] 1

[[2]]
[1] NA

[[3]]
[1] 2

[[4]]
[1] NA

[[5]]
[1] 2

Which all looks good. Now run it on your test stack:

> s31 = calc(s,threeones)

and check some locations. At position 6,6 there's a sequence of four ones starting at layer 3:

> s[6,6]
     layer.1 layer.2 layer.3 layer.4 layer.5 layer.6 layer.7
[1,]       3       0       1       1       1       1       3

and s31 at that point is....

> s31[6,6]
[1] 3
> 

There's probably a better implementation of my threeones function (its a bit fiddly and there's always nuisance edge casess) but I think the important thing is to show you how to run over pixels using calc.

Note that threeones isnt vectorised, and doesn't have to be for calc to work.

  • Very helpful, thank you! I'm not very good with writing functions so this definitely helps. I worked out a way using a for and while loop - this is probably much faster than that. – Scott Z Jul 14 at 16:22

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.