How can I get the image with less cloud coverage over a ROI in a time range with Sentinel 2?

Question

I am trying to compute a script to get the cleanest image of Sentinel 2 in a prescribed time range. Currently I tried with 2 approaches: a rough one and an automatic one, but I get two different results.

The rough one gives me the output that I want (without any clouds over the ROI), but the automatic one NO! Essentially it uses the QA60 band to calculate cloud coverage % over the ROI, but something is wrong (and I don't know why)

Any suggestion on how to fix the problem?

The script: https://code.earthengine.google.com/4ef624d82002ec759f633ad68992b6de

(I don't know if you can read my asset. In case: https://code.earthengine.google.com/?asset=projects/ee-matteogiomo04/assets/TEST_SET)

var Train = 
    /* color: #d63000 */
    /* shown: false */
    ee.Geometry.Polygon(
        [[[11.404322420267174, 46.32569478725022],
          [11.402605806497643, 46.2198529823167],
          [11.66559103598983, 46.28537560568817],
          [11.648081575540612, 46.34489524610421]]]),
    sentinel2 = ee.ImageCollection("COPERNICUS/S2_SR"),
    ROI = ee.FeatureCollection("projects/ee-matteogiomo04/assets/TEST_SET"),
    visParamsTrue = {"bands":["B4","B3","B2"],"min":0,"max":2500,"gamma":1.1};




// 1) Set some parameters

var start1 = ee.Date('2018-06-01')
var end1 = ee.Date('2018-07-30')

var ROI = ee.FeatureCollection(ROI).geometry()
print(ROI)


/////// FIRST ROUGH METHOD /////////

var image = sentinel2  
  .filterDate(start1, end1)
  .filterBounds(Train)
  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 10))
  .filter(ee.Filter.lt('DARK_FEATURES_PERCENTAGE', 20))
  .filter(ee.Filter.lt('CLOUD_SHADOW_PERCENTAGE', 20))
  .filter(ee.Filter.lt('MEAN_SOLAR_ZENITH_ANGLE', 30));

print('Less cloudy', image)
Map.addLayer(ROI, {}, 'Test Area');
Map.addLayer(image, visParamsTrue, 'Rough');
Map.setOptions('SATELLITE');
Map.centerObject(Train, 12);



///////// AUTOMATIC METHOD /////////

var images = sentinel2.filterBounds(ROI.centroid());    // filter in the centroid of the ROI

function calcCloudStats(img) {
   var imgPoly = ee.Algorithms.GeometryConstructors.Polygon( 
        ee.Geometry( img.get('system:footprint') ).coordinates() 
        )
   var intersection = ROI.intersection(imgPoly, ee.ErrorMargin(0.5))


    var areaImg = img.select('QA60').multiply(ee.Image.pixelArea())

    var stats = areaImg.reduceRegion({
       reducer: ee.Reducer.sum(),
       geometry: ROI,
       scale: 10,
       maxPixels: 1e9
    })

    var cloudPercent = ee.Number(stats.get('QA60')).divide(imgPoly.area())
    var coveragePercent = ee.Number(intersection.area()).divide(ROI.area())
    var cloudPercentROI = ee.Number(stats.get('QA60')).divide(ROI.area())

    img=img.set('1_CLOUDY_PERCENTAGE', cloudPercent)
    img=img.set('2_ROI_COVERAGE_PERCENT', coveragePercent)
    img=img.set('3_CLOUDY_PERCENTAGE_ROI', cloudPercentROI)

    return(img)
}

var cloudIndexes = images.map(calcCloudStats)

var image1 = cloudIndexes
  .filterDate(start1, end1)
  .sort('3_CLOUDY_PERCENTAGE_ROI')
  .first();

Map.addLayer(image1, visParamsTrue, 'Automatic');

Answer 1

You could also use the scene classification (SCL) band to do this:

// Define a function to mask out image clouds and shadows.
var cloudMask = function(img) {
  var qa = img.select("SCL");
  var cloud = qa.lt(8).and(qa.gt(3)).rename("percent_clear")
  var clear = cloud.reduceRegion(ee.Reducer.mean(), area, 20)
  return img.set(clear);
};

var s2masked = s2_col.map(cloudMask)
  .filterMetadata("percent_clear", "greater_than", 0.95)