在每个单元格的 R RasterStack 上实现 which.max()答案

【问题标题】：Implementing which.max() on an R RasterStack for each cell在每个单元格的 R RasterStack 上实现 which.max()
【发布时间】：2012-09-17 02:37:04
【问题描述】：

2012 年 9 月 17 日更新

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

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

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

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

# 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))
}

【问题讨论】：

我不认为你关于为什么它不起作用的结论是正确的。 ?which.max 说：“丢失和 NaN 值被丢弃。”您能否使用 dput() 提供该对象的一小部分？
哇。来吧。那是一个网页，而不是一个 FTP 服务器。您应该提供用于创建此对象的代码！我有一个解决方案，但无法用这种模糊程度进行测试。
请提供reproducible example，我们可以快速运行以测试您的问题的问题和解决方案。
所有，今晚我将尝试添加更多信息。我将把我尝试过的不同方法整合在一起，并产生错误。感谢您迄今为止的帮助。

标签： r max apply raster

【解决方案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))

【讨论】：

DWin。谢谢你的工作。我以为我提供了足够的信息。一个指向数据的链接（简单到可以下载）和用于加载数据的所有代码。您只需更改数据搜索路径即可。但是，今晚我将试一试您的解决方案并分享结果。我知道 calc(PRISM_stack,which.max) 对我不起作用。我会为你重现错误。
DWin。一个旁注。我遇到的一个问题是，当 RasterStack 的所有 12 层都是“NA”时，这导致了 which.max() 的字符（0）返回，这可能会导致我的问题。再次，我将一起获得更多信息晚上。谢谢！
我浏览到的那个页面有超过 48 个链接和很多年，我不清楚使用了哪些文件。所以你期望我们重复你的“手工”，甚至没有给出明确的指示如何去做。我期待你会构建像 sapply(filenames, function(x) x <- read.table(file=paste0(url, x)) ) 这样的代码
我已更新问题将重现我的错误所需的所有信息。
与我提交的 which.max2() 函数相比，您的 which.max2() 函数有什么好处（效率？）。我仍在学习 R 语言，希望能在这里获得一些见解。

【解决方案2】：

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

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

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

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

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

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

这个，在我原来的 RasterStack 中实现的工作正常。感谢大家的帮助，如果您有此解决方案的替代方案，请提出替代方案！

谢谢！凯尔

【讨论】：

这基本上是我半小时前发布的解决方案。