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