光栅:只有在其他RasterStack中没有NA的情况下才能在RasterLayer上计算

Question

我有一个由400个层组成的RasterStack s1，其中包含来自一个岛屿的数据。栅格的范围被裁剪到岛的范围，但由于其不规则的形状，只有大约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()

既然你要求一个内存安全的解决方案，你应该看看当你只分配一个单一的内核时，内存被分配了多少，然后估计你可以在多少个内核上运行计算，这样你就不会耗尽内存。

祝好运!希望能帮上忙。