栅格：仅在其他 RasterLayer 中不是 NA 时才在 RasterStack 上计算

Question

我有一个 RasterStack s1 由 400 层组成，其中包含来自岛屿的数据。栅格的范围被裁剪为岛屿的范围，但由于其形状不规则，只有大约 20% 的像素实际上是陆地面积并具有数据值；其他 80% 是水和NA。

我还有一个水陆掩码 lwm (RasterLayer)，其中土地编码为 1，水编码为 NA。

我想在s1 上进行不同类型的基于单元格的计算，但注意到这些需要很长时间才能完成。为了加快速度，计算应该只对陆地区域进行，而水域应该总是NA。在伪代码中：

for each cell:
    if cell is land
        do calculation
    if cell is water
        return(NA)

要求是memory-safety。

这里有一些示例数据来说明问题：

library(raster)
# generate data
lwm <- raster(nrow = 5, ncol = 5)
lwm[] <- c(rep(NA, 10), rep(1, 5), rep(NA, 10))

r1 <- raster(nrow = 5, ncol = 5)
r1[] <- runif(ncell(r1)) * 10
r2 <- raster(nrow = 5, ncol = 5)
r2[] <- runif(ncell(r2)) * 10
s1 <- stack(r1, r2)
s1 <- mask(s1, lwm)

# this works, but all NA-values on water are also unnecessarily evaluated 
calc(s1, function(x) {sum(!is.na(x))})

Answer 1

经过一番尝试，我终于找到了一个非常适合我的解决方案，并且从整体处理时间中节省了相当多的时间：

library(raster)
# generate data
lwm <- raster(nrow = 50, ncol = 50)
lwm[] <- 1

# replace 80% with NA values
lwm[sample(1:ncell(lwm), round(0.8 * ncell(lwm)))] <- NA

r1 <- raster(lwm)
r1[] <- runif(ncell(r1)) 
r1_list <- replicate(400 , r1)
s1 <- stack(r1_list)
s1 <- mask(s1, lwm)

# this works, but all NA-values on water are also unnecessarily evaluated 
system.time(r_sum1 <- calc(s1, function(x) {sum(x)}))
#user  system elapsed 
#0.14    0.00    0.14 

## new approach:

# stack land-water-mask with RasterStack 
s1_lwm <- stack(lwm, s1)

# function to check if first element of vector is NA; if yes, return NA; if no, do calculation
fun1 <- function(y) {
  if (!is.na(y[1])) {
    y = y[-1]
    return(sum(y))
  } else {
      return(NA)
    }
}

system.time(
  r_sum2 <- calc(s1_lwm, fun = fun1)
)
# user  system elapsed 
# 0.4     0.0     0.4 

# results are identical 
identical(r_sum1[], r_sum2[])

Answer 2

这是一个棘手的问题，不幸的是我没有一个直接的解决方案给你。

您可以对岛屿进行多次裁剪（即 2-3）以最小化 NA 值，并在每个裁剪的栅格上分别进行计算并拼接结果。

或者另一种选择是进行并行计算，这将显着加快处理速度：

  #initialize cluster
  #number of cores to use for clusterR function (max recommended: ncores - 1)
  beginCluster(3)

  #calculation
  result <- clusterR(s1, calc, args=list(fun=function(x) {sum(!is.na(x))}))

  #end cluster
  endCluster()

由于您要求使用内存安全的解决方案，因此您应该查看当您只有一个核心时分配了多少 RAM，然后估计可以在多少核心上运行计算，这样您就不会耗尽内存.

祝你好运！希望对您有所帮助。