Raster：仅在其他 RasterLayer 中不适用时才在 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

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

您可以对岛屿进行多次裁剪（即 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，然后估计你可以运行计算多少个核心，所以你赢了't 运行内存不足。

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

Answer 2

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

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[])

Raster：仅在其他 RasterLayer 中不适用时才在 RasterStack 上计算

Raster: Calculation on RasterStack only if not NA in other RasterLayer

if-statement

r

r-raster