Ga naar inhoud

SAR


Sentinel-1 composite Mato Groso, Brazil. (Source: Source: BELSPO)

Dit onderdeel vormt een introductie tot het verwerken en visualiseren van SAR-data in Google Earth Engine. SAR staat voor Synthetic Aperture Radar en is een van de meest gebruikte radar-systemen binnen remote sensing. SAR biedt hoge-resolutie radarbeelden en wordt toegepast in een breed scala van onderzoek en monitoring.

Voordelen van SAR ten opzichte van multispectrale beelden

  • Wolkpenetratie: SAR-golven penetreren door wolken, waardoor ze geschikt zijn voor landbedekkingsmonitoring in alle weersomstandigheden. Dit is vooral nuttig in tropische gebieden tijdens regenperiodes.

  • Kronendakpenetratie: Afhankelijk van de golflengte kan SAR gedeeltelijk door het bladerdak dringen, wat betere schattingen van biomassa en vegetatiestructuur mogelijk maakt. Langere golflengtes penetreren doorgaans dieper.

  • Actieve sensortechnologie: SAR is een actieve vorm van remote sensing. Het systeem genereert zelf energie, waardoor observaties tijdsonafhankelijk zijn en zowel overdag als 's nachts uitgevoerd kunnen worden.

SAR wordt vaak toegepast in bosmonitoring, overstromingskartering en rampenbeheer. Voor meer achtergrondinformatie en praktijkvoorbeelden wordt The SAR Handbook (NASA, 2019) aanbevolen.

Onderstaande voorbeelden demonstreren het gebruik van SAR-data in Google Earth Engine met behulp van Sentinel-1.

Sentinel-1: Achtergrondinformatie

Over Sentinel-1


Sentinel-1 is een ESA Copernicus-missie bestaande uit twee satellieten die 180° tegenover elkaar in een baan om de aarde draaien. De waarnemingen worden uitgevoerd in de C-band (ca. 5,405 GHz), met een herhalingsfrequentie van minstens elke 12 dagen wereldwijd. De missie richt zich op:

  • Monitoring van zee-ijs en ijsbergen
  • Ondersteuning bij humanitaire hulp in crisistijd
  • Marien milieubeheer
  • Monitoring van landbouw en bosbouw

Voorbeeld 1 - Kartering van overstromingen met Sentinel-1

Achtergrond

SAR-data biedt een betrouwbare en snelle manier om de impact van overstromingen te analyseren. Doordat het onafhankelijk is van weersomstandigheden, kan cruciale informatie onmiddellijk worden afgeleid. In dit voorbeeld analyseren we overstromingen in Sindh, Pakistan, een van de zwaarst getroffen gebieden tijdens de overstromingen van 2022.

  • Studiegebied: Provincie Sindh, Pakistan .
    var ROI = ee.Geometry.Polygon( [[[66.14354335941121, 29.212308658016006], [66.14354335941121, 23.425656978642976], [71.50487148441121, 23.425656978642976], [71.50487148441121, 29.212308658016006]]], null, false);

1. Filter de Sentinel-1-collectie

Werk met de VV-polarisatie en selecteer de 'descending' orbiterichting.

var S1_VV = ee.ImageCollection('COPERNICUS/S1_GRD')
    .filter(ee.Filter.eq('instrumentMode', 'IW'))
    .filterBounds(ROI)
    .select('VV')
    .filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'))
            .map(function(image) {
              var edge = image.lt(-30.0);
              var maskedImage = image.mask().and(edge.not());
              return image.updateMask(maskedImage);
            });

    // Richting van S1-overtoch is van belang: hier enkel de 'descending' richting nemen
    var desc = S1_VV.filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'));

2. Selecteer beelden vóór en na de overstromingen

Gebruik beelden van augustus-september 2021 als referentie en vergelijk deze met beelden van augustus-september 2022. Filter S1_VV dus voor beide jaren. Hoeveel beelden bevatten de collectie voor beide jaren?

