Several clouds yet in COPERNICUS/S2_SR despite restrict filters

[Leprechault]

I'd like to know if there is yet something to do not to download COPERNICUS/S2_SR images with cloud or shadows. In my reproducible example:

# Packages
library(tidyverse)
library(rgee)
library(sf)
library(raster)
ee_Initialize(drive=TRUE)

# Function for remove cloud and shadows ------------------------------------------
getQABits <- function(image, qa) {
  # Convert decimal (character) to decimal (little endian)
  qa <- sum(2^(which(rev(unlist(strsplit(as.character(qa), "")) == 1))-1))
  # Return a single band image of the extracted QA bits, giving the qa value.
  image$bitwiseAnd(qa)$lt(1)
}
s2_clean <- function(img) {
  # Select only band of interest, for instance, B2,B3,B4,B8
  img_band_selected <- img$select("B[2-4|8]")

  # quality band
  ndvi_qa <- img$select("QA60")

  # Select pixels to mask
  quality_mask <- getQABits(ndvi_qa, "110000000000")

  # Mask pixels with value zero.
  img_band_selected$updateMask(quality_mask)
}

# Define a Region of interest
roi <-ee$Geometry$Point(-52.19032,-30.25413)$buffer(500)

# Sentinel-2 MSI dataset into the Earth Engine’s public data archive ------------              
s2 <- ee$ImageCollection("COPERNICUS/S2_SR")

# Select S2 images ---------------------------------------------------------------
s2_roi  <- s2$
  filterBounds(roi)$
  filter(ee$Filter$lte("CLOUDY_PIXEL_PERCENTAGE", 1))$
  filter(ee$Filter$date(as.character(as.Date("2019-12-04")), as.character(as.Date("2020-05-03"))))$
  map(s2_clean)

s2_roi_add_area <- s2_roi$map(
  function(img) {
    img$set("area", img$clip(roi)$geometry()$area())
  }
)

# Get the dates and IDs of the selected images ------------------------------------
area <- floor(ee_utils_py_to_r(roi$area(maxError=1)$getInfo()))
ic_date_gt_area <- s2_roi_add_area$filterMetadata("area", "greater_than", area)

nimages <- ic_date_gt_area$size()$getInfo()
nimages 

ic_date <- ee_get_date_ic(s2_roi_add_area)
ic_date

# Download the results
s2_ic_local <- ee_imagecollection_to_local(
  ic = s2_roi_add_area,
  scale = 10,
  region = roi,
  add_metadata = TRUE)
# 

# Plot the image collected
# Cloud 1
plotRGB(stack("20191204T133221_20191204T133222_T22JCM.tif"),  r = 3, g = 2, b = 1, stretch = "lin")

#Cloud 2
plotRGB(stack("20200502T133231_20200502T133228_T22JCM.tif"),  r = 3, g = 2, b = 1, stretch = "lin")

It's clear the clouds in 2019-04-12 and 2020-05-02 dates, despite the filters filter(ee$Filter$lte("CLOUDY_PIXEL_PERCENTAGE", 1) and floor(ee_utils_py_to_r(roi$area(maxError=1)$getInfo())) applied. Please, any help with it?

[csaybar]

Maybe using the SLC band and COPERNICUS_S2_CLOUD_PROBABILITY can give you better results. You would have to:

1) Convert SLC and COPERNICUS_S2_CLOUD_PROBABILITY to a binary mask (0: clear, 1:mask). 2) Merge both products using AND operator. 3) Create your own "CLOUDY_PIXEL_PERCENTAGE" using ee$Image$reduceRegion. 4) Filtering your ic according to this new metadata. 5) if the results are unsatisfactory, try applying dilation to the merged mask product (step 2).

[Leprechault]

Thanks very much, @csaybar, but look a little bit complicated. Please, do you have any examples in R? I don't find any examples in R syntax from the web.

[csaybar]

Hi @Leprechault use rgeeExtra to have a more R-friendly interface to rgee.

library(rgeeExtra)
library(rgee)

# 1. Define your ROI
ROI <- ee$Geometry$Rectangle(
  coords = c(440030, 4568820, 445120, 4573910), 
  proj = "EPSG:32719", 
  geodesic = FALSE
)

# 2. load S2 collection
ic <- ee$ImageCollection("COPERNICUS/S2_SR") %>% 
  ee$ImageCollection$filterBounds(ROI) %>% 
  ee$ImageCollection$filterDate("2021-09-01", "2022-01-01")

add_new_property <- function(img) {

  # 3. Use sen2cor cloud mask
  sen2cor_cloudmask <- img[["SCL"]]
  bmask01 <- sen2cor_cloudmask == 8 | sen2cor_cloudmask == 9 | sen2cor_cloudmask == 10

  # 4. Use sen2cloudless cloud mask
  # sen2cloudless recommends 0.4 (as I remenber), but it varies according the ROI.
  bmask02 <- img[["MSK_CLDPRB"]] > 0.5

  # 5. Merge both cloud masks
  bmask <- bmask01 | bmask02

  # 6. Estimate the new CLOUDY_PIXEL_PERCENTAGE
  newproperty <- ee$Image$reduceRegion(
    image = bmask,
    geometry = ROI, # PUT HERE YOUR ROI!
    reducer = ee$Reducer$mean(),
    bestEffort = TRUE
  ) %>% 
    # The rgeeExtra results have the band name "layer".
    # It works as well the raster R package!
    ee$Dictionary$get("layer") %>% 
    ee$Number()

  # 7. Set the property
  img$set("newCLOUDY_PIXEL_PERCENTAGE", newproperty)
}

new_ic <- ic$map(add_new_property)

ncloud_coverage <- unlist(new_ic$aggregate_array("newCLOUDY_PIXEL_PERCENTAGE")$getInfo())
nold_coverage <- unlist(new_ic$aggregate_array("CLOUDY_PIXEL_PERCENTAGE")$getInfo())

[Leprechault]

So nice @csaybar, thanks again!! Problem solved