在栅格堆栈的各层之间同步 NA

Question

我正在尝试开发一个函数来“同步”光栅堆栈各层之间的 NA，即确保对于堆栈的任何给定像素，如果一个图层具有 NA，则所有图层都应设置为 NA对于那个像素。

这在组合来自不同来源的栅格以进行物种分布建模时特别有用，因为某些模型无法正确处理 NA。

我找到了两种方法来做到这一点，但我发现它们都不令人满意。其中之一需要使用函数getValues，因此不适用于非常大的堆栈或内存不足的计算机。另一个更安全，但速度要慢得多。因此，我在这里问是否有人有改进我的尝试的想法？

这里有两种可能性：

1. 使用 getValues()

   syncNA1 <- function (x) 
   {
     val <- getValues(x)
     NA.pos <- unique(which(is.na(val), arr.ind = T)[, 1])
     val[NA.pos, ] <- NA
     x <- setValues(x, val)
     return(x)
   }
   

2. 使用 calc()

   syncNA2 <- function(y)
   {
     calc(y, na.rm = T, fun = function(x, na.rm = na.rm)
     {
       if(any(is.na(x)))
       {
         rep(NA, length(x))
       } else
       {
         x
       }
     })
   }
   

现在演示它们各自对同一堆栈的计算时间：

> system.time(
+ b1 <- syncNA1(a1)
+ )
   user  system elapsed 
   3.04    0.15    3.20 
> system.time(
+ b2 <- syncNA2(a1)
+ )
   user  system elapsed 
   5.89    0.19    6.08 

非常感谢您的帮助，

鲍里斯

Answer 1

对于名为“s”的堆栈，我将首先使用calc(s, fun = sum) 计算一个遮罩层，该遮罩层记录至少一个堆栈层中具有 NA 值的所有单元格的位置。然后mask() 将允许您将该遮罩应用到堆栈中的每一层。

这是一个例子：

library(raster)

## Construct reproducible data! (Here a raster stack with NA values in each layer) 
m <- raster(ncol=10, nrow=10)
n <- raster(ncol=10, nrow=10)
m[] <- runif(ncell(m))
n[] <- runif(ncell(n)) * 10
m[m < 0.5] <- NA
n[n < 5] <- NA
s <- stack(m,n)

## Synchronize the NA values
s2 <- mask(s, calc(s,fun = sum))

## Check that it worked
plot(s2)

Answer 2

我不知道速度，但您可以尝试转换为数组，加载 NA，然后再转换回来。伪代码：

xarray<-as.array(xstack)
ind.na<-which(is.na(xarray),array.ind=TRUE)
for(j in nrow(ind.na) ) {
    xarray[ind.na[j,1],ind.na[j,2],]<-NA
   }
nastack<-raster(xarray)

我没有验证那里的索引选择是否正确，也没有验证我是否正确转换回raster stack，但我希望你能明白。

编辑：我运行了一个时间测试，栅格为 1000x1000，但其他方面为 Josh 创建的。

 microbenchmark(josh(s),syncNA1(s),syncNA2(s),times=5)
Unit: milliseconds
       expr       min        lq    median        uq        max
    josh(s)  774.2363  789.1653  800.2511  806.5364   809.9087
 syncNA1(s)  652.3928  659.8327  692.3578  695.8057   743.9123
 syncNA2(s) 7951.3918 8291.7917 8604.2226 8606.3432 10254.4739
 neval
     5
     5
     5

Answer 3

我最终在 syncNA1 和 Josh 的解决方案之间构建了一个混合函数。 如果计算机没有足够的 RAM，此函数是内存安全的，但如果计算机有足够的 RAM，则可以更快地处理：

synchroniseNA <- function(x)
{
  if(canProcessInMemory(x, n = 2))
  {
    val <- getValues(x)
    NA.pos <- unique(which(is.na(val), arr.ind = T)[, 1])
    val[NA.pos, ] <- NA
    x <- setValues(x, val)
    return(x)
  } else
  {
    x <- mask(x, calc(x, fun = sum))
    return(x)
  }
}

但是，我凭经验确定data.frame 使用的内存量是光栅文件大小的两倍（对于canProcessInMemory() 的n 参数），但我不确定我是否正确.