Code 
//Voor floods (2021)
var VV_2021 = desc.filterDate('2021-08-20','2021-09-15')
// Na floods (2022)
var VV_2022 = desc.filterDate('2022-08-20','2022-09-15')

// Bekijken van de geselecteerde collecties:
print('Tiles selected: Before Flood', VV_2021)//5 beelden
print('Tiles selected: After Flood', VV_2022)//12 beelden
  • Hierna maak je 2 "VV"-gepolariseerde beelden aan via een mean()-reducer: Voor en Na. Tevens zorgen we voor een 'speckle'-filter dat het peper-zout effect typerend aan radar wat reduceert. Hierna kunnen we het beeld visualiseren.
Code 
//------------------------------  DISPLAY PRODUCTS  ----------------------------------//
// Before and after flood SAR mosaic
Map.centerObject(ROI,8);
//Voor floods (2021)
var imgVV_2021 = VV_2021.mean().clip(ROI)
// Na floods (2022)
var imgVV_2022 = VV_2022.mean().clip(ROI)
// Peper-zout effect filteren (Speckle Filter)
var smoothing_radius = 50;
imgVV_2021 = imgVV_2021.focal_mean(smoothing_radius, 'circle', 'meters');
imgVV_2022 = imgVV_2022.focal_mean(smoothing_radius, 'circle', 'meters');
Map.addLayer(imgVV_2021,{min:-23,max:0}'VV August 2021')
Map.addLayer(imgVV_2022,{min:-23,max:0},'VV August 2022')

3. Visualiseer en bereken het verschil

Genereer een ratio van de beelden vóór en na de overstroming. Gebruik een drempelwaarde (bijvoorbeeld 1.25) om overstroomde gebieden te detecteren. Op basis van deze ratio kunnen we het overstroomde gebied in kaart brengen. Aangezien water doorgaans een lagere 'backscatter' heeft, zal dit resulteren in een nieuwe laag waar de overstroomde gebieden een hogere waarde hebben (voor (2021) = geen water = hogere backscatter, na (2022) = overstroomd = lagere backscatter).

  • Bepaal een treshold-waarde (via trial/error), (neem minstens 1.25, via de .gt()-functie). Het resultaat is een binair raster.
Code 
//------------------------------- Berekening overstromingsgebied -------------------------------//
// Calculate the difference between the before and after images
var difference = imgVV_2022.divide(imgVV_2021);
var threshold = 1.25; // Vooropgesteld na trial/error
var difference_binary = difference.gt(threshold);

4. Optimaliseer de overstromingslocaties

Gebruik aanvullende datasets (zoals de Global Surface Water-laag van JRC) om permanente waterlichamen en seizoensgebonden overstromingen te maskeren. Elimineer ruis door alleen gebieden met een minimale pixelconnectiviteit te behouden. Hiervoor kun je gebruik maken van connectedPixelCount().

          // Voeg de JRC-laag toe. Deze bevat informatie over seizoenale overstromingen, zodat pixels van gekende overstromingsgebieden gemaskeerd worden (deze worden immers niet meegerekend tot het rampgebied). Dit doen we dooe het weerhouden van de pixels die minstens 10 maanden per jaar overstroomd zijn.
          var swater = ee.Image('JRC/GSW1_0/GlobalSurfaceWater').select('seasonality');
          var swater_mask = swater.gte(10).updateMask(swater.gte(10));

          // Overstromingslaag waarbij vaste waterlichamen (water gedurende > 10 maaanden/jaar) een waarde van 0 krijgen
          var overstroming_mask = difference_binary.where(swater_mask , 0);
          // Definitief overstroomd gebied zonder pixels in eeuwige waterlichamen
          var overstroming = overstroming_mask.updateMask(overstroming_mask);

          // Bereken connectiviteit van pixels om diegene te elimineren die verbonden zijn met 8 of minder "buurpixels"
          // Deze operatie vermindert ruis van het overstromingsgebied 
          var connecties = overstroming.connectedPixelCount();    
          var flooded = overstroming.updateMask(connecties.gte(8));

