R raster::calc 使用 na.rm = FALSE 计算分位数

Question

我使用 raster::calc 计算跨不同层的每个单元格的分位数，但我不理解 na.rm = FALSE 时的行为，如下例所示。

让我们创建一个示例栅格并从随机单元格中删除 5 个值。

library(raster)

r <- raster::raster(nrow = 2, ncol = 2)
r[] <- 1:4

s <- raster::stack(r, r*2, r * 3, r * 4, r * 5)
s[]

set.seed(1)
s[][sample(1:4, 1), sample(1:5, 1)] <- NA
s[][sample(1:4, 1), sample(1:5, 1)] <- NA
s[][sample(1:4, 1), sample(1:5, 1)] <- NA
s[][sample(1:4, 1), sample(1:5, 1)] <- NA
s[][sample(1:4, 1), sample(1:5, 1)] <- NA
s[]

如果我删除 NA，下面的代码就可以工作！

fun <- function(x) {quantile(x, probs = 0.50, na.rm = TRUE)}
p <- raster::calc(s, fun)
p[]

但是，如果我想排除至少有一个 NA 的单元格，则代码不起作用！

fun <- function(x) {quantile(x, probs = 0.50, na.rm = FALSE)}
p <- raster::calc(s, fun)

我期待一个包含 4 个 NA 的向量，但上面的代码却抛出了这个错误：

Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) : 
  cannot use this function

谁能帮我理解为什么会这样？我应该怎么做才能得到我所期望的行为？

Answer 1

我认为错误信息很简单，可能与光栅包无关。基本思想是，如果将 quantile 函数应用于具有任何 NA 的值，quantile 函数将返回错误消息。

考虑以下示例。

# Set na.rm = TRUE
quantile(c(1, NA, 3, 4), probs = 0.50, na.rm = TRUE)

50% 
  3

# Set na.rm = FALSE
quantile(c(1, NA, 3, 4), probs = 0.50, na.rm = FALSE)

Error in quantile.default(c(1, NA, 3, 4), probs = 0.5, na.rm = FALSE) :
   missing values and NaN's not allowed if 'na.rm' is FALSE

当设置na.rm = FALSE 时，第二个例子刚刚返回了一个错误。将quantile 应用于栅格时也是如此。 na.rm 必须为真。

更新

为了说明如何在某些单元格为 NA 时应用分位数函数，我稍微修改了 OP 中的示例数据集。

s <- raster::stack(r, r*2, r * 3, r * 4, r * 5)
s[]

set.seed(1)
s[][sample(1:4, 1), sample(1:5, 1)] <- NA
s[][sample(1:4, 1), sample(1:5, 1)] <- NA
s[][sample(1:4, 1), sample(1:5, 1)] <- NA
s[]

     layer.1 layer.2 layer.3 layer.4 layer.5
[1,]       1       2       3       4      NA
[2,]       2      NA       6       8      10
[3,]       3       6       9      12      NA
[4,]       4       8      12      16      20

查看最后一行，其中没有任何NA。

然后我们可以创建一个函数。如果某个位置的任何层有任何NA，此函数将返回NA。否则，它将计算分位数。

# Design a function
quantile_fun <- function(x, probs = 0.50){
  if (anyNA(x)){
    return(NA)
  } else {
    return(quantile(x, probs = probs))
  }
}

之后，我们可以使用calc来应用这个函数

p <- raster::calc(s, quantile_fun) 
p[]
[1] NA NA NA 12