从 rasterstack 中删除具有许多 NA 的栅格图层

Question

我有一个 RasterStack 有几千层。但是，有些层有很多 NA，所以我会通过设置阈值来排除这些层。我使用 for 循环正确地做到了，但这非常慢。我试图用 calc 来做，但我的功能失败了。这是我的试验，如果有任何提示可以加快处理速度，我将不胜感激。

library(raster)
lst<-stack(r1,r2,r3,r4) # 
lst_new<-stack()
for (i in 1: nlayers(lst)){
   # total number of cells without NA
   no_NA<-length(lst[[i]][!is.na(lst[[i]])]) #
   if(no_NA >= 14652){ # 97% 
     l<-lst[[i]]
     lst_new<-stack(lst_new,l)
     }
   }
  #This code works OK but slow for big rasterstack. So I tried the
  # following using calc function    
  remove.badL<-function(x){
          no_NA<-length(x[is.na(x)])
          if(no_NA >= 14652){
          return(x)
          }
         }
  lst_new<-calc(lst,fun=remove.badL) 
# this is the error I got

错误（函数（类，fdef，mtable）：无法为签名“"RasterBrick"、"NULL"”

的函数“writeValues”找到继承方法

如果有任何建议，我将不胜感激。谢谢

Answer 1

这是一种方法，使用 is.na、cellStats 和条件 RasterStack 子集。

首先，让我们创建一些示例数据：

library(raster)
s <- stack(replicate(10, raster(matrix(runif(1e6), 1e3))))
s[s > 0.95] <- NA # insert some NAs

我们可以 return 每层 NA 个细胞的数量：

cellStats(is.na(s), sum)

有了这些知识，我们可以在子集操作中使用这些计数：

thr <- 14652
s2 <- s[[which(cellStats(is.na(s), sum) < thr)]]

层数少于 thr（此处为 14652）NA 个单元格将保留在新堆栈中，s2，而那些多于 NA 的单元格将保留在新堆栈中被抛弃。

将所有这些应用到您的数据，您应该能够使用：

lst_new <- lst[[cellStats(is.na(lst), sum) < 14652]]

Answer 2

在处理庞大的数据集时，

cellStats 可能并不总是最佳选择。例如，将@jbaums 样本数据扩展到 n = 100 层，在我的机器上需要相当长的时间。

## sample data, n = 100
library(raster)

set.seed(10)
s <- stack(replicate(100, raster(matrix(runif(1e6), 1e3))))
s[s > 0.95] <- NA

## set na limit (e.g., 5% of all cells)
limit <- 0.05 * ncell(s)


### cellStats -----

system.time(
  id1 <- cellStats(is.na(s), sum) < limit
)
# user  system elapsed 
# 28.794   0.253  29.050

您可以例如使用并行 foreach 手动创建指示少量缺失数据的索引向量，而不是使用 cellStats。

### parallel version -----

## open parallel backend
library(doParallel)
cl <- makeCluster(detectCores() - 1)
registerDoParallel(cl)

## loop over layers in parallel
system.time(
  id2 <- foreach(i = unstack(s), .packages = "raster", 
                 .combine = "c") %dopar% {
                   sum(is.na(i[])) < limit
                 }
)
# user  system elapsed 
# 0.337   0.005   3.802

如您所见，后一种方法的执行速度要快得多，同时 returns 结果相同。

## similarity check
identical(as.logical(id1), id2)
[1] TRUE

接下来唯一要做的就是关闭并行后端

## deregister parallel backend
stopCluster(cl)

并根据导出的索引向量创建 s 的子集。

## data subset
s[[which(id2)]]

从 rasterstack 中删除具有许多 NA 的栅格图层

Remove raster layers with many NAs from rasterstack

r

raster