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 函数 return 会出现错误消息。

考虑以下示例。

# 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 时，第二个示例只是 return 出错了。这与将 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，此函数将 return 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

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

R raster::calc calculating quantile with na.rm = FALSE

r

quantile

na

r-raster

更新