5. Bereken het overstroomde oppervlak

Op basis van het verkregen flooded-beeld kunnen we de oppervlakte van het totale overstroomde gebied bepalen (in hectare) (Zie practicum 4).

Code 
// Calculate flood extent area
// Create a raster layer containing the area information of each pixel 
var flood_pixelarea = flooded.select('VV')
  .multiply(ee.Image.pixelArea());

// Som van de overstroomde pixels
// Gebruik de 'bestEffort: true' om de berekening te vereenvoudigen, voor een  
// accurater resultaat, zet bestEffort naar "false" en zet een hogere waarde voor 'maxPixels'. 
var flood_stats = flood_pixelarea.reduceRegion({
  reducer: ee.Reducer.sum(),              
  geometry: ROI,
  scale: 10, // resolutie Sentinel-1 pixel
  //maxPixels: 1e9,
  bestEffort: true
  });

// Converteer de oppervlakte van het overstromingsgebied naar hectare. rond ag.
var flood_area_ha = flood_stats
  .getNumber('VV')
  .divide(10000)
  .round();

6. Visualiseer het resultaat

// Difference layer
Map.addLayer(difference,{min:0,max:2},"Difference Layer",0);

// Flooded areas
Map.addLayer(flooded,{palette:"0000FF"},'Flooded areas');

print('Oppervlakte overstromingsgebied',flood_area_ha)
Bronvermelding

Inspiratie voor bovenstaande code via UN-Spider.

Link naar het volledige script:

https://code.earthengine.google.com/7768fbbfb0181a55ac19bd70e2ceec55

Opdracht

Oplossing

Via dit scriptje: https://code.earthengine.google.com/70bf32b42c7f233120352fc5fc1134f5

Voorbeeld 2 - Near-Real Time Forest Monitoring (NRTM)

Over NRTM

Near-Real Time Forest Monitoring of NRTM monitort bosbedekking in gebieden met veel illegale houtkap. Hierbij wordt gebruik gemaakt van de meest recente satellietdata om veranderingen in het kronendak snel te detecteren. Dit kan bijdragen aan het tijdig onderscheppen en stoppen van illegale houtkapactiviteiten.

Visualisatie van ontbossing

  • Open een nieuw scriptje in Earth Engine.

  • Zoek naar Sentinel-1 in de zoekbalk en lees kort de achtergrondinformatie door.


  • Sentinel-1 beelden bevatten verschillende polarisaties: HH, HV, VV, VH. Uit de theorieles weten we dat bijvoorbeeld 'VV' staat voor verticaal gepolariseerd signaal uit en verticaal gepolariseerd resultaat ontvangen. Selecteer de VV- en VH-polarisaties, en filter op een studiegebied (Region of Interest = ROI). In dit voorbeeld richten we ons op een tropisch regenwoudgebied ten oosten van het Brokopondo-meer in Suriname. Gebruik onderstaande code om de collectie in te laden en te filteren.


Code 
// Uit de Catalogus te kopiëren:
// Load Sentinel 1 C band SAR Ground Range collection (log scale, VV, descending)

// VV-polarisatie
var VV_coll = ee.ImageCollection('COPERNICUS/S1_GRD')
    .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
    .filter(ee.Filter.eq('instrumentMode', 'IW'))
    .filter(ee.Filter.eq('orbitProperties_pass','DESCENDING'))
    .select('VV')
    .map(function(image) {
      var edge = image.lt(-30.0);
      var maskedImage = image.mask().and(edge.not());
      return image.updateMask(maskedImage);
    });

// VH-polarisatie (extra)
var VH_coll = ee.ImageCollection('COPERNICUS/S1_GRD')
    .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))
    .filter(ee.Filter.eq('instrumentMode', 'IW'))
    .filter(ee.Filter.eq('orbitProperties_pass','DESCENDING'))
    .select('VH')
    .map(function(image) {
      var edge = image.lt(-30.0);
      var maskedImage = image.mask().and(edge.not());
      return image.updateMask(maskedImage);
    });

