Popis zpracování dat
Tento webový portál nabízí tři samostatné mapové aplikace – Trendy lesů, Historické snímky Landsat a Aktuální snímky Sentinel-2, které využívají shodná vstupní data – husté časové řady pozorování optických satelitních systémů Landsat 5 a 8 (provozováno NASA a USGS) a Sentinel-2 (provozováno Evropskou vesmírnou agenturou - ESA), liší se pouze ve způsobu vizualizace. Oproti např. leteckým datům mají satelitní pozorování nespornou výhodou v tom, že zájmové území mapují systematicky, typicky v řádu dní až týdnů, a to v širokém rozsahu vlnových délek od viditelného až po infračervené záření (mluvíme o tzv. multispektrálních senzorech). K dispozici tak máme bohatou časovou řadu pozorování, která jde v případě systému Landsat až do první poloviny 70.let (pro ČR existují první systematicky snímaná data od jara 1984) a v případě dat Sentinel-2 do roku 2015. To umožňuje získat ucelený pohled nejen o aktuálním stavu porostů, ale i o trendech jejich zdravotního stavu. Nevýhodou je naopak jejich typicky nižší prostorové rozlišení v řádu 10–30 m, které neumožňuje pracovat na úrovni jednotlivých stromů, ale spíše lesních celků. Problematická až donedávné doby byla náročnost zpracování takto velkého objemu dat a to i s ohledem k tomu, že značná část pozorování byla nějakým způsobem ovlivněna oblačností, která se na snímcích vyskytovala v momentě jejich pořízení.
Při interpretaci satelitních dat jsme typicky pracovali s jednotlivými snímky pořízenými v konkrétní čas, což neumožňovalo využít plného potenciálu bohatých časových řad satelitních snímků a jejich dlouhé historie. Toto se změnilo s nástupem technologie tzv. cloudového uložení a zpracování dat, kdy jsou veškerá data uložena v ucelených celcích na vzdálených výpočetních serverech a umožňují tak efektivní zpracování velkého objemu dat. Příkladem takovéto služby pro data dálkového průzkumu Země je nástroj Google Earth Engine, na kterém jsou založeny popisované mapové aplikace.
Cloudové zpracování satelitních dat v prostředí Google Earth Engine
Google Earth Engine (GEE) je nástroj firmy Google pro efektivní zpracování heterogenních datových sad produkovaných pestrou sadou senzorů DPZ a z nich odvozených produktů. Kromě zdrojových dat - snímků různých družicových systémů (Landsat, Sentinel-1,-2,-3,-5, MODIS) obsahuje rovněž datové sady o klimatu (satelitní pozorování, simulace klimatických modelů), krajinném pokryvu (např. databázi CORINE Land Cover) či topografii (digitální modely terénu). Veškeré popisované datové sady mají globální pokryv a jsou v prostředí GEE uloženy ve formě tzv. ImageCollection – balíku všech dostupných pozorování v čase. Zpracování a vizualizace těchto ImageCollections je možná pomocí jednoduchého webového editoru kódu programovacího jazyka JavaScript, který nabízí množství funkcí pro zpracování dat – typicky filtraci (odstraň snímky, které nevyhovují kritériím), maskování (odstraň části obrazu, které nevyhovují zadaným kritériím - typicky množství oblačnosti), redukci (ze sady snímků vyprodukuj jeden pro zobrazení na mapě), či publikaci (vlastní API, které nabízí různé nástroje pro tvorbu mapových aplikací). To umožňuje globální analýzy časových řad satelitních snímků, které by dříve nebylo možné provést. Jako příklad pro lesnictví může sloužit aplikace Global Forest Change.
Ukázka vývojového prostředí Google Earth Engine, který byl použit pro kompletní zpracování a publikaci mapových aplikací
Pro úspěšnou interpretaci satelitních dat většího území (zde Česká republika) je nutné splnit následující kritéria: 1) pracovat s kvalitními vstupními daty - v případě optických satelitních dat se snímky bez vlivu oblačnosti, které byly pořízeny ve vhodný čas snímkování a 2) zvolit odpovídající algoritmus pro zpracování těchto kvalitních vstupních dat. Způsob výpočtu bezoblačných mozaik pro časovou řadu satelitních pozorování senzorů Landsat a Sentinel-2 a způsob jejich interpretace pro vizualizaci trendů zdravotního stavu je popsán dále.
Tvorba bezoblačných mozaik
Satelitní snímky Landsat jsou pořizovány jednou za 16 dní, snímky dvojice družic Sentinel-2 poté jednou za pět dní. Vzhledem k tomu, že oba senzory jsou optické, jsou snímky často ovlivněny aktuálním stavem atmosféry (množství aerosolů a vodní páry) a přítomností mraků či stínů mraky vržených. Rovněž vegetace má své přirozené fenologické trendy, které jsou výrazné obzvláště u listnatých porostů, které limitují datum pořízení typicky na vegetační sezónu (zde 15.5. – 15.9.). To vše výrazně omezuje množství reálně využitelných snímků. U staršího satelitního systému Landsat 5 (1984–2012) je datový archiv neúplný a neumožňuje tak produkovat kvalitní mozaiky pro celé území ČR v ročním intervalu. Z výše uvedeného jsme se rozhodli produkovat 1) bezoblačné mozaiky ze systému Landsat v tříletém kroku (data od roku 1984 do 2020) pro studium dlouhodobých trendů zdravotního stavu a 2) ze systému Sentinel-2 v ročním kroku (data od roku 2015) pro hodnocení krátkodobého trendu zdravotního stavu posledních tří let.
Toho je v prostředí GEE dosaženo následujícím způsobem (následuje ukázka kódu pro jednu mozaiku):
var l5_tm = ee.ImageCollection("LANDSAT/LT05/C01/T1_SR").filterMetadata('CLOUD_COVER_LAND', 'less_than', 50).filterMetadata('IMAGE_QUALITY', 'greater_than', 5).filter(ee.Filter.dayOfYear(150, 250));
Tento kód definuje datovou sadu (ImageCollection) jako atmosféricky korigované snímky senzoru Landsat-5, které mají méně než 50% oblačnost, interní indikátor kvality je vyšší než 5 a které jsou pořízeny vždy ve vegetační sezóně.
function mask_tm_l4_l5(image) {
var cloudShadowBitMask = (1 << 3);
var cloudsBitMask = (1 << 5);
var qa = image.select('pixel_qa');
var qa2 = image.select('B1');
var qa3 = image.select('B3');
var qa4 = image.select('B4');
var qa5 = image.select('sr_atmos_opacity');
var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
.and(qa.bitwiseAnd(cloudsBitMask).eq(0)).and(qa2.gt(100)).and(qa2.lt(1500)).and(qa3.gt(100)).and(qa4.gt(500)).and(qa4.lt(8000)).and(qa5.lt(300));
return image.updateMask(mask);
}
Tento kód definuje funkci pro maskování pixelů ve filtrované ImageCollection snímků Landsat-5 pomocí aplikace kombinované masky, která využívá masky oblačnosti každého snímku, očekávatelných hodnot odrazivosti v několika spektrálních kanálech a maximální množství aerosolů v atmosféře.
var l5_tm_mask = l5_tm.map(mask_tm_l4_l5);
Tento kód aplikuje kombinovanou masku, která odstraní nekvalitní pixely z následných analýz
var l5_hq_1984_1986 = l5_tm_mask.filterDate('1984-05-15', '1986-09-15').median().clip(NUTS1);
Tento kód vypočítá finální bezoblačnou mozaiku ze sady filtrovaných a maskovaných snímků jako prostou mediánovou hodnotu odrazivosti pixelů mezi 15.5. a 15.9. vybraných let.
Výsledkem jsou sady tříletých (v případě dat systému Landsat) a ročních (v případě dat systému Sentinel-2) bezoblačných mozaik zobrazujících odrazivost povrchu v několika kanálech viditelného (400–700nm), blízkého infračerveného (700–1600 nm) a středního infračerveného záření (1600–2500 nm). Ty je možno mezi sebou jednoduchým způsobem vizuálně porovnat pro různé formy vizualizace, k čemuž slouží zde prezentované aplikace Historické snímky Landsat a Aktuální snímky Sentinel-2, a dále hodnotit krátkodobé a dlouhodobé trendy zdravotního stavu lesů.
Poznámka na závěr: Zejména u staršího satelitního systému Landsat-5, jenž obsahuje výrazně méně snímků pro zpracování kvalitních bezoblačných mozaik, se může stát, že i tříletá syntéza pozorování pro určité území neumožní získat bezoblačná data. V tomto případě je oblast zamaskována, toto není chyba ale vlastnost dat. Výrazně kvalitnější pozorování získáváme pro data Landsat-8 (2014 a dále) a Sentinel-2 (od roku 2017). Věříme, že i přesto jsou historická satelitní pozorování, které jdou téměř 40 let do historie, hodnotná.
Hodnocení krátkodobých a dlouhodobých trendů zdravotního stavu lesů
Kvalitní bezoblačné mozaiky satelitních multispektrálních dat v daném časovém kroku je možno interpretovat jako trendy zdravotního stavu lesů následujícím způsobem. Nejprve musíme zvolit vhodnou interpretaci vstupních hodnot odrazivostí v jednotlivých obrazových kanálech, která by byla citlivá na sledovaný fenomén – zde množství vody v korunách stromů. Z předchozích studií, zejména certifikované metodiky „Hodnocení zdravotního stavu lesních porostů v České republice pomocí satelitních dat Sentinel-2“ víme, že touto dostatečně citlivou obrazovou transformací je tzv. Wetness komponenta transformace obrazu Tasseled Cap. Jedná se o jednoduchou matematickou kombinaci všech vstupních obrazových kanálů, která využívá empiricky odvozené koeficienty z předchozích studií. Hodnocení krátkodobých a dlouhodobých trendů zdravotního stavu je založeno na hodnocení sklonu přímky lineárního modelu y = ax+b (komponenta a) u časové řady pozorování (osa x – roky pozorování, osa y – hodnota Wetness). Kladná hodnota koeficientu a značí nárůst vitality (množství vody v pletivech). Mírný pokles souvisí s přirozenými fenologickými trendy růstu porostu, zatímco prudký pokles indikuje rozpad porostu a těžbu. V aplikaci prezentujeme dva typy trendů – dlouhodobý trend 1984–2020 odvozený z dat Landsat a krátkodobý trend 2018–2020 z dat Sentinel-2. Obě tyto vrstvy zobrazují jiný typ informace – zatímco dlouhodobé trendy jsou vhodné pro studium např. regenerace porostů v čase a zachycují průměrný trend posledních 36 let, krátkodobé trendy zachycují změnu zdravotního stavu lesů posledních tří let, tedy v současných podmínkách primárně kalamitní kůrovcovou těžbu, trend chřadnutí porostů a regeneraci na plochách nedávno vytěžených.
Ukázka vývoje trendu zdravotního stavu lesního porostu od roku 1984 do současnosti. Na počátku 80. let byla plocha po nedávné těžbě opět zalesněna (nárůst hodnoty zdravotního stavu – množství vody v porostu), tento trend se však změnil mezi lety 2017 až 2020, kdy došlo k napadení porostu lýkožroutem smrkovým a prudkému poklesu zdravotního stavu. Zatímco dlouhodobý trend má kladnou hodnotu značící celkový trend nárůstu biomasy a regenerace původní holiny, krátkodobý trend posledních 3 let indikuje rozpad porostu.
Analýzy trendů Wetness komponenty a její vizualizace v prostředí GEE je provedena následovně:
function add_vi_tm_l4_l5 (img) {
var TCwet = img.expression('1*(0.1509*float(b("B1")) + 0.1973*float(b("B2"))+ 0.3279*float(b("B3")) + 0.3406*float(b("B4")) - 0.7112*float(b("B5")) - 0.4572*float(b("B7")))').rename(['tcw']);
return img.addBands([TCwet]);
}
Funkce, která přidá do bezoblačné mozaiky novou vrstvu – hodnoty komponenty Wetness.
var landsat_wetness_collection =
ee.ImageCollection([l5_1984.select('constant','tcw'),l5_1987.select('constant','tcw'),l5_1990.select('constant','tcw'),l5_1993.select('constant','tcw'),l5_1996.select('constant','tcw'),l5_1999.select('constant','tcw'),l5_2002.select('constant','tcw'),l5_2005.select('constant','tcw'),l5_2008.select('constant','tcw'),l5_2011.select('constant','tcw'),l8_2014.select('constant','tcw'),l8_2017.select('constant','tcw'),l8_2020.select('constant','tcw')])
var sentinel_wetness_collection = ee.ImageCollection([s2_2018_tcw.select('constant','tcw'),s2_2019_tcw.select('constant','tcw'),s2_2020_tcw.select('constant','tcw')]);
Definice nových ImageCollection pro bezoblačné mozaiky Landsat a Sentinel-2 obsahujících komponenty Wetness (kanál tcw) pro daný rok (kanál constant).
var fit_landsat = landsat_wetness_collection.reduce(ee.Reducer.linearFit());
var fit_sentinel = sentinel_wetness_collection.reduce(ee.Reducer.linearFit());
Výpočet lineárního modelu y = ax + b pro časovou řadu komponenty Wetness pro data Landsat a Sentinel-2.
mapPanel.addLayer(fit_landsat.select('scale').mask(olil).clip(NUTS1),{min: -25, max: 25, palette: ['#ff0000', 'yellow' , '#1d910f','#1d910f']},"Trendy zdravotního stavu 1984-2020",0,1);
mapPanel.addLayer(fit_sentinel.select('scale').mask(olil).clip(NUTS1),{min: -200, max: 200, palette: ['#ff0000', 'yellow' , '#1d910f','#1d910f']},"Trendy zdravotního stavu 2018-2020",1,1);
Vizualizace kanálu scale (komponent a – sklon přímky lineárního modelu) pro data Landsat a Sentinel-2.
Vizuálně je možno krátkodobé a dlouhodobé trendy hodnotit v aplikaci Trendy lesů, která také nabízí možnost výpočtu trendu pro uživatelem zvolený pixel. Po kliknutí na vybrané místo v hlavním mapovém okně se trend vykreslí dynamicky v grafu, který je vygenerován v pravé straně mapové aplikace. Graf je možno zvětšit a exportovat.
Poznámka na závěr: V grafu je možné pro některé oblasti získat hodnoty v daném roce –9999, toto jsou oblasti s vysokou mírou oblačnosti, které nebyly pro výpočet trendu použity.
Pokud jste dočetli až sem a máte nějaké další otázky nebo byste chtěli získat zdrojové JavaScript soubory pro GEE, kontaktujte nás.