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?

  • You have over 5 billion cells in that stack. Perhaps you haven't waited long enough. How long does it take on a single layer? Multiply by 768. Quick test on my laptop means it would take about 6 minutes. How's your patience?
    – Spacedman
    Dec 12, 2019 at 20:40
  • @Spacedman yeah it has 6.8 million cells. I waited for more than 3 hours! I guess it shouldn't take that long.
    – SA Khan
    Dec 13, 2019 at 0:25
  • It has 6.8 million cells in each of the 768 layers making 5 billion in total. It might be quicker to do each layer separately since then not so much will be loaded into memory at once... Hmmm...
    – Spacedman
    Dec 13, 2019 at 0:31
  • @Spacedman ah! It makes sense now. Thank you.
    – SA Khan
    Dec 13, 2019 at 0:46

1 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"
[1] "raster"

So when I created it s1 was a RasterStack:

> class(s1)
[1] "RasterStack"
[1] "raster"

I do the replacement and...

> s1[s1 %in% c(22,33,44)] = NA
> class(s1)
[1] "RasterBrick"
[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.

  • Thank you for the illustration of the problem. I will create a for loop and apply over the raster stack.
    – SA Khan
    Dec 13, 2019 at 0:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.