4

Using this code, I am able to get a cloud-free composite, but it still includes shadows in my image.

var collection = ee.ImageCollection("COPERNICUS/S2")
  .filterBounds(ee.Geometry.Rectangle(-85.14404296875, 46.9502622421856, 
  -71.87255859375, 41.60722821271717));
var S2_cloud = collection.min();

How do I create a Sentinel-2 cloud free, cloud-shadow free composite or scene on Earth Engine (GEE)?

2
  • 1
    Can you include screenshots of areas in question?
    – MaryBeth
    Feb 21 '18 at 20:57
  • You can find an example for cloud masking Sentinel-2 in the Code Editor Examples. Scripts tab > Examples > Cloud Masking Feb 26 '18 at 2:54
7

As far as I know, Sentinel 2 has no computed cloud-shadow mask, so I ommit that part. But for the rest, I think you should at least mask out clouds, and filter a little.

I've made public a repo I have to perform this: users/fitoprincipe/geetools

First the core:

var computeQAbits = function(image, start, end, newName) {
    var pattern = 0;

    for (var i=start; i<=end; i++) {
        pattern += Math.pow(2, i);
    }

    return image.select([0], [newName]).bitwiseAnd(pattern).rightShift(start);
};

var sentinel2 = function(image) {

  var cloud_mask = image.select("QA60");
  var opaque = computeQAbits(cloud_mask, 10, 10, "opaque");
  var cirrus = computeQAbits(cloud_mask, 11, 11, "cirrus");
  var mask = opaque.or(cirrus);

  return image.updateMask(mask.not());
}

You can use the function sentinel2 in a map function.

to your question:

