从具有多层的栅格对象中删除异常值

Question

我有一个空间分辨率非常高的 NetCDF 文件：

dimensions : 2160, 4320, 9331200, 172  (nrow, ncol, ncell, nlayers)
resolution : 0.0833333, 0.0833333  (x, y)
extent     : -180, 179.9998, -90, 89.99993  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 

最终目标是计算所有 172 层的平均值（即时间维度）。在此之前，我想逐个单元格地从时间序列中删除异常值。例如使用公式：

x ± 2.5 * std(x)

，其中 x 是每个单元格中数据的时间序列。

我尝试将 Raster 对象转换为数据框，但考虑到分辨率，它非常慢。我还尝试使用 raster 包中的函数

clamp

和

reclassify

，但是用于子集的值对于整个 Raster 对象是固定的：

x <- clamp(x, lower = 0, upper = 2.5)

我正在寻找的是根据每个单元格的时间序列（即上面的动态下限和上限范围）删除异常值的功能。

有没有一种方法可以直接使用 Raster 格式，使用 Raster 或 Terra 包？

有一个类似的问题here但在那种情况下Raster对象没有多层（对应于时间维度）。

Answer 1

示例数据（在问 R 问题时请始终包含一些数据）

library(terra)
s <- rast(ncol=10, nrow=10, nlyr=30)
set.seed(1)
values(s) <- rnorm(size(s), 10)
s[1] <- NA

获取均值和标准差，并使用 ifel

mn <- mean(s)
sdv <- stdev(s) * 1.5
lower <- mn-sdv
upper <- mn+sdv

x <- ifel(s < lower, NA, ifel(s > upper, NA, s))

检查单元格的结果

i <- 57
round(s[i], 2)
#  lyr.1 lyr.2 lyr.3 lyr.4 lyr.5 lyr.6 lyr.7 lyr.8 lyr.9 lyr.10 lyr.11 lyr.12 lyr.13
#1  9.63    11   7.6  8.83 10.51  7.67    10 12.02  9.23  10.96  11.21  10.12  10.61
#  lyr.14 lyr.15 lyr.16 lyr.17 lyr.18 lyr.19 lyr.20 lyr.21 lyr.22 lyr.23 lyr.24 lyr.25
#1   8.83  10.26  11.73   9.49   8.73   9.02  11.75   7.92   9.95    9.6    9.5   9.46
#  lyr.26 lyr.27 lyr.28 lyr.29 lyr.30
#1   9.72  11.38  10.61   7.35  10.36

cbind(lower[i], upper[57])
#      mean    mean
#1 7.987305 11.6816

round(x[i], 2)
#  lyr.1 lyr.2 lyr.3 lyr.4 lyr.5 lyr.6 lyr.7 lyr.8 lyr.9 lyr.10 lyr.11 lyr.12 lyr.13
#1  9.63    11    NA  8.83 10.51    NA    10    NA  9.23  10.96  11.21  10.12  10.61
#  lyr.14 lyr.15 lyr.16 lyr.17 lyr.18 lyr.19 lyr.20 lyr.21 lyr.22 lyr.23 lyr.24 lyr.25
#1   8.83  10.26     NA   9.49   8.73   9.02     NA     NA   9.95    9.6    9.5   9.46
#  lyr.26 lyr.27 lyr.28 lyr.29 lyr.30
#1   9.72  11.38  10.61     NA  10.36

或使用

app

f <- function(x) {
    m <- mean(x) + c(-1,1) * sd(x)
    x[x < m[1] | x > m[2]] <- NA
    x
}
y <- app(s, f)

应该也可以做到

 z <- clamp(s, lower, upper)

我会把它添加到“terra”