Analyse van bosbedekking over meerdere jaren

  • In dit voorbeeld bekijke we voor 3 tijdstippen de verandering in bosbedekking van dit gebiedje; 2019, 2021, 2023. Gebruik de maand augustus als referentieperiode en bereken de gemiddelde VV-polarisatie (in augustus) per jaar.
Code 
// Filteren op basis van locatie en tijdstip    
var VV_2019= VV_coll.filterBounds(ROI).filterDate('2019-08-01','2019-08-31').mean()
var VV_2021= VV_coll.filterBounds(ROI).filterDate('2021-08-01','2021-08-31').mean()
var VV_2023= VV_coll.filterBounds(ROI).filterDate('2023-08-01','2023-08-31').mean()
  • Visualiseer afzonderlijke beelden: Bekijk de verschillen tussen de beelden van de geselecteerde jaren.
Code 
var VV_Param = {"opacity":1,"bands":["VV"],"min":-12.695602621422937,"max":-2.5938492158251467,"gamma":1}; 
// Afzonderlijke beelden mappen: visueel weinig verschil te zien
Map.addLayer(VV_2019,VV_Param,'VV_2019',0)
Map.addLayer(VV_2021,VV_Param,'VV_2021',0)
Map.addLayer(VV_2023,VV_Param,'VV_2023',0)
  • Maak een composiet van de drie jaren: Combineer de beelden van 2019, 2021 en 2023 tot een RGB-composiet. Wat valt je op? Tip: om de afzonderlijke beelden (VV_2019,VV_2021 en VV_2023) samen te voegen tot één image kun je de .addBands() functie. Analyseer je resultaat.
Code 
//Afzonderlijke jaren samenvoegen

var VV_yearly = ee.Image(VV_2019).addBands(VV_2021).addBands(VV_2023);

Map.addLayer(VV_yearly, {min: -20, max: -5},'VV_2019_2021_2023')

Detectie van bosverlies

  • Bereken het verschil tussen de jaren 2023 en 2019: Pixels waar bos in 2019 aanwezig was maar in 2023 verdwenen is, worden negatief. Uit trial-and-error kun je bij grove benadering stellen de verschilwaarden van (VV_2023 -VV_2021) tussen -6 en -15 overeenkomen met ontbossing. Voer dus volgende stappen uit:
    1. Neem het verschil tussen 2023 en 2019.
    2. Rond de waarden af naar boven met de .ceil()-functie.
    3. Reclassificeer de waarden tussen -6 en -15 naar een waarde 1 met de .remap() functie. (Zie ook tresholding)
Code 
//Tussen 2 jaren: verschil nemen (2023 - 2019)
var VV_20192023 = VV_2023.subtract(VV_2019)

//Verschilbeeld illustreren
Map.addLayer(VV_20192023,{min:-13, max:10},'Diff_2021_2013')

// REMAP: om verlies te bepalen (-6 tot -15 komt grofweg overeen met ontbossing)

var VV_20192023= VV_20192023.ceil(); // afronden naar boven
var loss_20192023= VV_20192023.remap([-6,-7,-8,-9,-10,-11,-12,-13,-14,-15],[1,1,1,1,1,1,1,1,1,1])
Map.addLayer(loss_20192023,{palette: 'red'},'Loss_2019_2023')
  • EXTRA: Visualiseer nu ook 2 Sentinel-2-beelden (2019 en 2023) en controleer je gevonden ontboste gebieden.

Opdracht 2: Mangrove monitoring met SAR

  • Bekijk aan de hand van bovenstaand voorbeeld ook de wijzigingen in het mangrovegebied langs de volledige kustlijn in Suriname, tussen 2020-2023. Wat valt je op? Welk type ontbossing is dit?

  • Benader ook de aanwas van mangrove tussen deze jaren. De verschilwaarden tussen 7 en 15 kun je als grove grenzen nemen als wat 'aanwas' van mangrove is.

Oplossing

https://code.earthengine.google.com/49137e3139ed24346b3b5234333274b1