Using min(), which may not be the best approach to generate a composite, (but that's another story so I am not diggin on it), this would be a complete 'solution':

var s2mask = require('users/fitoprincipe/geetools:cloud_masks').sentinel2;

var AOI = ee.Geometry.Rectangle(-85.14404296875, 46.9502622421856, 
  -71.87255859375, 41.60722821271717);

var collection = ee.ImageCollection("COPERNICUS/S2")
  .filterBounds(AOI)
  .filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', 70)
  .map(s2mask());

var image = ee.Image(collection.min());
Map.addLayer(image, {bands:['B8', 'B12', 'B4'], min:0, max:3000});
Map.centerObject(AOI);
2

I am quite new to Google Earth Engine but here is the solution I came up nonetheless. Instead of using min() function to composite image, I use percentile. This should reduce shadow captured using min() function while the overall image should not be too choppy (If in the time interval, there is a section where there is no clouds)

One thing I want to improve in the future is to somehow smooth out the mosaic effect on RGB image. But currently, I am out of idea. Another one is rgb_p might not represent the true color as each channel is selected based on its percentile separately. Better method is to group RGB and select it whole

The code are the following

// Setup
var s2mask = require('users/fitoprincipe/geetools:cloud_masks').sentinel2;
var cloud_threshold = 80; // Filter image from collection to have cloud less than 80 percent
var rgb_threshold = 0.23; // Cut-off, Probably unmasked cloud
var rgb_percentile = 35; // Select RGB at 35th percentile
var ndvi_percentile = 90; // Select NDVI at 90th percentile

// Define helpers
function get_collection(start_date, end_date) {  
    var collection = ee.ImageCollection("COPERNICUS/S2")
        .filterDate(start_date, end_date)
        .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', cloud_threshold))
        .map(s2mask())
        .map(function (image) {
            return image.divide(10000);
        });
    return collection.map(function(image) {
        var ndvi = image.normalizedDifference(['B8', 'B4']).rename('NDVI');
        image = image.addBands(ndvi.toFloat());
        return image.toFloat();
    })
}

function get_final(collection) {
    var final = collection.max(); // Baseline value
    // Select value
    collection = collection.map(function(img) {
        return img.updateMask(
            img.select("B2").lt(rgb_threshold)
            .bitwiseAnd(img.select("B3").lt(rgb_threshold))
            .bitwiseAnd(img.select("B4").lt(rgb_threshold)))
    });
    var rgb_p = collection.select(["B2", "B3", "B4"])
                              .reduce(ee.Reducer.percentile([rgb_percentile]))
    rgb_p = rgb_p.select([
                            "B2_p" + rgb_percentile,
                            "B3_p" + rgb_percentile,
                            "B4_p" + rgb_percentile,
                        ], ["B2", "B3", "B4"])
    var ndvi_p = collection.select(["NDVI"])
                                .reduce(ee.Reducer.percentile([ndvi_percentile]))
                                .select([
                                    "NDVI_p" + ndvi_percentile
                                ], ["NDVI"])
    // Replace old bands
    final = final.addBands({
        srcImg: rgb_p.select("B2").toFloat().rename("B2"),
        overwrite: true
    });
    final = final.addBands({
        srcImg: rgb_p.select("B3").toFloat().rename("B3"),
        overwrite: true
    });
    final = final.addBands({
        srcImg: rgb_p.select("B4").toFloat().rename("B4"),
        overwrite: true
    });
    final = final.addBands({
        srcImg: ndvi_p.select("NDVI").toFloat().rename("NDVI"),
        overwrite: true
    });
    return final.toFloat();
}

// Operation
var collection = get_collection("2018-01-01", "2018-03-01");
var final = get_final(collection).toFloat();

// Visualize
Map.addLayer(final.select(["B4", "B3", "B2"]), {
    min: 0,
    max: 0.3
}, "RGB");
Map.addLayer(final.select(["NDVI"]), {
    min: -1,
    max: 1,
    palette: ['blue', 'white', 'green']
}, "NDVI");

For shadow-free image, I am not sure. But I think would go something like this. For every pixel threshold outlier first and then select bands using percentile

As of right now in Google Earth Engine, in my opinion it is much better to use "S2_SR" image collections as its corrected for atmospheric correction

If I am missing something, please feel free to comment. Thank you

Update 1

I no longer use that function. Instead I use hollstein_s2 decision tree to classy shadow region and new S2_CLOUD_PROBABILITY which is processed by GEE itself. The update code are as follows

var dt = require('users/fitoprincipe/geetools:decision_tree');

function hollstein_S2_shadow(img) {
    // Ref: https://github.com/fitoprincipe/geetools-code-editor/blob/master/cloud_masks
    var difference = function (a, b) {
        var wrap = function (img) {
            return img.select(a).subtract(img.select(b))
        }
        return wrap
    }
    var ratio = function (a, b) {
        var wrap = function (img) {
            return img.select(a).divide(img.select(b))
        }
        return wrap
    }
    //1
    var b3 = img.select('B3').lt(3190)
    //2
    var b8a = img.select('B8A').lt(1660)
    var r511 = ratio('B5', 'B11')(img).lt(4.33)
    //3
    var b3_3 = img.select('B3').lt(5250)
    var s37 = difference('B3', 'B7')(img).lt(270)
    //4
    var r15 = ratio('B1', 'B5')(img).lt(1.184)
    var s911 = difference('B9', 'B11')(img).lt(210)
    var s911_2 = difference('B9', 'B11')(img).lt(-970)

    var dtf = dt.binary({
        1: b3,
        21: b8a,
        22: r511,
        31: s37,
        34: b3_3,
        41: s911_2,
        42: s911,
        46: r15
    }, {
        'shadow-1': [
            [1, 1],
            [21, 1],
            [31, 1],
            [41, 0]
        ],
        'shadow-2': [
            [1, 1],
            [21, 1],
            [31, 0],
            [42, 0]
        ],
        'shadow-3': [
            [1, 0],
            [22, 0],
            [34, 1],
            [46, 0]
        ],
    }, 'hollstein')
    var results = dtf(img)
    return img.updateMask(results.select("shadow").not())
}

function get_S2_SR_clean(criterion, maximum_cloud_prob, mode, masked_shadow, rgb_percentile) {
    // Helper functions
    mode = mode || "mosaic";
    masked_shadow = masked_shadow || true;
    rgb_percentile = rgb_percentile || 20;

    function maskClouds(img) {
        var clouds = ee.Image(img.get('cloud_mask')).select('probability');
        var isNotCloud = clouds.lte(maximum_cloud_prob);
        return img.addBands(clouds).updateMask(isNotCloud);
    }

    function maskEdges(s2_img) {
        return s2_img.updateMask(
            s2_img.select('B8A').mask().updateMask(s2_img.select('B9').mask()));
    }

    function maskShadow(s2_img) {
        return hollstein_S2_shadow(s2_img)
    }

    function get_percentile_rename(img, bands, p) {
        // Assume p is list with size of 1
        return img.select(bands).reduce(ee.Reducer.percentile(p)).rename(bands);
    }
    var allBands = [
        "B1", "B2", "B3", "B4", "B5", "B6", "B7", "B8",
        "B8A", "B9", "B11", "B12", "AOT", "WVP", "SCL",
        "TCI_R", "TCI_G", "TCI_B", "MSK_CLDPRB", "MSK_SNWPRB",
        "QA10", "QA20", "QA60", "probability"
    ];

    function getBandsAnythingBut(excludeBands) {
        var bandNames = [];
        for (var i = 0; i < allBands.length; i++) {
            if (excludeBands.indexOf(allBands[i]) == -1) {
                bandNames.push(allBands[i])
            }
        }
        return bandNames;
    }
    // Operations
    // -> Import
    var S2_SR = ee.ImageCollection("COPERNICUS/S2_SR");
    var S2_Clouds = ee.ImageCollection("COPERNICUS/S2_CLOUD_PROBABILITY");
    // -> Filtering operations
    S2_SR = S2_SR.filter(criterion).map(maskEdges);
    S2_Clouds = S2_Clouds.filter(criterion);
    // -> Join images
    var S2_SR_Cloud_Mask = ee.Join.saveFirst('cloud_mask').apply({
        primary: S2_SR,
        secondary: S2_Clouds,
        condition: ee.Filter.equals({
            leftField: 'system:index',
            rightField: 'system:index'
        })
    });

    // var s2CloudMasked = ee.ImageCollection(S2_SR_Cloud_Mask).map(maskClouds).mosaic();
    var s2CloudMasked = ee.ImageCollection(S2_SR_Cloud_Mask).map(maskClouds);
    if (masked_shadow) {
        s2CloudMasked = s2CloudMasked.map(maskShadow);
    }
    if (mode == "mosaic") {
        s2CloudMasked = s2CloudMasked.mosaic();
    } else if (mode == "median") {
        s2CloudMasked = s2CloudMasked.median();
    } else if (mode == "rgb_percentile") {
        var tmp = get_percentile_rename(s2CloudMasked, ["B4", "B3", "B2"], [rgb_percentile]);
        s2CloudMasked = s2CloudMasked.select(getBandsAnythingBut(["B4", "B3", "B2"])).addBands(tmp)
    }
    return s2CloudMasked;
}

I have use this composite method in Smoothing/interpolating across images in an ImageCollection to remove missing data, feel free to look at the cleaned up code for example usage

1

http://www.gisandbeers.com/generar-imagenes-satelite-sin-nubes/

Landsat and sentinel Cloudness mosaic

Cloudness mask is similar as obtained with Principe`s Script, but using median instead min().

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.