Ga naar inhoud

Cloud masking

1. Cloud Masking van Landsat data

Wolkbedekking is een grote barrière tijdens het analyseren en verwerken van (spectrale) satellietbeelden. Recente satellietdata komen vaak met automatische classificaties van de wolkbedekking, waardoor deze relatief eenvoudig uit het beeld verwijderd kan worden.

Earth engine bevat naast deze standaard ‘cloud masks’ ook algoritmes om de wolken en wolkschaduw te verwijderen uit het beeld. Een keten van filter, maskeer en reduce strategiën kan de aanwezigheid van wolken minimaliseren. Volgende 3 stappen kunnen onderdeel zijn van deze keten:

  1. Filteren op maximaal wolkenpercentage
  2. Cloudmasking algoritme toepassen om wolken te verwijderen
  3. Reduceren met een mediane reducer van de resterende collectie

1.1 Filteren van de ImageCollection op wolkbedekking

Een eerste optie is om een beeldcollectie te filteren (zie voorgaand) op wolkbedekking, waardoor enkel de beelden binnen een paalde range van wolkenpercentages worden weerhouden:

// Filteren van de Landsat 9 collectie tot beelden met maximaal 40% wolkbedekking
L9 = L9.filter(ee.Filter.lt('CLOUD_COVER', 40))

ee.filter()

Bij de voorgaande filters gebruikten we de redelijk eenvoudige functies .filterBounds() en .filterDate(), twee standaardfilters om op respectievelijk locatie (van een geometrie) en datum te filter.

De.filter() methode wordt gebruikt om te filteren op eender welke metadata-eigenschap die een beeldcollectie bevat. Binnen de filter-functie wordt vervolgens een Earth Engine filter-functie toegepast, die aangeeft op welke manier er gefilterd moet worden. Deee.Filter.lt('CLOUD_COVER', 40) = 'less than', weerhoudt dus alle beelden met minder dan 40% wolkbedekking.

Gebruik de Docs om het gebruik van deze functie verder te bekijken.

Beschikbare metadata kan steeds in de Earth Engine Catalog worden geraadpleegd. Voor Landsat 9 is dit: https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC09_C02_T1_L2#image-properties.

1.2 Cloud Mask functies

Als 2e stap kunnen de overgebleven wolken/wolkschaduwen per beeld worden ‘geknipt’ (cloudmask) door deze pixels naar een waarde 0 om te zetten. Voor Landsat 8 en 9 wordt dit gedaan op basis van het CFMask algoritme, dat automatishe detecties van wolken, wolkschaduwen en sneeuw of ijs nastreeft. Dit resultaat zit vervat in de metadata van elk beeld als binair raster, wat gebruikt kan worden om de wolken weg te knippen.

Toepassen van de cloudmask voor Landsat 9:

// TO DO: LAAD DE L9 collectie in, filter op locatie en gewensta data (bv Gent, zomer 2024)

var L9 = ... // Vul hier de start van je script aan.

// 1. Voeg een extra filter in obv wolkbedekking ('Cloud_cover')
L9 = L9.filter(ee.Filter.lt('CLOUD_COVER', 40)) 

// Definieren van CloudMask-functie (gegeven)

function maskL89sr(image) {
  // Gebaseerd op de QA-waarde, wat de uitkomst is van het FMASK algoritme
  // Develop masks for unwanted pixels (fill, cloud, cloud shadow).
  // De waarden voor de eerste 4 QA_pixelbanden moet gelijk zijn aan 0 (dus geen wolken/schaduw aanwezig)
  var qaMask = image.select('QA_PIXEL').bitwiseAnd(parseInt('11111', 2)).eq(0);
  var saturationMask = image.select('QA_RADSAT').eq(0);

  // Apply masks.
  return image.updateMask(qaMask).updateMask(saturationMask);
}

// Pas de functie over elk beeld binnen de collectie toe met de .map-functie:
var L9_masked = L9.map(maskL89sr);

// Eerste beeld uit collectie zonder Cloudmask:
Map.addLayer(L9.first(), trueColor, 'L9 - 1e beeld - No Cloudmask')

