Přeskočit obsah

Ukázka práce v Google Earth Engine

Cíl cvičení

  • Seznámit se cloudovým nástrojem Google Earth Engine
  • Umět detekovat spáleniště pomocí dat DPZ

Základní přehled

Google Earth Engine je cloudová platforma pro geoprostorovou analýzu, která uživatelům umožňuje vizualizovat a analyzovat satelitní snímky naší planety.

Code Editor

Code Editor je webové vývojové prostředí Google Earth Engine využívající programovací jazyk JavaScript.

Ukázky JavaScript kódu

print('Hello world!');

Zobrazení dat MODIS:

var dataset = ee.ImageCollection('MODIS/061/MOD09GA')
  // filter by dates
  .filterDate('2023-04-01', '2023-06-01');

// Select the bands you want to visualize
var visParams = {
  bands: ['sur_refl_b01', 'sur_refl_b04', 'sur_refl_b03'],
  min: 0,
  max: 5000,
  gamma: 1.4,
};

// add to the map window
Map.addLayer(dataset, visParams, 'MODIS');

Zobrazení dat Landsat 9 s oblačností menší než 10 % a s přeškálováním:

var dataset = ee.ImageCollection('LANDSAT/LC09/C02/T1_L2')
  // filter by dates
  .filterDate('2023-07-01', '2023-08-31')
  // filter by cloud cover
  .filter(ee.Filter.lessThan('CLOUD_COVER', 10));

// Applies scaling factors
function applyScaleFactors(image) {
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0);
  return image.addBands(opticalBands, null, true)
              .addBands(thermalBands, null, true);
}

dataset = dataset.map(applyScaleFactors);

// Select the bands you want to visualize
var visParams = {
  bands: ['SR_B4', 'SR_B3', 'SR_B2'],
  min: 0.0,
  max: 0.3,
};

// add to the map window
Map.addLayer(dataset, visParams, 'Landsat 9');

Zobrazení dat Sentinel-2 s oblačností menší než 10 % a jejich oříznutí na zájmové území:

var sentinel_col = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
  // filter by dates
  .filterDate('2023-05-01', '2023-06-30')
  // filter for the drawn rectangle
  .filterBounds(geometry)
  // filter by cloud cover
  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 10));

// Define a function to clip each image in the collection to the geometry
var clipToGeometry = function(image) {
  return image.clip(geometry);
};

// Use the map function to apply the clip function to each image in the collection
var clipped_col = sentinel_col.map(clipToGeometry);

// Select the bands you want to visualize
var visParams = {
  bands: ['B4', 'B3', 'B2'],
  min: 0.0,
  max: 3000,
};

// add to the map window
Map.addLayer(clipped_col, visParams, 'Sentinel-2');

Spoustu dalších příkladů lze nalézt přímo v prostředí Google Earth Engine v části Scripts → Examples.

Detekce spálenišť pomocí Google Earth Engine

V rámci tohoto cvičení využijeme Google Earth Engine pro detekování oblastí zasažených požárem. Konkrétní oblastí našeho zájmu bude národní park České Švýcarsko, kde došlo k požáru v období 23. 7. 2022 - 12. 8. 2022 (v německé části požár trval od 25. 7. 2022 - 19. 8. 2022). Práci začneme tím, že si vytvoříme nový soubor, do kterého budeme náš kód psát. Nový soubor vytvoříme v části Scripts → NEW → File. Protože jsme si ale nevytvořili repozitář, do kterého se soubor uloží, Google Earth Engine nás nejprve vyzve k vytvoření právě nového repozitáře. Teprve poté můžeme vytvořit nový soubor.

Vytvořený soubor najdeme ve Scripts → Owner, kde si ho i můžeme otevřít.

Výběr zájmové oblasti a získání dat

Zájmové území si označíme přidáním bodu do mapy. Bod přidáme pomocí funkce Přidat značku, která se nachází v levé horní části mapového okna. Bod umístíme přibližně mezi obce Hřensko a Mezná, tedy na místo, kde požár zhruba vypukl.

Pomocí funkce Edit layer properties si můžeme vytvořený bod přejmenovat z geometry např. na point.

Nyní si můžeme vyhledat potřebná data. Konkrétně budeme pracovat s daty Sentinel-2. Bohužel v roce 2022 nebylo v daném místě mnoho bezoblačných scén, proto nastavíme datum na obdomí mezi 1. 4. 2022 - 31. 10. 2022. Jako limit pro oblačnost nastavíme 30 %.

// Filter the Sentinel-2 image collection
var S2 = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED') // import data
                  .filterBounds(point) // spatial filter
                  .filterDate('2022-04-01', '2022-10-31') // temporal filter
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30)); // metadata filter

// Print data to the Console
print("Scenes:", S2);

// Add True color RGB composite to the map
Map.addLayer(S2.first(),
            {min:0,max:3000,bands:"B4,B3,B2"}, 
            "RGB");

Pokud bychom si chtěli zobrazit jiný než první prvek z ImageCollection, je potřeba si ImageCollection převést na List. To můžeme udělat přidáním následujícího kódu.

