I'm trying to perform a time-series analysis of Sentinel-2 satellite imagery over Water, and to do that I've been muddling through using functions to perform analysis on images over an Image Collection.
I previously got further with my analyses with individual images using .mosaic()
, but can't use that for the time-series. Using this earlier question and the GEE Guides I've been able to run NDVI and NDWI over a test area, create individual masks based on these, and 'Add' (Combine/Overlay/Union?) them into a new layer that has all the areas covered by both masks.
The issue I've run into is that when I have the final combined Mask Layer, the unmasked imagery doesn't match either the combined Mask, or the individual Masks. I suspect the imagery with the mask applied is an entirely different image from the original also, as much of the clipped image includes clouds in new areas which the original L1C image doesn't have.
For instance:
The Original Image The Combined NDVI/NDWI Mask (i.e. Areas of Water matching '1', red on the Mask) The same Combined Mask applied to the imagery (L1C Applied, in White-Gray), compared to the L1C Mask above
So, what is causing a different Mask area to appear instead? What section of my code is causing this error, and is there a fix or easier alternative to the way that I have to apply a mask to each individual image in the Image Collection?
I have included the code below, as well as a GEE Link. I've tried my best to organise and explain what each section does, but as I'm very new to GEE and .js, there's probably something evidently wrong here and there. I'm happy to provide more clarification on anything.
// ### IMPORTED GEOMETRY ###
var ShannonGeom =
/* color: #d6cfcf */
/* shown: false */
/* displayProperties: [
{
"type": "rectangle"
}
] */
ee.Geometry.Polygon(
[[[-9.061457413652322, 52.9400316020557],
[-9.061457413652322, 52.6131383784731],
[-8.26623608491697, 52.6131383784731],
[-8.26623608491697, 52.9400316020557]]], null, false);
// ### FILTER PARAMETERS ###
//These apply the parameters to stop repetitive data entry and searching for the correct line to enter them.
// CLOUD % format '00' , numbers only without quotes
var MIN_CLOUD_PERCENT = 50;
// START DATE & END DATE format: 'YYYY-MM-DD', include the quotes
var START_DATE = '2022-01-01';
var END_DATE = '2022-04-01';
// BOUNDS can either be a drawn geometry shape within GEE, or an imported SHP file
// See: https://developers.google.com/earth-engine/guides/table_upload#upload-a-shapefile
var BOUNDS = ShannonGeom;
// ZOOM is based on GEE's 1-24 level system. Larger number = Larger Zoom
var ZOOM = 10;
// PARAMETERS END - You don't need to change anything below this line to make the script function
// ### MISC ###
//Prints out the Parameters set
print('Filtering available Sentinel 2 imagery between ' + START_DATE + ' & ' + END_DATE + ' with ' + MIN_CLOUD_PERCENT + '% cloud over the area of interest.');
// Centre based on the Geometry (Region of Interest, Zoom Level)
Map.centerObject(BOUNDS, ZOOM, print('Map centered on Region of Interest'));
// Sets the default Map to Terrain mode (Roadmap overlain with hillsahde)
Map.setOptions("TERRAIN");
// ### EEPALETTES & LAYER VISUALISATION ###
/* Required for most layer visualisations
* See https://github.com/gee-community/ee-palettes for more information
*
* IF IT FAILS TO LOAD the PALETTES, LOAD THIS URL FIRST, THEN REFRESH THE PAGE:
* (https://code.earthengine.google.com/?accept_repo=users/gena/packages)
*/
var palettes = require('users/gena/packages:palettes');
// NDWI palette
var NDWIPalette = palettes.cmocean.Ice[7].reverse();
// NDVI palette
var NDVIPalette = palettes.colorbrewer.RdYlGn[10];
// Truecolour (R-G-B) Visualisation
var rgbVis = {
min: 0,
max: 0.35,
bands: ['B4', 'B3', 'B2'],
};
// NDWI Visualisation
var NDWIVis = {
min: -1,
max: 1,
bands: ['NDWI'],
palette: NDWIPalette
};
// NDVI Visualisation
var NDVIVis = {
min: -1,
max: 1,
bands: ['NDVI'],
palette: NDVIPalette
};
// NDVI Mask Visualisation
var NDVIMaskVis = {
min: 0, // Land Areas
max: 1, // Other Areas
bands: ['NDVI_Mask'],
palette: ['cccccc','088300'],
opacity: 0.65
};
// NDWI Mask Visualisation
var NDWIMaskVis = {
min: 0, // Land and Non-Water Areas
max: 1, // Water Areas
bands: ['NDWI_Mask'],
palette: ['cccccc','0000ff'],
opacity: 0.65
};
// L1C Water/Veg Mask Visualisation
var L1CMaskVis = {
min: 0, // Land and Non-Water Areas
max: 1, // Water Areas
bands: ['L1CMask'],
palette: ['cccccc','f90000'],
opacity: 0.65
};
// ### CLOUD MASKING ###
// Sentinel 2 Cloud Masking Function using the 60m Cloud Mask Band
/* Function to mask clouds using the Sentinel-2 QA band
* @param {ee.Image} image Sentinel-2 image
* @return {ee.Image} cloud masked Sentinel-2 image
*/
function maskS2clouds(image) {
var qa = image.select('QA60');
// Bits 10 and 11 are clouds and cirrus, respectively.
var cloudBitMask = 1 << 10;
var cirrusBitMask = 1 << 11;
// Both flags should be set to zero, indicating clear conditions.
var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(qa.bitwiseAnd(cirrusBitMask).eq(0));
return image.updateMask(mask).divide(10000);
}
print('Sentinel 2 Cloud Mask Function Complete');
// ### IMAGE COLLECTIONS ###
//Load and Map L1C imagery with the filter parameters applied
/*
* Load Sentinel-2 'Harmonized' Top Of Atomsphere (L1C) data
* Dataset details: https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S2_HARMONIZED
* HARMONIZED makes sure scenes after 25 January 2022 have the same DN ranges as older L1C scenes.
* Harmonised L1C Data is available from Sentinel 2 Launch (2015-06-23) onwards.
*/
var S2_L1C = ee.ImageCollection('COPERNICUS/S2_HARMONIZED')
// Filter by Date Period (YYYY-MM-DD)
.filterDate(START_DATE, END_DATE)
/*
* Pre-filter to get less cloudy granules
* 'Default' is aiming for 10% cloud
* Dependent on availability of cloud-free imagery in the time period set
* Longer periods will take longer to load
*/
.filter(ee.Filter.lte('CLOUDY_PIXEL_PERCENTAGE', MIN_CLOUD_PERCENT))
// Select only image tiles that fall within the Geometry (Region of Interest) to reduce processing time
.filterBounds(BOUNDS)
// Applies the S2 Cloud Masking Function to each image in the IC
.map(maskS2clouds)
// Clips each image in the IC by the Bounds to reduce processing time further
.map(function(image) {
return image.clip(BOUNDS);
});
print('Time, Date, Bounding, and Cloud Tile Filtering parameters set for Imagery');
// Add the pre-clipped IC's to the map
Map.addLayer(S2_L1C,rgbVis,'L1C');
// ### FUNCTIONS ###
// Add NDWI band to IC
var addNDWI = function(image) {
return image.addBands(image.normalizedDifference(['B3', 'B8']).rename('NDWI'));
};
// Add an NDVI band to IC
var addNDVI = function(image) {
return image.addBands(image.normalizedDifference(['B8', 'B4']).rename('NDVI'));
};
// ### L1C MASK FUNCTIONS ###
// Function to mask out NDWI (L1C)
var WaterMaskL1C = function(image) {
var NDWI = image.select(['NDWI']);
return image.addBands(NDWI.gte(0.1).rename('NDWI_Mask'));
};
// Function to mask out NDVI (L1C)
var VegMaskL1C = function(image) {
var NDVI = image.select(['NDVI']);
return image.addBands(NDVI.lte(0).rename('NDVI_Mask'));
};
// Function to combine both L1C Masks
var CombinedL1CMask = function(image) {
var NDWI_Mask = image.select(['NDWI_Mask']);
var NDVI_Mask = image.select(['NDVI_Mask']);
// Adds the 2 masks together, and then clamps the output to be binary to keep it suitable for masking
return image.addBands(NDWI_Mask.add(NDVI_Mask).clamp(0,1).rename('L1CMask'));
};
// Function to update the clipping of imagery to just the water areas, as defined by the combined Mask
//.rename(['B1_Msk','B2_Msk','B3_Msk','B4_Msk','B5_Msk','B6_Msk','B7_Msk','B8_Msk','B8A_Msk','B9_Msk','B10_Msk','B11_Msk','B12_Msk','QA10_Msk','QA20_Msk','QA60_Msk','NDWI_Msk','NDVI_Msk','NDVI_Mask2','NDWI_Mask2','L1CMask2'])
var L1CMasking = function(image) {
var Mask = image.select(['L1CMask']);
return image.updateMask(Mask.eq(1));
};
// Applies all the functions to respective image collections
var S2_L1C_Func = S2_L1C.map(addNDWI).map(addNDVI).map(VegMaskL1C).map(WaterMaskL1C).map(CombinedL1CMask);
// Add the individual new bands to the Map
// NDVI
Map.addLayer(S2_L1C_Func, NDVIVis, 'L1C NDVI');
// NDWI
Map.addLayer(S2_L1C_Func, NDWIVis, 'L1C NDWI');
// NDVI Mask
Map.addLayer(S2_L1C_Func, NDVIMaskVis, 'L1C NDVIMask');
// NDWI Mask
Map.addLayer(S2_L1C_Func, NDWIMaskVis, 'L1C NDWIMask');
// Combo Mask
Map.addLayer(S2_L1C_Func, L1CMaskVis, 'L1C Mask');
/* Creating a new var to run the mask apply function, because otherwise for some reason it retroactively
* breaks all the earlier separate ND VI/WI masks
*/
var S2_L1C_Masked = S2_L1C_Func.map(L1CMasking);
// Clipped hopefully?!
Map.addLayer(S2_L1C_Masked, rgbVis, 'L1C MaskApplied');
I've also attempted updating the masks to the imagery separately as shown below, but instead that just brings up another, different and incorrect, masked image (evident as clouds are included within the masked L1C image that aren't in the original L1C image, and even the tides have changed in the Estuary on the bottom right corner!).
/* Instead of combining the two Masks prior to applying them to the IC, this updates the IC's Mask by the two separate layers
* together. This is only applied to the "image" bands (B1-B12), not the other newly created bands.
*/
var MaskBandsL1C = function(image) {
return image.select(['B1','B2','B3','B4','B5','B6','B7','B8','B9','B10','B11','B12','NDVI_Mask','NDWI_Mask'])
.updateMask(image.select('NDVI_Mask'))
.updateMask(image.select('NDWI_Mask'));
};
// Combo Mask
var ICMaskL2A = S2_L2A_Func.map(MaskBandsL2A);
var ICMaskL1C = S2_L1C_Func.map(MaskBandsL1C);
Map.addLayer(ICMaskL2A, rgbVis, 'L2A Masked');
Map.addLayer(ICMaskL1C, rgbVis, 'L1C Masked');
The alternative (and incorrect) result of the second masking example
So the mask is 'working' in that it is providing the areas fitting my conditions, but it is bringing up an image from at least a different time of day (at worst, an image from a different day altogether) - could it also be possible to have the masked image(s) aligned by system:time_start
as I will have to use that to create a Time Series chart?
Thank you for your time in providing any explanation as to what is happening.
Answer Update
I posted this similar question earlier to the sister GIS Stack Exchange (I'm not sure if that is bad etiquette, I just thought asking here afterwards might give better reach), and Olga Danylo was able to explain the issue I was having. If anyone else is having a similar problem, I recommend you read their answer and give them a +1 if they helped you too.
In short:
The masking is working, it's just that when I .map
a layer without a specific filtered date or image, any result is pulled up and displayed, and it just so happens that my separate masks appeared to match up until the final combined Mask layer. Adding a system:time_start
property to the first cloud mapping function and then filtering all the Map.addLayer
's to the same date (i.e. Map.addLayer(S2_L1C_Func.filterDate('2022-03-30', '2022-04-01')
gave the map results I wanted. Olga gives a better explanation than I, so I really recommend just reading their answer on the link above.