在 R RasterStack 上为每个单元实现 which.max()

Question

更新 09/17/2012

这是一段使用自包含数据重现我的问题的代码：

请记住，我拥有的实际数据维度很大......

尺寸：3105、7025、21812625、12（nrow、ncol、ncell、nlayers）

我需要的是每一行的最大值的索引，层上的列。所有 NA 都应该返回 NA 并且多个最大副本应该返回第一个最大索引（或其他东西，必须一致）

# Create a test RasterStack

require(raster)

a <- raster(matrix(c(11,11,11,
                     NA,11,11,
                     11,11,13),nrow=3))

b <- raster(matrix(c(12,12,12,
                     NA,12,12,
                     40,12,13),nrow=3))

c <- raster(matrix(c(13,9,13,
                     NA,13,13,
                     13,NA,13),nrow=3))

d <- raster(matrix(c(10,10,10,
                     NA,10,10,
                     10,10,10),nrow=3))

corr_max <- raster(matrix(c(13,12,13,
                            NA,13,13,
                            40,12,13),nrow=3))

stack <- stack(a,b,c,d)


which.max2 <- function(x, ...)which.max(x)

# stackApply method
max_v_sApp <- stackApply(stack,rep(1,4),which.max2,na.rm=NULL)

# calc method
max_v_calc <- calc(stack,which.max)

希望这提供了足够的信息。

更新：

这可能有效......现在测试：

which.max2 <- function(x, ...){
  max_idx <- which.max(x)   # Get the max
  ifelse(length(max_idx)==0,return(NA),return(max_idx))
}

Answer 1

这是对解决方案的猜测。这不是因为 which.max 不“支持” na.rm 参数，只是因为它已经假设了它并且只有“有空间”用于数据参数。在从帮助页面中提取的一个小测试用例上进行了测试，但没有在您的数据上进行测试。您可以使用以下任一方法：

require(raster)
 which.max2 <- function(x, ...) which.max(x)           # helper function to absorb "na.rm"
 wsa <- stackApply(PRISM_stack, rep(1,12), fun=which.max2, na.rm=NULL)

显然这种方法不需要帮助函数来剥离 na.rm：

calc(PRISM_stack, which.max)

由于单元格中所有 NA 的新问题，这两种方法似乎都成功了：

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

就像这样：

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