1

我有一个全球 fAPAR 数据的栅格块b ，它是 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非常相似。考虑到大小问题，我还想知道是否需要将我的数据划分为单独的年份并在多个文件上运行某些东西。

4

1 回答 1

2

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

你可以循环多年（我注意到你认为他们都有 365 天）i

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)
于 2021-09-09T14:22:57.660 回答