// Convert ImageCollection to List
var S2List = S2.toList(S2.size());
print(S2List);

// Get second image from the list
Map.addLayer(ee.Image(S2List.get(1)),
            {min:0,max:3000,bands:"B4,B3,B2"}, 
            "RGB 2");

Mezi vrstvami v mapovém okně lze jednoduše přepínat pomocí nástroje Layers v pravé části mapového okna.

Při zavření a opětovném otevření Google Earth Engine není mapa vycentrována na naše území. Přidáním následujícího kódu do skriptu se mapa vycentruje na zobrazenou scénu.

// Center the map to the image
Map.centerObject(S2.first());

Maskování oblačnosti

Mohli jsme si všimnout, že na některých scénách je oblačnosti poměrně dost. To pro jakékoliv analýzy není zrovna nejlepší, a je proto potřeba oblačnost zamaskovat. Existuje více možností, jak scény maskovat, a lze k tomu použít i různé vrstvy obsažené v produktech Sentinel-2. Zde je např. ukázka maskování oblačnosti pomocí vrstvy QA60, což je maska oblačnosti s rozlišením 60 m. Bohužel pro data od ledna 2022 není tato vrstva stále k dispozici. Lze ale využít i jiné vrstvy, jako např. MSK_CLDPRB, což je pravděpodobnost výskytu oblačnosti v pixelu s rozlišením 20 m. Přidáme si tedy do kódu funkci, která zamaskuje oblačnost pomocí vrstvy MSK_CLDPRB, zvolíme nějaký threshold určující, kolikaprocentní pravděpodobnost oblačnosti budeme brát v potaz, a aplikujeme masku oblačnosti na naše data. Můžeme poté porovnat zamaskovaná a nezamaskovaná data.

// Filter the Sentinel-2 image collection
var S2 = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED') // import data
                  .filterBounds(point) // spatial filter
                  .filterDate('2022-04-01', '2022-10-31') // temporal filter
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30)); // metadata filter

// Cloud mask based on MSK_CLDPRB layer
function mask_clouds(image) {
  // Select the MSK_CLDPRB band.
  var cloudProbability = image.select('MSK_CLDPRB');

  // Create a mask for clouds using a threshold.
  var cloudMask = cloudProbability.lt(50);

  // Update the image mask using the cloud mask.
  return image.updateMask(cloudMask);
}

// Apply on image collection
var S2_masked = S2.map(mask_clouds);

// Add True color RGB composites to the map
Map.addLayer(S2.first(),
            {min:0,max:3000,bands:"B4,B3,B2"}, 
            "RGB");
Map.addLayer(S2_masked.first(),
            {min:0,max:3000,bands:"B4,B3,B2"}, 
            "RGB masked");

Porovnání dat před a po použití masky oblačnosti

Kdybychom se chtěli podívat znovu i na jiné scény, převedeme jako v předchozím příkladu ImageCollection na List. Správně by se měly maskovat i stíny tvořené oblačností, ale to si když tak zkuste později samostatně.

Výpočet NBR indexu

Index NBR, neboli Normalized Burned Ratio, je index určený ke zvýraznění spálenišť v rozsáhlých požárních oblastech. Vzoreček pro NBR je následující:

NBR = (NIR - SWIR)/(NIR + SWIR)

NBR index pro celou ImageCollection spočítáme přidáním následujícího kódu.

// Compute NBR index
function add_NBR(image) {
  return image.addBands(image.normalizedDifference(['B8', 'B12']).rename('NBR'));
}
var S2_nbr = S2_masked.map(add_NBR);
print(S2_nbr);

Následně si ho můžeme přidat i do mapového okna. Co se týče oblačnosti, tak před požárem nejlépe vychází scéna z 23. 4. 2022 (tj. druhá scéna v ImageCollection), a po požáru pak nejspíše scéna z 25. 9. 2022 (tj. pátá scéna v ImageCollection). Proto si budeme znázorňovat tyto dvě.

// Convert ImageCollection to List
var S2_nbr_list = S2_nbr.toList(S2_nbr.size());

// Add layers to map
Map.addLayer(ee.Image(S2_nbr_list.get(1)),
            {min:0,max:3000,bands:"B4,B3,B2"}, 
            "RGB pre-fire");
Map.addLayer(ee.Image(S2_nbr_list.get(4)),
            {min:0,max:3000,bands:"B4,B3,B2"}, 
            "RGB post-fire");
Map.addLayer(ee.Image(S2_nbr_list.get(1)),
            {min: -1, max: 1, bands: ["NBR"]}, 
            "NBR pre-fire");
Map.addLayer(ee.Image(S2_nbr_list.get(4)),
            {min: -1, max: 1, bands: ["NBR"]}, 
            "NBR post-fire");

Na obrázku níže můžeme porovnat NBR index před a po požáru. Tmavé plochy znázorňují místa, kde k požáru došlo.

Porovnání NBR indexu před a po požáru

Dále si vytvoříme graf znázorňující časovou řadu NBR. Časovou řadu budeme ale znázorňovat jen pro území, kde k požáru opravdu došlo. Proto si toto území nebo jeho část vyznačíme polygonem, který vytvoříme pomocí nástroje Nakreslit tvar.

