raster calc 函数查找最大索引（即最近的图层）满足标准

Question

对于 rasterStack 的每个单元格，我想找到值超过固定阈值的最新层。层按时间顺序堆叠，因此这对应于最大索引。最终我想知道 1）该图层的年份（取自图层名称），以及 2）该年份的值。

我为此编写了一个函数，但得到了错误的结果。我以一种我认为不会改变它的方式修改了函数；我现在得到正确的结果。我的问题是为什么这些函数会产生不同的结果。

设置示例：

library(raster)
library(rasterVis)

###   Example raster stack: 
set.seed(123123)
r1 <- raster(nrows = 10, ncols = 10)
r2 <- r3 <- r4 <- r1
r1[] <- rbinom(ncell(r1), 1, prob = .1)
r2[] <- rbinom(ncell(r1), 1, prob = .1)
r3[] <- rbinom(ncell(r1), 1, prob = .1)
r4[] <- rbinom(ncell(r1), 1, prob = .1)
rs <- stack(r1, r2, r3, r4)
names(rs) <- paste0("yr", 1:4)

这些是我认为相同的功能...是否有理由不使用向量作为中间步骤？

###  Function to find index of last event:
# v1; this is wrong
findLast <- function(x) {
  fire.ind <- ifelse(any(x > 0), which(x > 0), NA)
  max.ind <- max(fire.ind)
}

# v2; this one gives correct answer
findLast2 <- function(x) {
  max.ind <- ifelse(any(x > 0), max(which(x > 0)), NA)
}

testFind <- calc(rs, findLast)
freq(testFind)

testFind2 <- calc(rs, findLast2)
all.equal(testFind, testFind2)

显示示例输入和不同的结果：

# plot:
s2 = stack(rs, testFind, testFind2)
levelplot(s2, pretty = TRUE)

获取我想要的最终层的代码：

###  Most recent year:
nameFromInd <- function(x) {
  yr <- as.integer(gsub(".*(\\d.*).*", "\\1", names(rs)[x]))
}
testYr <- calc(testFind2, nameFromInd)

###  Value in most recent year:
testYrValue <- stackSelect(rs, testFind2)

对我在这里看不到的内容有任何见解吗？我没有尝试过提高速度的替代方案，但欢迎提出任何建议，因为我将在一个非常大的数据集上执行此操作。

sessionInfo()
R version 3.2.4 Revised (2016-03-16 r70336)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)

Answer 1

ifelse 可能有点令人惊讶，因为它返回的值与第一个参数的“形状”相同，长度为 1。所以它只返回which(x > 0) 结果的第一个值。在使用 calc 之前，请务必检查您使用的功能。

x <- c(-1:2, 1:-1)
ifelse(any(x > 0), which(x > 0), NA)
#[1] 3
ifelse(any(x > 0), max(which(x > 0)), NA)
#[1] 5

ifelse 是一个非常复杂的函数，我认为您可以通过以下方式避免它：

r <- calc(rs, function(x) max(which(x > 0)))
y <- stackSelect(rs, r)

（并忽略警告）。或者压制他们：

options(warn=-1)
r <- calc(rs, function(x) max(which(x > 0)))
options(warn=0)
testYrValue <- stackSelect(rs, r)

你也可以像这样结合calc和stackSelect做的事情

f1 <- function(x) {
    if (any(x>0)) {
        i <- max(which(x > 0))
        cbind(i, x[i])
    } else {
        cbind(NA, NA)
    }
}

rr1 <- calc(rs, f1)

或者快捷方式：

f2 <- function(x) {
    i <- max(which(x > 0))
    cbind(i, x[i])
}

rr2 <- calc(rs, f2)

或者

f3 <- function(x) {
    i <- max(which(x > 0))
    z <- cbind(i, x[i])
    z[!is.finite(z)] <- NA
    z
}

rr3 <- calc(rs, f3)