Replacing cloud values with NA in a raster stack

Question

I am removing clouds from a raster stack of 768 pixel_qa bands which are provided by USGS as part of surface reflectance product. The details of raster stack (qa) are:

> qa
class      : RasterStack 
dimensions : 2520, 2714, 6839280, 768  (nrow, ncol, ncell, nlayers)
resolution : 30, 30  (x, y)
extent     : 205755, 287175, 3950955, 4026555  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=43 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 

Before I use this qa raster stack to mask NDVI stack, I use the following code to replace cloud values to NA:

> qa[qa %in% c(224, 480, 992)] <- NA 

where 224 is cloud pixel value for Landsat 4-7 and 480 & 992 for Landsat 8.

The code runs forever! What could be wrong here?

Answer 1

Doing it a layer at a time in a for loop seems much faster. Here's an example.

Make a 500x500 raster with values from 1 to 200:

> set.seed(999); r = raster(matrix(sample(200, 500*500,TRUE),500,500))

Make that into a stack of five of them:

> s = stack(r,r,r,r,r)

Make two copies:

> s1 = s
> s2 = s

The first one, replace all the values with a single replacement:

> system.time({s1[s1 %in% c(22,33,44)] <- NA})
   user  system elapsed 
  1.716   0.000   1.730 

The second one, loop over nlayers(s2)

> system.time({for(i in 1:nlayers(s2)){s2[[i]][s2[[i]] %in% c(22,33,44)] <- NA}})
   user  system elapsed 
  0.117   0.000   0.118 

Ten times faster. The resulting values are identical:

> identical(s1[],s2[])
[1] TRUE

Although the raster objects aren't. Hmm...

> s1 = s
> class(s)
[1] "RasterStack"
attr(,"package")
[1] "raster"

So when I created it s1 was a RasterStack:

> class(s1)
[1] "RasterStack"
attr(,"package")
[1] "raster"

I do the replacement and...

> s1[s1 %in% c(22,33,44)] = NA
> class(s1)
[1] "RasterBrick"
attr(,"package")
[1] "raster"

Its now a RasterBrick. This might explain the slowness - R is creating another object the same size and possible now you have two copies in memory. Doing the layer at a time method doesn't change the object to a RasterBrick so might be more efficient in its memory handling.