Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It only takes a minute to sign up.
Sign up to join this community
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)?
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
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.
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);
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
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
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().
