栅格砖对于功能来说太大：每个网格单元中最大值的 20% 阈值

Question

我有一个栅格砖 b 全球 fAPAR 数据，它是 35 年的日常数据（即大）。

目标

有没有一种方法可以让我循环遍历这个积木并一次在 365 层（天）上运行一个函数？最终目标是运行一个函数来识别每个网格单元的值超过该网格单元每年最大值的 20% 阈值的层（日期）。

所以代码需要：

1. 计算每个网格单元在一年中的最大值
2. 返回每个网格单元超过该最大值 20% 的图层（日期）

代码

以下代码查找每年的最大值，但由于其大小（每日、全球、0.5 度数据），它不适用于我的数据。关键可能是一次处理 365 层，而不是整块砖（如果可能的话）。

# Example 10-year brick:
b <- brick(nrow=108,ncol=132,nl=3650)
values(b) <- runif(ncell(b)*nlayers(b))

i <- rep(1:nlayers(b), each=365)
# if incomplete last set of days:  
i <- i[1:nlayers(b)]
x <- stackApply(b, i, max)

我的目标与this 和this 非常相似。考虑到大小问题，我还想知道是否需要将我的数据划分为单独的年份并在多个文件上运行一些东西。

Answer 1

这可能是因为光栅评估内存需要错误。您可以使用rasterOptions 设置较低的内存限制。

您可以使用i 循环多年（我注意到您认为它们都有 365 天）

yrs <- unique(i)
out <- list()
for (j in yrs) {
  bb <- b[[which(i==j)]]
  out[[j]] <- max(bb) 
}
out <- stack(out)

您也可以尝试使用terra 执行此操作（替换raster 的新包）

library(terra)
x <- rast(b)
z <- tapp(x, i, max)

我认为这没有很好的定义：

返回每个网格单元超过最大值 20% 的图层（日期）

因为可以有很多这样的日期。下面，out1 获取所有内容，out1 获取第一个匹配项：

library(terra)
x <- rast(nrow=108,ncol=132, nl=3650)
set.seed(1)
values(x) <- runif(size(x))
i <- rep(1:nlayers(b), each=365)
i <- i[1:nlayers(b)]

zz <- tapp(x, i, function(i) 0.2 * max(i))

yrs <- unique(i)
out1 <- out2 <- list()
for (j in yrs) {
    xx <- x[[which(i==j)]]
    rr <- zz[[j]] < xx
    out1[[j]] <- rr
    out2[[j]] <- which.lyr(rr)
}

out1 <- rast(out1)
out2 <- rast(out2)