Cloud Masking with Sentinel 2

Overview

Objectives

Create a cloud mask for a Sentinel-2 scene, and merge with another scene. Perform a histogram match before merging both scenes. The location for this practice, is Cape Town, South Africa.

Data

For this practice, use the following files:

• Raster: S2A_MSIL2A_20200714T081611_N0214_R121_T34HBH_20200714T142158, S2B_MSIL2A_20200907T081609_N0214_R121_T34HBH_20200907T124303

Procedure

1. Install and/or load necessary packages

As for every new session of R, you should set your working directory and load the necessary libraries. Use ipak to install required packages.

ipak <- function(pkg){
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg)) 
    install.packages(new.pkg, dependencies = TRUE)
  sapply(pkg, require, character.only = TRUE)
}


packages <- c("rgdal", "raster","dplyr", "magrittr", "RStoolbox")
ipak(packages)

##     rgdal    raster     dplyr  magrittr RStoolbox 
##      TRUE      TRUE      TRUE      TRUE      TRUE

2. Load and stack rasters

Set working directory.

setwd("YOURDRIVE:\\yourpath\\yourfolder")

Import image bands. For this practice, we are only working with bands: Blue, Green and Red. Depending of the objectives of you work you want to do afterwards, you should adapt the bands and products.

# Import bands for the July image
# the path given is relative to your working directory (or your script if wd was not set)
# b2 is the blue band, b3 the green band and b4 the red band

july_b2 <- raster("data\\14july\\S2A_MSIL2A_20200714T081611_N0214_R121_T34HBH_20200714T142158\\S2A_MSIL2A_20200714T081611_N0214_R121_T34HBH_20200714T142158.SAFE\\GRANULE\\L2A_T34HBH_A026430_20200714T084249\\IMG_DATA\\R20m\\T34HBH_20200714T081611_B02_20m.JP2")

july_b3 <- raster("data\\14july\\S2A_MSIL2A_20200714T081611_N0214_R121_T34HBH_20200714T142158\\S2A_MSIL2A_20200714T081611_N0214_R121_T34HBH_20200714T142158.SAFE\\GRANULE\\L2A_T34HBH_A026430_20200714T084249\\IMG_DATA\\R20m\\T34HBH_20200714T081611_B03_20m.JP2")

july_b4 <- raster("data\\14july\\S2A_MSIL2A_20200714T081611_N0214_R121_T34HBH_20200714T142158\\S2A_MSIL2A_20200714T081611_N0214_R121_T34HBH_20200714T142158.SAFE\\GRANULE\\L2A_T34HBH_A026430_20200714T084249\\IMG_DATA\\R20m\\T34HBH_20200714T081611_B04_20m.JP2")

# Import bands for the September image

sept_b2 <- raster("data\\07sept\\S2B_MSIL2A_20200907T081609_N0214_R121_T34HBH_20200907T124303\\S2B_MSIL2A_20200907T081609_N0214_R121_T34HBH_20200907T124303.SAFE\\GRANULE\\L2A_T34HBH_A018308_20200907T084431\\IMG_DATA\\R20m\\T34HBH_20200907T081609_B02_20m.JP2")

sept_b3 <- raster("data\\07sept\\S2B_MSIL2A_20200907T081609_N0214_R121_T34HBH_20200907T124303\\S2B_MSIL2A_20200907T081609_N0214_R121_T34HBH_20200907T124303.SAFE\\GRANULE\\L2A_T34HBH_A018308_20200907T084431\\IMG_DATA\\R20m\\T34HBH_20200907T081609_B03_20m.JP2")

sept_b4 <- raster("data\\07sept\\S2B_MSIL2A_20200907T081609_N0214_R121_T34HBH_20200907T124303\\S2B_MSIL2A_20200907T081609_N0214_R121_T34HBH_20200907T124303.SAFE\\GRANULE\\L2A_T34HBH_A018308_20200907T084431\\IMG_DATA\\R20m\\T34HBH_20200907T081609_B04_20m.JP2")

Since the July bands are the ones that have a high cloud cover, the cloud mask will be performed on this scene. Therefore, load the “MSK_CLDPRB_20m.JP2” band of the QI DATA from the month of July.

# Import the mask with the cloud probability values (from the July image)
july_mask <- raster("data\\14july\\S2A_MSIL2A_20200714T081611_N0214_R121_T34HBH_20200714T142158\\S2A_MSIL2A_20200714T081611_N0214_R121_T34HBH_20200714T142158.SAFE\\GRANULE\\L2A_T34HBH_A026430_20200714T084249\\QI_DATA\\MSK_CLDPRB_20m.JP2")

Stack the bands and visualize them

# Stack
july_st <- stack(july_b2, july_b3, july_b4)
sept_st <- stack(sept_b2, sept_b3, sept_b4)

# Plot
par(mfrow=c(1,2))
# Use x11() to display plot in separate window
plotRGB(july_st, r=3,g=2,b=1, stretch="lin", main = "July scene")
plotRGB(sept_st, r=3,g=2,b=1, stretch="lin", main = "September scene")

Figure 1. True Color composition of July and September scenes. Cape Town, South Africa.

3. Determine appropriate cloud probability values and set them to NA

First open the image of the mask in QGIS, open the properties of the image, go to symbology, for “Render type” select “Paletted/Unique values” and click on “Classify”. Click on “OK” to close the window.

By using “Identify Features” (Ctrl + Shift + I) you can Identify the cloud probability value of the location you click on. Use this to determine the values for clouds. (Import the TCI image for comparison).

In this case, we have decided that all values that equal to 0 and 100 should be masked out. However, you might choose otherwise. Check the pixel values.

# set determined values of the mask to NA
# example:

july_mask_combi <- july_mask
july_mask_combi[july_mask == 0 |july_mask == 100] <- NA

plot(july_mask_combi) 

Figure 2. Mask values.

4. Mask out the clouds

# perfom the masking of NA values
july_masked <- raster::mask(july_st, july_mask_combi)

# visualize
plotRGB(july_masked, r=3,g=2,b=1, stretch="lin")

Figure 3. Masked values of July scene.

Export raster

Export your images and load them in QGIS. Compare.

writeRaster(july_masked, "results\\july_masked.tif", overwrite=TRUE)
writeRaster(sept_st, "results\\sept_st.tif", overwrite=TRUE)

5. Merge masked image, to fill in missing values

Once the image is free of clouds, it can be used for further analysis. Depending on the analysis to perform, the missing values can be compensated with values from another scene. In this case, we will use a cloud free scene from September of the same year to merge it with the masked July scene. Before performing a merge, it is important to apply a histogram matching on the image. Otherwise, the new acquired values will not be in the same range of the receiving image.

Histogram matching

# Use histMatch from RStoolbox. See: https://www.rdocumentation.org/packages/RStoolbox/versions/0.2.6/topics/histMatch
july_m_hist <- histMatch(july_masked, sept_st, xmask = NULL)

Histogram comparison

opar <- par(mfrow = c(1, 3), no.readonly = TRUE)
sept_st[,1:50] <- NA
redLayers <- stack(sept_st, july_masked, july_m_hist)[[c(1,4,7)]]
names(redLayers) <- c("September scene", "July scene masked", "July scene masked & matched")

hist(redLayers) 

Figure 4. Histogram comparison of the July masked scene and the September Rasterstack.

Merge raster

st_merged <- raster::merge(sept_st, july_m_hist)

Visualize the merged image

plotRGB(st_merged, r=3,g=2,b=1, stretch="lin")

Figure 5. July scene after cloud masking and merge.