Remove raster layers with many NAs from rasterstack

Question

I've a RasterStack with several thousand layers. However, some of the layers have a lot of NAs so I would to exclude those layers by setting a threshold. I did it correctly using a for loop, but this is very slow. I was trying to do it with calc, but my function failed. Here are my trials and I would appreciate any hint in order to make the processing faster.

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

Error in (function (classes, fdef, mtable) : unable to find an inherited method for function ‘writeValues’ for signature ‘"RasterBrick", "NULL"’

I would appreciate any suggestions. Thanks

fdetsch fdetsch · Accepted Answer · 2016-03-01T17:33:52

cellStats might not always be the best option when dealing with huge datasets. Extending @jbaums sample data to n = 100 layers, for instance, takes quite some time on my machine.

## 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

Instead of using cellStats, you could e.g. create the index vector that indicates low amounts of missing data manually using a parallelized foreach.

### 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

As you can see, the latter approach performs considerably faster and, at the same time, returns identical results.

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

The only thing left to do then is close the parallel backend

## deregister parallel backend
stopCluster(cl)

and create a subset of s based on the derived index vector.

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

Remove raster layers with many NAs from rasterstack

2 Answers