2

Can someone please tell me how will I discard the clouds from Landsat 8 images using R? I have used the code

LC08_20151128_new <- cloudMask(LC08_20151128, threshold = 0.8, blue = "layer.1", tir = "layer.9", buffer = NULL, plot = FALSE)

I read that for landsat 8 layer 1 and layer 9 might be useful to detect the clouds. Is it correct? I got the image below after runninng the code

LC08_20151128_new

And then I masked the image with the cloudMask

LC08_20151128_crop <- mask(LC08_20151128,LC08_20151128_new, maskvalue = NA)

But it seems that the code retains the area with the cloud cover instead of removing it. Please see the RGB image below

 LC08_20151128_crop

2

Below code provides a reproducible example using the Landsat 8 scene 'LC08_L1TP_195025_20130707_20170503_01_T1' (Collection 1 Level-1) which is roughly centered on Frankfurt am Main, Germany, and was downloaded from USGS Earth Explorer. In order to discard all cloudy pixels (ie set them to NA), which I assume is your overall goal, specify either

  • maskvalue = 1 or
  • maskvalue = NA, inverse = TRUE

inside your raster::mask() call. Have a look at the function's help pages for further details.

library(satellite)
library(RStoolbox)

## sample data
lc08 = list.files("LC08_L1TP_195025_20130707_20170503_01_T1", 
                  pattern = ".TIF$", full.names = TRUE)
lc08 = sortFilesLandsat(lc08)
lc08 = stack(lc08[c(1:7, 9:11)]) # stack without band 8 (15 m)

## clip image
rext = extent(c(530824, 563455, 5518663, 5545203))
lc08 = crop(lc08, rext)

## mask clouds and visualize
cmsk = cloudMask(lc08, threshold = .4, blue = 2, tir = 9) # bands 2, 10
wocl = mask(lc08, cmsk$CMASK, maskvalue = 1)

library(mapview)
viewRGB(wocl, r = 4, g = 3, b = 2)

viewRGB

  • Many thanks for your reply. Do I need to crop the original image with woc1? Or should I continue my analysis with with it? I have used the code you suggested but it did not remove the cloud – D.Banerjee Aug 3 '17 at 13:41
  • In my example, wocl contains the original image with clouds removed (see resulting image created by mapview::viewRGB), which is likely what you are aiming at. In case the clouds haven't been captured in your particular case, make sure to have a look at plot(cmsk) and adjust threshold to the range of NDTCI accordingly (see ?cloudMask). – fdetsch Aug 3 '17 at 13:52
  • Many thanks for your help. This is really very useful. – D.Banerjee Aug 3 '17 at 22:14

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.