CODEWITHSUNDEEP
rrasterr-raster
Marc
Marc

Reputation: 661

Find layer names for the second highest values in a raster stack in R

Following this question (Find second highest value on a raster stack in R), how can one find, for each raster stack xy coordinate, the name of the layer holding the second highest value?

I'm able to find the name (layer number) of the layer containing the highest value with a "which.max()" function:

set.seed(123)
require(raster)
r1 <- raster(nrows = 10, ncols = 10)
r2 <- r3 <- r4 <- r1
r1[] <- runif(ncell(r1))
r2[] <- runif(ncell(r1)) + 0.2
r3[] <- runif(ncell(r1)) - 0.2
r4[] <- runif(ncell(r1))
rs <- stack(r1, r2, r3, r4)

which.max.na <- function(x, ...) ifelse(length(x) == sum(is.na(x)), 0, which.max(x))

m1 <- calc(rs, which.max.na)

plot(m1)

However, how can I obtain a raster with the names (layer numbers) containing the second highest values?

I tried the solution in (How to find second highest value and corresponding layer name in a raster stack in R):

m2 <- calc(rs, fun=function(x, na.rm) x[order(x, decreasing=T)[2]]) & calc(rs, fun=function(x, na.rm) order(x, decreasing=T)[2])

plot(m2)

but without success as plot(m2) shows..

Upvotes: 0

Views: 358

Answers (1)

www
www

Reputation: 39154

Here is an approach modifying the which.max.na function to report the second highest index. Notice that I added sum(!is.na(x)) == 1 to let the function report 0 when there is only one non-NA value.

which.second.max.na <- function(x, ...) 
  ifelse(length(x) == sum(is.na(x)) | sum(!is.na(x)) == 1, 0, 
         which.max(`[<-`(x, which.max(x), NA)))

m2 <- calc(rs, which.second.max.na)

We can print the first few values to see if which.max.na and which.second.max.na are working.

head(values(m1))
[1] 2 1 4 2 1 2

head(values(m2))
[1] 4 3 2 1 2 3

head(values(rs))
       layer.1   layer.2    layer.3     layer.4
[1,] 0.2875775 0.7999890 0.03872603 0.784575267
[2,] 0.7883051 0.5328235 0.76235894 0.009429905
[3,] 0.4089769 0.6886130 0.40136573 0.779065883
[4,] 0.8830174 1.1544738 0.31502973 0.729390652
[5,] 0.9404673 0.6829024 0.20257334 0.630131853
[6,] 0.0455565 1.0903502 0.68024654 0.480910830

It seems like for the first six values from the RasterStack, both functions work as expected.

Finally, be careful that if there are ties in the RasterStack, these two functions could be problematic.

Upvotes: 1

Related Questions