// Eerste beeld uit collectie mét Cloudmask:
Map.addLayer(L9_masked.first(), trueColor, 'L9 - 1e beeld - met Cloudmask')
Resulterend is een ImageCollectie met dezelfde beelden, maar waaruit de wolken gemaskeerd zijn (mask toegepast). Echter kunnen wel sommige wolkenranden nog zichtbaar zijn, die de mask-functies hebben gemist.

Opdracht 3.3

Maak de L9_masked collectie aan, en neem hiervan een .median() reducer. Visualiseer dit beeld. Merk je een verbetering in vergelijking met de voorgaande .median()-gereduceerde beelden, zonder de cloudmask?

De .map()-functie

In bovenstaand voorbeeld werd de cloudmask-functie toegepast door gebruik te maken van .map(). De .map() wordt steeds gebruikt om een functie (die op afzonderlijke beelden dient toegepast te worden, zoals maskL9sr) toe te passen over elk beeld binnen een ImageCollection afzonderlijk. Het is als het ware een veel efficiënte manier dan de aangemaakte functie te itereren via een for-loop.

2. Cloud Mask van Sentinel-2 data

Om wolken en wolkenschaduwen in Sentinel-2 beelden in Google Earth Engine te verwijderen, kun je het beste de Sentinel-2 Cloud Probablity dataset gebruiken. Deze dataset bevat een cloud probability score waarmee je pixels met wolken effectief kunt identificeren en maskeren.

Voor deze methode maak je eerst 2 functies aan:

  • getS2_SR_CLOUD_PROBABILITY: dat zowel de "originele" Sentinel-2 (harmonized) Surface Reflectance collectie oproept als de 'S2 cloud probability' collectie. Deze functie heeft geen geen parameters nodig en geeft een ImageCollection als resultaat waarbij de de 'cloud probability laag' toegevoegd wordt aan de Sentinel-2 beelden.

  • maskClouds: dat op basis van de toegevoegde 'cloudprobability'-laag een cloudmask aanmaakt en toepast.

Om dus tot een collectie te komen waarbij de 'cloudmask' is toegepast, gebruik je dus onderstaande code:

var getS2_SR_CLOUD_PROBABILITY = function () {
    var innerJoined = ee.Join.inner().apply({
         primary: ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED"),
         secondary: ee.ImageCollection("COPERNICUS/S2_CLOUD_PROBABILITY"),
         condition: ee.Filter.equals({
         leftField: 'system:index',
         rightField: 'system:index'
                  })
                });
     var mergeImageBands = function (joinResult) {
         return ee.Image(joinResult.get('primary'))
                  .addBands(joinResult.get('secondary'));
         };
     var newCollection = innerJoined.map(mergeImageBands);
         return ee.ImageCollection(newCollection);
         };


// Mask out clouds en scaling factor Sentinel-2
   var maskClouds = function(image) {
      var cloudProbabilityThreshold = 40;
      var cloudMask = image.select('probability').lt(cloudProbabilityThreshold);
      return image.updateMask(cloudMask)
          .multiply(0.0001) // Herschaling naar 0 - 1 range (Surface Reflectance)
                  .copyProperties(image, ["system:time_start"]);  // Tijdstip en beeldeigenschappen worden behouden
        };

var S2_coll = getS2_SR_CLOUD_PROBABILITY()
                .map(maskClouds)

Vervolgens kun je met S2_coll verder werken (filters toepassen, reducing, ..).

maskS2clouds

In het voorbeeldscriptje in de Sentinel-2 data catalogus van Earth Engine staat ook een voorbeeldscriptje gebruik makend van de functie maskS2clouds:

/**
 * 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).copyProperties(image, ["system:time_start"]);
}

Deze functie werkt niet voor alle Sentinel-2 beelden na Januari 2022. Gezien na deze datum de preprocessing algoritmen wat gewijzigd zijn, zijn ook bitmasks aangepast, waardoor bovenstaande code niet meer zorgt voor het gewenste resultaat.

Om dit op te lossen, maak je best gebruik van de cloudmaskprocedure met de S2_Cloud_Probability collectie zoals hierboven beschreven. Deze methode, die gebruik maakt van externe wolkbedekkingsdetectie-algoritmen (mooi woord voor scrabble), is nauwkeurigere resolutie.