1

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 '19 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 '19 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 '19 at 0:31
  • @Spacedman ah! It makes sense now. Thank you. – SA Khan Dec 13 '19 at 0:46
0

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.

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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