查找 R 中栅格堆栈中第二高值的图层名称

Question

在这个问题（在 R 中查找栅格堆栈上的第二高值）之后，如何找到每个栅格堆栈 xy 坐标的具有第二高值的图层的名称？

我能够使用"which.max（）"函数找到包含最高值的层的名称（层数）：

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)

但是，如何获取名称（图层编号）包含第二高值的栅格？

我尝试了解决方案（如何在 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)

但正如plot(m2)所显示的那样没有成功..

Answer 1

这是一种修改 which.max.na 函数以报告第二高索引的方法。请注意，我添加了 sum(!is.na(x)) == 1，让函数在只有一个非 NA 值时报告0。

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)

我们可以打印前几个值以查看which.max.na和which.second.max.na是否正常工作。

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

对于 RasterStack 中的前六个值，这两个函数似乎都按预期工作。

最后，请注意，如果 RasterStack 中存在连接，则这两个函数可能会出现问题。