-1

2012 年 9 月 17 日更新

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

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

尺寸：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))
}
4

2 回答 2

5

这是一个解决方案的猜测。这不是因为 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))
于 2012-09-17T04:07:56.073 回答
0

所以，这是我发现的最后一个问题和我的解决方案。

我遇到的问题 which.max() 是它如何处理所有 NA 的向量

>which.max(c(NA,NA,NA))
integer(0)
>

当 stackApply() 函数尝试将此值写入新的 RasterLayer 时，它会失败。在函数应该返回 NA 的地方，它返回长度 = 0 的整数（0）。

解决方案（我的解决方案）是为 which.max() 编写一个包装器

which.max.na <- function(x, ...){
   max_idx <- which.max(x)
   ifelse(length(max_idx)==0,return(NA),return(max_idx))
}

这个，在我原来的 RasterStack 上实现的效果很好。感谢所有人的帮助，如果您有这些解决方案，请提出替代方案！

谢谢！凯尔

于 2012-09-18T04:05:28.100 回答