Ga naar inhoud

SAR

Synthetic Aperture Radar


Sentinel-1 composite Mato Grosso, 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 radarsystemen binnen remote sensing. SAR biedt hoge-resolutiebeelden en wordt toegepast in een breed scala aan monitoringtoepassingen.

Voordelen van SAR ten opzichte van multispectrale beelden

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

  • (Beperkte) kronendakpenetratie: Afhankelijk van de golflengte kan SAR gedeeltelijk door het bladerdak dringen, wat betere schattingen van biomassa en vegetatiestructuur mogelijk maakt. Langere golflengtes (bijv. L-band) penetreren doorgaans dieper dan C-band (Sentinel-1).

  • Actieve sensortechnologie: SAR is een actieve vorm van remote sensing. De sensor zendt zelf energie uit en meet het teruggekaatste signaal (backscatter). Observaties zijn daardoor onafhankelijk van zonlicht en kunnen zowel overdag als ’s nachts gebeuren.

Typische toepassingen zijn bosmonitoring, overstromingskartering en rampenbeheer. Voor meer achtergrondinformatie en praktijkvoorbeelden wordt The SAR Handbook (NASA, 2019) aanbevolen.

Onderstaande voorbeelden tonen hoe je SAR-data van Sentinel-1 in Google Earth Engine kunt gebruiken.

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. Omdat radarbeelden onafhankelijk zijn van wolken en zonlicht, kunnen we zelfs tijdens extreme weersomstandigheden data verkrijgen.

In dit voorbeeld analyseren we overstromingen in Sindh, Pakistan, een van de zwaarst getroffen regio’s 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

We werken met de VV-polarisatie (gevoelig voor glad water én verticale structuren) en selecteren 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'))   // Enkel descending behouden
    .map(function(image) {
      // Randmasker: hele donkere randen (artefacten) uitsluiten
      var edge = image.lt(-30.0);
      var maskedImage = image.mask().and(edge.not());
      return image.updateMask(maskedImage);
            });

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

Bij overstroming worden eerder droge of vegetatie-bedekte gebieden plots waterig → de VV-backscatter daalt sterk.

Een intuïtieve maat is de ratio vóór/na de overstroming:
- Voor (2021): geen water → hoge backscatter
- Na (2022): overstroomd → lage backscatter
⇒ ratio VV_2021 / VV_2022 wordt groot in overstroomde zones

Bepaal dus een treshold-waarde (via trial & error), (neem minstens 1.25, via de .gt()-functie). Het resultaat is een binair raster.

Code 
//------------------------------- Berekening overstromingsgebied -------------------------------//
// Ratio vóór (droger) / na (natter)
// Let op: we werken op dB-beelden; voor een echt fysisch correcte ratio zou je naar lineair
// kunnen converteren, maar voor dit practicum volstaat deze benadering.

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().

          // JRC Global Surface Water: 'seasonality' = # maanden/jaar met water
          var swater = ee.Image('JRC/GSW1_0/GlobalSurfaceWater').select('seasonality');
          var swater_mask = swater.gte(10).updateMask(swater.gte(10));

          // Overstromingslaag waarbij permanente waterlichamen (water gedurende > 10 maanden/jaar) de waarde van '0' krijgen
          var overstroming_mask = difference_binary.where(swater_mask , 0);

          // Definitief overstroomd gebied (zonder permanente waterlichamen)
          var overstroming = overstroming_mask.updateMask(overstroming_mask);

          // Connectiviteit: kleine “spikkels” wegfilteren
          // Deze operatie vermindert ruis van het overstromingsgebied 
          var connecties = overstroming.connectedPixelCount();    
          var flooded = overstroming.updateMask(connecties.gte(8));

5. Bereken het overstroomde oppervlak

Nu berekenen we de oppervlakte van het overstroomde gebied (in hectare).

Code 
// Pixeloppervlak (m²) voor overstroomde pixels
var flood_pixelarea = flooded.select('VV')
  .multiply(ee.Image.pixelArea());

// Som van de overstroomde pixels
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)  // Omrekenen naar ha
  .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