[GIS] Cloudmask in Landsat 8 images

landsat 8r

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

Best Answer

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

Related Question