Nyní pomocí následujícího kódu vytvoříme graf.

// Create Time series chart
var TSChart = ui.Chart.image.series({
    imageCollection: S2_nbr.select(['NBR']),
    region: polygon,
    reducer: ee.Reducer.mean(),
    scale: 10,
}).setOptions({
    title: 'NBR time-series of the burned areas'
});

print(TSChart);

Vizualizace spálenišť

Dále si spáleniště nějak vhodně vizualizujeme. Lépe než kombinace B4-B3-B2 se pro spáleniště jeví kombinace B12-B11-B9 (tedy kombinace SWIR pásem a pásma Water vapour). Zároveň si vyselektujeme druhou a pátou scénu z ImageCollection do vlastních proměnných.

// Select pre-fire and post-fire images
var pre_fire_image = ee.Image(S2_nbr_list.get(1));
var post_fire_image = ee.Image(S2_nbr_list.get(4));

Map.addLayer(pre_fire_image,
            {min:0,max:3000,bands:"B12,B11,B9"}, 
            "pre_fire_image");
Map.addLayer(pre_fire_image,
            {min:-1,max:1,bands:"NBR"}, 
            "pre_fire_NBR");
Map.addLayer(post_fire_image,
            {min:0,max:3000,bands:"B12,B11,B9"}, 
            "post_fire_image");
Map.addLayer(post_fire_image,
            {min:-1,max:1,bands:"NBR"}, 
            "post_fire_NBR");

Zvýraznění spáleniště pomocí RGB kombinace B12-B11-B9

Dále si spočítáme tzv. Burn severity, což je veličina udávající závažnost poškození požárem, a vypočítá se podle následujícího vzorce:

dNBR = NBRPrefire - NBRPostfire

Čím vyšší hodnota, tím větší poškození požárem.

dNBR vypočítáme v Google Earth Engine pomocí následující kódu. Zároveň si i vyznačíme místa se střední až vysokou hodnotou poškození.

// Calculate the dNBR
var dNBR = pre_fire_image.select('NBR').subtract(post_fire_image.select('NBR')).rename('dNBR');

// High and very high burned severity
var modAndHighBS = dNBR.gt(0.27);

Map.addLayer(dNBR,
            {min:-1,max:1,bands:"dNBR"}, 
            "dNBR Image");
Map.addLayer(modAndHighBS,
            {min:-1,max:1,bands:"dNBR"}, 
            "Moderate and high burned severity");

Na závěr si můžeme vytvořit jednoduchou Swipe aplikaci pro přehledné znázorňování stavu před a po požáru.

// Application
// Create two map windows
var MapLeft = ui.Map();
var MapRight = ui.Map();

MapLeft.addLayer(pre_fire_image,
            {min:0,max:3000,bands:"B12,B11,B9"}, 
            "pre_fire_image");

MapLeft.addLayer(pre_fire_image,
            {min:-1,max:1,bands:"NBR"}, 
            "pre_fire_NBR");

MapRight.addLayer(post_fire_image,
            {min:0,max:3000,bands:"B12,B11,B9"}, 
            "post_fire_image");

MapRight.addLayer(post_fire_image,
            {min:-1,max:1,bands:"NBR"}, 
            "post_fire_NBR");

// Create a swipe option
var splitPanel = ui.SplitPanel({
  firstPanel: MapLeft,
  secondPanel: MapRight,
  orientation: "horizontal",
  wipe: true
})

// Clear map and add the swipe
ui.root.clear()
ui.root.add(splitPanel)

// Link the map windows
ui.Map.Linker([MapLeft,MapRight]);

// Center it around our selected point
MapRight.centerObject(point, 10)

MapLeft.setControlVisibility({all: true});
// MapRight.setControlVisibility({all: true});

// Add map description
var leftDescription = ui.Panel({
  widgets: [  
    ui.Label({value: "Pre-fire",
style: {textAlign:'centasdsadered',fontSize: '15px', color: '484848', fontWeight: 'bold'}})],
  style: {position: 'bottom-left'},
  layout: null,
});
MapLeft.add(leftDescription);

var rightDescription = ui.Panel({
  widgets: [  
    ui.Label({value: "Post-fire",
style: {textAlign:'centasdsadered',fontSize: '15px', color: '484848', fontWeight: 'bold'}})],
  style: {position: 'bottom-right'},
  layout: null,
});
MapRight.add(rightDescription);

// Create the panel widget
var panel = ui.Panel()
panel.style().set({ width: '350px',position: "middle-left"});
ui.root.insert(1,panel)

var title = ui.Label({value: "Detekce spálenišť pomocí GEE",
style: {textAlign:'centasdsadered',fontSize: '21px', color: '484848', fontWeight: 'bold'}});

panel.add(title)

var infotext1 = ui.Label({value:
"Aplikace znázorňující stav Českého Švýcarska před a po požáru v léte 2022." 
,style: {fontSize: '13px', color: '484848',textAlign:'justify'}});
panel.add(infotext1)

Inspirace pro toto cvičení čerpána z: GEE Workshop at Copernicus Forum 2022