R中raster::quantile计算NA分布不同的栅格堆栈返回错误值问题

问题描述 #

我尝试在R中使用raster::quantile计算栅格堆栈（raster stack）的分位数，单张栅格计算结果正常，但当堆栈内不同栅格的NA像元位置不一致时，得到的分位数结果存在错误。

问题复现 #

library(raster)
#> Loading required package: sp

# 构建栅格堆栈
r <- raster(system.file("external/test.grd", package="raster"))
r1 <- setValues(r, sample(100:2000, ncell(r), replace = TRUE))
s <- stack(r, r1)

plot(s)

# 直接对栅格堆栈计算分位数
quantile(s)
#>              0%      25%       50%       75%     100%
#> test.1 138.7071 293.9575  371.9001  501.0102 1736.058
#> test.2 100.0000 596.0000 1062.0000 1530.0000 2000.000

# 对堆栈内的栅格单独计算分位数
quantile(s[[1]])
#>        0%       25%       50%       75%      100% 
#>  138.7071  293.9575  371.9001  501.0102 1736.0580

quantile(s[[2]])
#>   0%  25%  50%  75% 100% 
#>  100  591 1071 1530 2000

从结果可见，对栅格堆栈整体调用quantile得到的第二张栅格的25%、50%分位数，与单独计算该栅格分位数的结果不一致。若堆栈内所有栅格的NA值位置完全相同，则不会出现该问题，示例如下：

quantile(stack(r1, sqrt(r1)))
#>        0%       25%        50%        75%       100%
#> test  100 585.00000 1045.50000 1537.00000 2000.00000
#> layer  10  24.18677   32.33419   39.20459   44.72136

quantile(r1)
#>     0%    25%    50%    75%   100% 
#>  100.0  585.0 1045.5 1537.0 2000.0

quantile(sqrt(r1))
#>       0%      25%      50%      75%     100% 
#> 10.00000 24.18677 32.33419 39.20459 44.72136

问题原因 #

raster::quantile处理栅格堆栈时默认采用按像元对齐剔除NA的逻辑：只要某一像元位置在任意一层栅格中为NA，所有层的该位置都会被统一剔除后再计算分位数。而单独计算单张栅格分位数时，仅剔除当前层的NA值，因此当不同层NA位置不一致时，两种计算方式的样本集不同，最终得到的分位数结果存在偏差。

解决方案 #

方案1：遍历堆栈每层单独计算分位数 #

如果需要得到和单张栅格计算一致的结果，可以遍历每层单独计算后再合并结果，示例代码如下：

# 遍历栅格堆栈每层计算分位数，合并为矩阵输出
stack_quantile <- t(sapply(1:nlayers(s), function(i) quantile(s[[i]])))
rownames(stack_quantile) <- names(s)
print(stack_quantile)

方案2：预处理统一NA位置 #

如果业务逻辑要求必须对堆栈整体计算，可以先对栅格堆栈做NA值统一处理：所有层的NA位置统一为所有层NA的并集，再调用quantile函数计算即可。

内容的提问来源于stack exchange，提问作者Jason