NmetS - Metodiky schválené příslušným orgánem státní správy, do jehož kompetence daná problematika spadá


Metodika je výsledkem projektu: Lesní mikroklima v čase a prostoru: reálné dopady změny klimatu na vybraná chráněná území


Číslo projektu: SS06010011

Začátek projektu: 04/2023

Konec projektu: 03/2026

Řešitel projektu: Botanický ústav AV ČR, v. v. i., Zámek 1, Průhonice, 25243, Česká republika

Důvěrnost a dostupnost: S – Úplné a pravdivé údaje o projektu nepodléhající ochraně podle zvláštních právních předpisů.

Verze: 1.1

Vydaná: 2024 Aktualizace: 22.5.2025 (Update myClim na verzi 1.40)

Testováno s R 4.4.1, myClim 1.40, RStudio 2025.05.0






Tento projekt je financován se státní podporou Technologické agentury ČR a Ministerstva životního prostředí ČR v rámci Programu Prostředí pro život.“

Informace o autorském týmu:

RNDr. Josef Brůna, Ph.D., Řešitel

Mgr. Tereza Klinerová, členka řešitelského týmu

Mgr. Martin Kopecký, Ph.D., člen řešitelského týmu

Mgr. Martin Macek, Ph.D., člen řešitelského týmu

Mgr. Matěj Man, člen řešitelského týmu

Mgr. Anna Růžičková, členka řešitelského týmu

doc. Ing. Jan Wild, Ph.D., člen řešitelského týmu
























Tento projekt je financován se státní podporou Technologické agentury ČR a Ministerstva životního prostředí ČR v rámci Programu Prostředí pro život.“

Úvod

Mikroklimatická měření jsou dnes díky vývoji automatizovaných autonomních měřicích přístrojů široce dostupná, ale jejich zpracování může být náročné vzhledem k velkým objemům dat a nestandardizovaným výstupním formátům dat i nestandardizovaným způsobům měření v závislosti na konkrétní výzkumné otázce. Proto bývají v ochraně přírody (např. při zhodnocení vlivu klimatu na předměty ochrany) využívány pouze ojediněle. Cílem této metodiky je usnadnit využití mikroklimatických dat v praxi ochrany přírody a zajistit reprodukovatelnost výsledků, a to za pomoci praktických ukázek a ověřených postupů. Cílovou skupinou uživatelů jsou především pracovníci správ NP a CHKO, ale též pracovníci Agentury ochrany přírody a krajiny.

První část metodiky se věnuje odvození biologicky relevantních proměnných zejména na základě časových řad z mikroklimatických stanic TOMST, ale i z jiných zdrojů, tak, aby bylo možné využít i další kompatibilní data. Hlavní motivací je fakt, že přestože klasická meteorologická data již standardizovaná jsou a nabízí řadu srovnatelných proměnných dostupných z různých zdrojů, tak data z kompaktních mikroklimatických stanic zdaleka tak standardizovaná nejsou a analýza dat získaných z jednotlivých stanic může být na první pohled náročná. Na příkladech založených na reálných datech ukážeme hlavní problémy při zpracování dat a doporučíme osvědčené postupy. Druhá část metodiky se věnuje praktickému využití odvozených proměnných pro vyhodnocení dat a také využití mikroklimatických a klimatických map.

Metodika popisuje postupy pro získání a zpracování mikroklimatických měření z autonomních dataloggerů (instalace, správa dat, kontrola kvality dat, zpracování a vizualizace naměřených hodnot, základní analýzy dat) kompletně v open-source programovém prostředí R. Věříme, že absence vazeb na komerční software usnadní její využití. Výhodou zpracování dat v prostředí R je rovněž reprodukovatelnost postupů a výsledků, protože veškeré operace s daty zůstávají zdokumentovány ve zdrojovém kódu.

Jak pracovat s metodikou

Tato metodika existuje jako samostatný dokument (ve formátu .pdf), ale je též celá dostupná online včetně všech zdrojových souborů použitých v příkladech zpracování dat. Celý projekt pro program RStudio - https://www.rstudio.org/ a jednotlivé úlohy v něm uvedené tak mohou uživatelé spustit a případně upravit kód pro vlastní potřebu na svém zařízení.

Formát .Rmd kombinuje vlastní text dokumentu a spustitelný kód v jazyce R, pro vlastní analýzy lze jednotlivé části kódu přenést do nového skriptu ve formátu R. Zároveň jsou dostupné samotné kódy jednotlivých částí metodiky pro prostředí R bez okolního textu (formát .R).

Vzhledem k dynamické povaze knihoven pro prostředí R bude online verze metodiky průběžně aktualizována po dobu implementace projektu, aby byla kompatibilní s aktuálním vývojem knihovny myClim (Man et al. 2023), případně dalších důležitých knihoven a zdrojů dat. Drobné odchylky se i tak mohou vyskytnout a způsobit chyby v běhu kódu, v případě potřeby se na nás neváhejte obrátit.

Online verze se všemi vzorovými daty je dostupná zde: https://gitlab.ibot.cas.cz/vysledky/metodika-mikroklima-ppz

Odtud je třeba si stáhnout z webu pomocí tlačítka Code - zde zvolte zip. Přímý odkaz je: https://gitlab.ibot.cas.cz/vysledky/metodika-mikroklima-ppz/-/archive/main/metodika-mikroklima-ppz-main.zip

Stažení repozitáře ve formátu ZIP.
Stažení repozitáře ve formátu ZIP.

Následně obsah extrahovat na disku a otevřít soubor metodika-mikroklima-ppz.Rproj v programu RStudio.

V okně RStudia lze pak otevírat jednotlivé soubory Rmd (Obrázek 2, v části soubory vpravo dole) a postupně pouštět kód v hlavním okně pomocí zelených šipek v pravém rohu každého okna kódu:

Doporučujeme zobrazit s formátováním pomocí tlačítka Visual. V části s proměnnými v paměti lze najít aktuálně načtené proměnné a prozkoumat jejich obsah. Orientaci usnadňuje osnova, kterou lze zapnout v pravé části hlavního okna tlačítkem Outline. Ve spodní části vlevo je výstup konzole, kde je vidět průběh jednotlivých příkazů a mohou se zde zobrazit některé méně podstatné výstupy, případně hlášky o chybách.

Prostředí R studia s popisky částí, barevný styl RStudia i rozložení oken se může lišit.
Prostředí R studia s popisky částí, barevný styl RStudia i rozložení oken se může lišit.

Hlavní výstupy (obrázky) jsou pak vždy pod každým blokem kódu (Obrázek 3). Jednotlivé kapitoly je třeba pouštět postupně odshora, následující bloky kódu často využívají proměnné vytvořené v předchozích blocích.

Příklad výstupu jednoho bloku kódu, červeně označeno tlačítko spuštění bloku kódu.
Příklad výstupu jednoho bloku kódu, červeně označeno tlačítko spuštění bloku kódu.

Nastavení prostředí a funkcí pro vykreslování tabulek

1. Workflow

V úvodní části se zaměříme na získání relevantních mikroklimatických dat a jejich přípravu pro vstup do následných analýz. Projdeme postupně celý proces od instalace čidel v terénu, přes stahování dat až po jejich kontrolu, čistění, kalibraci a spojování časových řad z více odečtů. Dále ukážeme, jak vést záznamy o historii čidel (např. odečty, výměny čidel v rámci lokality) a nastavit vhodnou strukturu pro skladování stažených dat.

Důraz klademe na práci s čidly systému TMS (TOMST) - https://tomst.com/web/cz/systemy/tms/, v současnosti celosvětově široce využívanými pro mikroklimatická měření (Lembrechts at al., 2022), ale metodiku je možné zobecnit i na jiné zdroje mikroklimatických dat - např. z instrumentálního měření systémy HOBO (Onset), iButton (Maxim/Dallas), Minikin (EMS Brno), EnvLoggers (ElectricBlue) a dalšími. Pro práci s mikroklimatickými daty využíváme R knihovnu myClim (Man et al., 2023), ve které jsme shromáždili zkušenosti z mnoha let práce s mikroklimatickými daty.

1.1. Instalace TMS čidel v terénu

Správná instalace TMS čidel v terénu je základním předpokladem pro získání relevantních mikroklimatických dat. Předchází jí výběr reprezentativních lokalit pro umístění čidel, který se odvíjí od konkrétní výzkumné otázky. Dodržení jednotných postupů zajistí srovnatelnost naměřených dat. Detailní informace shrnuje Metodika měření mikroklimatu pomocí mikroklimatických stanic systému TMS (Brůna et al., 2021), zde uvádíme pouze stručné shrnutí.

Thermologger

Čidlo Thermologger (TOMST) je autonomní datalogger měřící teplotu vzduchu. Osazuje se pomocí metrického závitu M5 do bílého stínění, chránícím čidlo před přímým slunečním zářením, zpravidla ve svislé poloze ve výšce 2 m nad povrchem země (nevyžaduje-li povaha výzkumné otázky jiné umístění). V lesním porostu lze umístit na severní stranu kmene stromu (Obr. 4 A a B) či na levnou a odolnou zahradnickou poplastovanou tyčku.

A) a B) Thermologger osazený na severní stranu kmene stromu ve výšce 2 m nad zemí. C) TMS-4 datalogger - umístění signalizační LED diody (při instalaci v terénu směřuje na jih) a správná hloubka zanoření spodní (zelené) části dataloggeru v půdě. D) Výzkumná plocha osazená TMS-4 dataloggerem a Thermologgerem umístěným na kmeni stromu. E) TMS-4 datalogger v ochranné kleci. F) Výzkumná plocha osazená TMS-4 dataloggerem a Thermologgerem umístěným na tyči.
A) a B) Thermologger osazený na severní stranu kmene stromu ve výšce 2 m nad zemí. C) TMS-4 datalogger - umístění signalizační LED diody (při instalaci v terénu směřuje na jih) a správná hloubka zanoření spodní (zelené) části dataloggeru v půdě. D) Výzkumná plocha osazená TMS-4 dataloggerem a Thermologgerem umístěným na kmeni stromu. E) TMS-4 datalogger v ochranné kleci. F) Výzkumná plocha osazená TMS-4 dataloggerem a Thermologgerem umístěným na tyči.

TMS-4 datalogger

TMS-4 datalogger (TOMST) byl navržen pro monitoring mikroklimatických proměnných relevantních pro bylinné patro (Wild et al., 2019). Datalogger TMS-4 měří teplotu ve třech pozicích - teplotu vzduchu v 15 cm nad zemí, teplotu povrchu půdy, teplotu půdy 8 cm pod povrchem půdy a vlhkost půdy. Čidlo se osazuje zpravidla do ochranné klece, která brání poškození čidla lesní zvěří a pádem větví (Obr. 4 E), signalizační LED diodou na jih, správnou hloubku zanoření spodní (zelené) části dataloggeru se senzorem půdní vlhkosti a teploty v půdě ukazuje Obr. 4 C. Pokud se čidlo instaluje v kamenité půdě, je pro vyhloubení otvoru v zemi vhodné použít pomocný kovový nástroj, aby nedošlo k poškrábání nebo úplnému poničení senzoru půdní vlhkosti (obvod na plošném spoji). Důležitá je kontrola, zda je celá spodní část čidla v dobrém kontaktu s půdou (čidlo drží ve vyhloubeném otvoru pevně, neviklá se) a případné vyplnění vzduchových kapes půdním substrátem. Ke stínění čidla před přímým slunečním zářením je možné použít jednoduché nebo dvojité stínění. Dvě varianty umístění čidel na výzkumné ploše ukazují Obr. 4 D (Thermologger na stromě) a Obr. 4 F (Thermologger na poplastované tyčce).

1.2. Stažení a správa dat

Rozhraní pro obsluhu čidel (Obr. 5) - stahování dat a základní pohled na data zajišťuje software Lolly od výrobce čidel TMS. Při první instalaci je třeba mít k počítači připojen TMD adaptér. Kompletní návod pro obsluhu softwaru naleznete na stránkách výrobce: https://tomst.com/web/wp-content/uploads/2023/06/Lolly-software-Handbook.pdf (v angličtině).

Z hlediska použití v praxi poukazujeme především na:

  • možnost nastavení různých intervalů měření dat

  • možnost nastavení různých formátů datumu a času (dbát na kompatibilitu s myClim)

  • možnost vytvoření bookmarku (záložky)

  • možnost odečtu dat od konkrétního datumu nebo bookmarku

  • možnost zobrazení přehledového grafu naměřených dat

Rozhraní softwaru Lolly - nastavení parametrů často potřebných v praxi. Kompletní návod pro obsluhu softwaru je zde: tomst.com/web/wp-content/uploads/2023/06/Lolly-software-Handbook.pdf.
Rozhraní softwaru Lolly - nastavení parametrů často potřebných v praxi. Kompletní návod pro obsluhu softwaru je zde: tomst.com/web/wp-content/uploads/2023/06/Lolly-software-Handbook.pdf.

Časové řady s naměřenými hodnotami jsou po vyčtení dataloggeru uloženy do CSV souboru pojmenovaného textovým řetězcem:

“data_” + sériové číslo + rok_měsíc_den + pořadové číslo odečtu + “.csv”

Soubor obsahuje tabulková data ve formě hodnot oddělených středníky. Jednotlivé sloupce představují postupně: pořadové číslo pozorování, datum a čas měření v UTC, kód časového pásma, teplota senzoru T1, teplota senzoru T2, teplota senzoru T3, půdní vlhkost, indikátor otřesu, chybový příznak. Pro loggery TOMST TMS-4 jsou T1 půdní senzor (-8cm), T2 přízemní senzor (+2cm), T3 vzdušný senzor (+15cm). Pro TOMST Thermologger je teplota uložena v poli T1, v ostatních nevyužitých polích je uložena hodnota -200.

Dále je spolu s daty uložen i  soubor “command_*.csv” obsahující výpis průběhu odečtu, který pro další práci nebudeme potřebovat, nicméně jej doporučujeme archivovat pro případné řešení problémů s odečtem s výrobcem.

Stažená data je vhodné ukládat na cloudové úložiště pro on-line zálohu dat (možnost nastavit cestu pro ukládání dat do synchronizované složky přímo v softwaru Lolly). Data je vhodné ukládat do složek tak, aby bylo následně snadné propojit terénní zápisky s úložištěm dat. Např. pokud odečítáme data 1 x ročně, pak pojmenovat složky podle daného roku a v terénních zápiskách organizovat poznámky také s příznakem daného roku. Kopii stažených souborů archivujeme pro úsporu místa nejlépe v komprimovaném archivu (.zip), protože .csv soubory jsou datově náročné.

V případě potřeby změnit formát datumu a času zpětně u již odečtených dat je možné využít software Lolly (neoficiální návod viz Obr. 6 A-C).

A) Nejprve nastavíme na kartě Options požadovaný formát datumu kompatibilní se softwarem myClim. B) Poté přejdeme na kartu Info a klikneme na “Show Data”. C) Tím aktivujeme prohlížeč souborů, navigujeme do složky, kde máme uložená data určená k převodu na jiný formát data a času, a klikneme na “Fix time”. Takto převedeme všechny .csv soubory v dané složce na formát data a času nastavený na kartě Options. Převedené soubory budou mít nově předponu conv_.
A) Nejprve nastavíme na kartě Options požadovaný formát datumu kompatibilní se softwarem myClim. B) Poté přejdeme na kartu Info a klikneme na “Show Data”. C) Tím aktivujeme prohlížeč souborů, navigujeme do složky, kde máme uložená data určená k převodu na jiný formát data a času, a klikneme na “Fix time”. Takto převedeme všechny .csv soubory v dané složce na formát data a času nastavený na kartě Options. Převedené soubory budou mít nově předponu conv_.

Sběr dat v terénu je vhodné provádět jednou až dvakrát ročně (jaro, podzim), aby v případě poruchy čidla nedošlo k větší ztrátě dat. Pro rutinizaci a snazší koordinaci odečtových prací je možné využít „terénní manuál“ uvedený v Příloze 1 výše zmíněné Metodiky měření mikroklimatu pomocí mikroklimatických stanic systému TMS (Brůna et al., 2021). Je nutné zaznamenávat průběh odečtu a všechny nestandardní události (vytažení čidla, spadlé stínění, ID nově nainstalovaného čidla atp.). Tyto poznámky jsou nezbytné pro správné procesování dat:

  1. pro kontrolu, zda máme odečty ze všech čidel a lokalit,
  2. pro správnou interpretaci dat a vyhodnocení případných chyb v datech.

Časté problémy a poruchy zjistitelné přímo v terénu

  • špatný kontakt spodní části čidla TMS-4 s půdou

Pokud je čidlo v půdě volně, hrozí v důsledku špatného kontaktu s půdou nepřesné měření půdní vlhkosti (čidlo měří jen v malém objemu půdy těsně přiléhajícím k senzoru). V grafu se problém projevuje náhlými propady vlhkosti, či celkově nízkými hodnotami (pod 1000 jednotek). Čidlo je třeba přeinstalovat. Problém poznamenáme a při procesování naměřených hodnot dbáme v případě půdní vlhkosti zvýšené obezřetnosti, případně nahradíme půdní vlhkosti NA hodnotami.

  • vytažené TMS-4 čidlo

Pokud je TMS-4 čidlo nalezeno volně ležící na povrchu půdy (a není viditelně mechanicky poškozené), v softwaru Lolly zkontrolujeme funkčnost čidla (na přehledovém grafu nejsou patrné výpadky měření, konstantní ani zjevně nesmyslné hodnoty měřených proměnných). Pokud je čidlo funkční, opětovně jej nainstalujeme. Při zpracování dat je ovšem třeba dohledat, kdy k vytažení čidla z půdy došlo, což lze poznat podle náhlého vyrovnání teplot T1, T2 a T3 a prudkého poklesu hodnot půdní vlhkosti. V daném období pak musíme všechny datové řady nahradit NA hodnotami (teplotu T2 měřenou na povrchu půdy lze při liberálním přístupu zachovat), podrobněji níže. R knihovna myClim také nabízí funkci pro automatickou detekci dat naměřených na vytaženém čidle mc_prep_TMSoffsoil().

  • spadlé stínění

Zvláště čidla TMS bez ochrany před zvěří trpí ztrátou stínění. Je třeba vrátit stínění na čidlo. V takovém případě nejsou naměřené teploty srovnatelné s hodnotami naměřenými na stíněných čidlech (zejména teplota 15 cm nad zemí a teplota na povrchu půdy). Stínění brání přehřívání čidla při dopadu přímého slunečního záření na senzor.

  • výpadky senzorů

Rozpoznáme na přehledovém grafu v Lolly softwaru. Jde-li o jednorázové události, nemusí být funkčnost čidla omezena a můžeme jej ponechat na místě. Pokud jsou výpadky rozsáhlejší, čidlo vyměníme za nové. Pokud chybí celý časový úsek, je možné, že se jedná o chybu přenosu dat při vyčítání. V takovém případě opakujeme odečet. V případě, že je nastavené vyčítání od bookmarku, je nutné tuto volbu dočasně změnit a vyčíst čidlo celé nebo od konkrétního data, protože se pravděpodobně již vytvořil nový bookmark.

  • posun času

Při každém odečtu se synchronizuje čas v počítači a v čidle. V případě, že výrazně nesouhlasí interní čas na čidle a v počítači (tj. rozdíl více než 10 minut uváděný v Lolly v poli “Delta:” - Obrázek 5B), je třeba tuto skutečnost poznamenat v terénních zápiscích a zjistit pravděpodobnou příčinu. Může jít o důsledek zpožďujících se vnitřních hodin dataloggeru, např. u čidel se slábnoucí baterií. Může jít také o chybu v důsledku špatně nastaveného systémového času počítače buď při předchozím odečtu nebo při aktuálním odečtu (pozn. čidla TMS uvádějí čas UTC, od středoevropského času se tedy liší o 1 hodinu, resp. o 2 hodiny od letního času).

1.3. Vedení záznamů

Nutnou podmínkou pro přípravu relevantních vstupních dat jsou záznamy o umístění čidel na lokalitě a poznámky z odečtů v terénu, které slouží ke správné interpretaci a vyhodnocení dat. Při návštěvách lokalit je proto potřeba evidovat veškeré změny stavu čidel. Tyto informace jsou pak nedílnou součástí následného zpracování. Nejjednodušším případem takové záznamové tabulky pro naše vzorová data může být např. tato:

# Načtení terénních poznámek z Excelu
zaznamy <- read.xlsx("./data/zaznamyR.xlsx")
Ukázka terénních poznámek
locality_id logger_type 2021 2021.pozn 2022 2022.pozn 2023 2023.pozn
CZ2_MASTALE HOBO_RH 20024371 ok, 30min 20024371 ok, 30min 20024371 výměna membrány
CZ2_HODDUB HOBO_RH 21061738 ok, 30min 21061738 ok, 30min 21061738 ok, 30min
CZ2_CEPICKA HOBO_RH 21061739 ok, 30min 21061739 ok, 30min 21061739 ok, 30min
CZ2_LUZNICE HOBO_RH 21061759 ok, 30min 21061759 ok, 30min 21061759 ok, 30min
CZ2_KLINOV HOBO_RH 21061889 ok, 30min 21061889 ok, 30min 21061889 ok, 30min
CZ2_HRUSUT HOBO_RH 21061891 ok, 30min 21061891 ok, 30min 21061891 ok, 30min
CZ2_BARINY HOBO_RH 21067836 ok, 30min 21067836 ok, 30min 21067836 výměna membrány
CZ2_KRKLAB HOBO_RH 21067844 ok, otočeno větrem 21067844 ok, otočeno větrem 21067844 ok, 30min
CZ2_BUKOVEC HOBO_RH 21067850 ok, 30min 21067850 ok, 30min 21067850 ok, 30min
CZ2_HODDUB Thermo 91171001 ok 91171001 ok 91171001 ok
CZ2_LUZNICE Thermo 91171050 ok 91171050 1x error, ponecháno 91171050 ok
CZ2_CEPICKA Thermo 91192386 ok 91192386 ok 91192386 ok
CZ2_BUKOVEC Thermo 91201102 ok 91201102 ok 91201102 ok
CZ2_MASTALE Thermo 91201109 ok 91201109 ok 91201109 ok
CZ2_HRUSUT Thermo 91201114 ok 91201114 ok 91201114 ok
CZ2_BARINY Thermo 91201121 ok 91201121 ok 91201121 ok
CZ2_KLINOV Thermo 91201144 ok 91201144 ok 91201144 ok
CZ2_KRKLAB Thermo 91201145 ok 91201145 ok 91201145 ok
CS_26 TMS4 94183144 ok 94183144 vykopnuté, zapíchnuto 94183144 vytažené, rozkousané, ale odečteno
CZ2_CEPICKA TMS4 94192751 ok 94192751 ok 94192751 výpadky vlhkosti, -5:02:21

Pro každý odečet (v případě našich vzorových dat odečet 1 x ročně) je potřeba pro jednoznačnou identifikaci zaznamenat sériové číslo dataloggeru (ID) (vytištěné na tubusu), unikátní identifikátor lokality a další poznámky k odečtu. V našem příkladu jsou na lokalitách umístěny stále stejné dataloggery, to však nemusí být pravidlem - pokud dojde k poškození či výměně dataloggeru v průběhu času, změní se tak ID. Díky těmto záznamům budeme později schopní správně přiřadit datové řady k jednotlivým lokalitám, i když dojde k výměně dataloggeru. Rovněž důrazně nedoporučujeme přejmenování souborů, které jsou identifikované sériovým číslem dataloggeru, protože sériové číslo dataloggeru slouží i k přiřazení kalibračních dat. Z poznámek k odečtu je pak zřejmé, že při zpracování dat bude potřeba se zaměřit na dataloggery na lokalitách CS_26, kde došlo k vytažení dataloggeru ze země, což znehodnotilo měření, a dále se zaměříme na lokality CZ2_CEPICKA a CZ2_LUZNICE, kde máme poznámku o výpadku senzoru vlhkosti.

1.4. Načtení a čištění surových mikroklimatických dat

Ukázka přípravy dat pro další analýzy ze surových dat z terénu. Jako příklad uvádíme data naměřená na 10 lokalitách napříč ČR. Na všech lokalitách byly umístěny TOMST TMS-4 dataloggery měřící vzdušnou teplotu ve výšce 15 cm a 2 cm nad povrchem půdy, 8 cm pod povrchem půdy a vlhkost půdy ve svrchní vrstvě půdy (0 - 15 cm). Na 9 lokalitách byl umístěn také Thermodatalogger měřící vzdušnou teplotu 200 cm nad zemí a HOBO U23 datalogger měřící vzdušnou teplotu a relativní vlhkost vzduchu ve výšce 150 cm nad zemí. Časová řada obsahuje měření přibližně od dubna 2021 do října 2023, ale všechny dataloggery neměřily stejnou dobu, v datech jsou chybějící úseky.

# Seznam potřebných knihoven
required_packages <-
  c("myClim", "dplyr", "tidyr", "stringr", "openxlsx", "ggplot2")

# Instalace a načtení knihoven
suppressWarnings(lapply(required_packages, function(pkg) {
  if (!require(pkg, character.only = TRUE))
    install.packages(pkg, repos = "https://mirrors.nic.cz/R/")
  library(pkg, character.only = TRUE)
}))

# Úklid proměnných v prostředí R
rm(required_packages)

Máme nainstalované a načtené potřebné knihovny pro práci s mikroklimatickými daty a nyní načteme údaje o souborech odečtených v terénu. Pro zpracování dat v knihovně myClim je důležité správně přiřadit jednotlivé soubory k lokalitám, správně zadat očekávaný formát dat (typ dataloggeru, ze kterého data pochází) a dále definovat, v jakém formátu jsou v souborech uvedené datumy a časy (zda je to např. 2023.12.31 12:00 nebo 31.12.2023 12:00:00). Tomuto definičnímu souboru říkáme files_table, můžeme jej vytvořit v kódu z naší tabulky záznamů nebo ručně v tabulkovém editoru a pak načíst ve formátu CSV.

## sestavení files_table z výše načtené excel tabulky "zaznamy"
## v kombinaci s cestami ke staženým souborů z loggerů
## Získání cest ke staženým souborům
soubory <- list.files("./data/example_data_prep_calc",
                      recursive = T,
                      full.names = T)
## Extrakce ID pro 2 typy loggerů z názvů souborů (TMS, HOBO)
jmena <- basename(soubory) # názvy souborů
logger_id <-
  str_match(jmena, "data_\\s*(.*?)\\s*_") # extrakce TOMST TMS ID
hob <-
  is.na(logger_id[, 2]) # které loggery nmají TMS ID, to jsou HOBO
logger_id[hob] <-
  substr(jmena[hob], 1, 8) # extrakce ID HOBO loggerů

## Extrakce roku odečtu ze jména složky, ve které se data nachází
odecty <- basename(dirname(soubory))

## Nyní známe cestu, ID loggeru a rok odečtu.
ft <- data.frame(path = soubory,
                 logger_id = logger_id[, 2],
                 odecty = odecty)

## Chybí nám však zatím údaj o lokalitě, typu loggeru a formátu data
## tyto informace získáme z tabulky terénních záznamů
## potřebujeme klíč, které loggery byly ve kterém roce na které lokalitě
## to jsou sloupce 3,5,7
zaznamy.dlouha <- pivot_longer(zaznamy,
  cols = c(3, 5, 7),
  names_to = "odecty",
  values_to = "logger_id")

## sloupce s poznámkou nyní nepotřebujeme
## získáme tabulku s klíčem lokalita-logger_type-odečet-logger_id
zaznamy.dlouha <- zaznamy.dlouha[,-c(3:5)]

## vytvoříme "files_table"
## spojení klíče lokalit a cesty k souborům podle logger ID a roku odečtu
files_table <-
  merge(ft, zaznamy.dlouha, by = c("logger_id", "odecty"))

## z údajů o typu loggeru odvodíme myClim datový typ
data_format <- files_table$logger_type
data_format[data_format == "Thermo"]  <- "TOMST"
data_format[data_format == "TMS4"]    <- "TOMST"
data_format[data_format == "HOBO_RH"] <- "HOBO"
files_table$data_format <- data_format

## zkoumáním sloupců s datem v našich souborech z loggerů
## zjistíme že soubory z HOBO mají konzistentní formát data
## u TOMST souborů se vyskytují 4 formáty data
## formát datumu závisí na nastavení výstupu z odečítacího software
## čím více osob a počítačů se na stahování dat podílí, tím více možných formátů
date_format <- files_table$data_format
date_format[date_format != "HOBO"] <-
  "%d.%m.%Y %H:%M:%S@%d.%m.%Y %H:%M@%Y.%m.%d %H:%M@%d.%m.%Y"
date_format[date_format == "HOBO"] <- "%d.%m.%Y %H:%M:%S"
files_table$date_format <- date_format

## Ponecháme pouze potřebné sloupce
files_table <-
  files_table[, c("path", "locality_id", "data_format", "date_format")]

Soubor files_table.csv musí dodržet danou strukturu (viz tab. 2):

Ukázka struktury tabulky definující příslušnost souborů k lokalitám a datový typ souborů
path locality_id data_format date_format
./data/example_data_prep_calc/2021/20024371.txt CZ2_MASTALE HOBO %d.%m.%Y %H:%M:%S
./data/example_data_prep_calc/2022/20024371.txt CZ2_MASTALE HOBO %d.%m.%Y %H:%M:%S
./data/example_data_prep_calc/2023/20024371.txt CZ2_MASTALE HOBO %d.%m.%Y %H:%M:%S
./data/example_data_prep_calc/2021/21061738.txt CZ2_HODDUB HOBO %d.%m.%Y %H:%M:%S
./data/example_data_prep_calc/2022/21061738.txt CZ2_HODDUB HOBO %d.%m.%Y %H:%M:%S
./data/example_data_prep_calc/2023/21061738.txt CZ2_HODDUB HOBO %d.%m.%Y %H:%M:%S
./data/example_data_prep_calc/2021/21061739.txt CZ2_CEPICKA HOBO %d.%m.%Y %H:%M:%S
./data/example_data_prep_calc/2022/21061739.txt CZ2_CEPICKA HOBO %d.%m.%Y %H:%M:%S
./data/example_data_prep_calc/2023/21061739.txt CZ2_CEPICKA HOBO %d.%m.%Y %H:%M:%S
./data/example_data_prep_calc/2021/21061759.txt CZ2_LUZNICE HOBO %d.%m.%Y %H:%M:%S

Nyní načteme vlastní data z mikroklimatických stanic do jednoho objektu třídy myClim (raw format):

## Načtení mikroklimatických dat do myClim
## přiřazení loggerů na lokality podle "files_table"
micro.data <-  mc_read_data(files_table = files_table,
               silent = F,
               clean = F)

## Načtení pomocí files_table připravené ručně, uložené do .csv souboru
micro.data <-  mc_read_data(files_table = "./data/files_table.csv",
               silent = F,
               clean = F)

## Provedeme základní strojové čištění a kontrolu dat, konzistence 
## časového kroku, duplicity, posloupnost časové osy, chybějící hodnoty
micro.data.clean <- mc_prep_clean(micro.data, silent = T)

Pro základní přehled o datech si vypíšeme informace o rozsahu dat, kroku měření a chybách v datech (chybějící hodnoty, duplicity, špatné řazení).

## Vypsání informací o loggerech: rozsah dat, krok měření, 
## duplicity, pořadí, chybějící data
info.tms <- mc_info_clean(micro.data.clean)
Výsledek kontroly dat v objektu ‘micro.data.clean’ - rozsah měření, počty duplicitních, chybějících a chronologicky špatně řazených záznamů.
locality_id serial_number logger_name start_date end_date step duplicities missing disordered rounded
CS_26 94183144 TMS_1 2021-04-09 2021-11-03 900 0 0 0 FALSE
CS_26 94183144 TMS_2 2021-11-03 2022-10-19 900 0 0 0 FALSE
CS_26 94183144 TMS_3 2022-10-19 2023-04-06 900 0 0 0 FALSE
CZ2_BARINY 21067836 HOBO_U23-001A_1 2021-04-26 2021-11-08 1800 0 0 0 FALSE
CZ2_BARINY 21067836 HOBO_U23-001A_2 2021-11-08 2022-10-24 1800 0 0 0 FALSE
CZ2_BARINY 21067836 HOBO_U23-001A_3 2022-10-24 2023-10-31 1800 0 0 0 FALSE
CZ2_BARINY 91201121 Thermo_1 2021-04-29 2021-11-08 900 0 0 0 FALSE
CZ2_BARINY 91201121 Thermo_2 2021-11-08 2022-10-24 900 33602 0 1 FALSE
CZ2_BARINY 91201121 Thermo_3 2022-10-24 2023-10-31 900 0 0 0 FALSE
CZ2_BARINY 94196117 TMS_1 2021-04-29 2021-11-08 900 0 0 0 FALSE

Pracujeme s daty z období podzim 2021 až podzim 2023 naměřenými na 10 lokalitách napříč ČR. Jde o data ze třech různých odečtových akcí. Výpis protokolu ze strojové kontroly z knihovny myClim mc_info_clean() ukazuje znepokojivý počet duplicitních hodnot a hodnot v nesprávném pořadí v některých souborech. Při bližším zkoumání jednotlivých souborů uvidíme, že se skutečně v .csv souborech u určitých dataloggerů opakují kusy časové řady. Jde o známou chybu vyskytující se u dat vyčtených pomocí starší verze stahovacího software Lolly. Takto “poškozené” soubory zde demonstrujeme záměrně. V datech žádné měření nechybí, naopak přebývá. Bez strojové kontroly mc_prep_clean() by bylo velmi náročné tento problém odhalit. Podobně je tomu s nesprávným pořadím za sebou jdoucích měření na časové ose.

Pokud v datech chybí jedno měření nebo je v datech několik málo duplicit, většinou jde o chybu, která vznikne při odečtu během sezóny, kdy se přenastaví čas, a řešení těchto chyb lze přenechat na automatické opravě myClim funkcí mc_prep_clean().

Větší počty chyb v datech většinou značí reálnou poruchu dataloggeru či odečítacího software nebo jiný problém a vyžadují hlubší inspekci a individuální manuální řešení na úrovni zdrojových .csv souborů - chybné hodnoty přímo ve zdrojovém souboru manuálně nahradíme za NA nebo využijeme myClim funkce mc_states_insert(), která umožňuje ke konkrétnímu časovému úseku měření daného senzoru přiřadit určitý příznak (tag) a všechny úseky označené určitým tagem posléze pomocí příkazu mc_states_replace() nahradit NA.

V našem případě je tedy potřeba soubory s vysokým počtem chyb otevřít v tabulkovém editoru a prozkoumat. Ukáže se, že nejde o chyby, které by bránily použití dat v analýze, můžeme zde důvěřovat strojové kontrole myClim, která duplicitní hodnoty odstraní a zajistí správnou posloupnost časové řady.

Kromě výpisu mc_info_clean() nám pomůže vizuální inspekce grafického zobrazení dat. Ať už pomocí statických obrázků, či interaktivně v aplikaci myClimGUI (viz dále). Pro celkový přehled o průběhu dat můžeme využít funkce mc_plot_lines() (Obr. 7) a mc_plot_raster() (Obr. 8). Obě funkce jsou optimalizované jak pro zápis výsledných obrázků na disk PC, tak pro prohlížení v R.

V tabulce terénních záznamů vidíme problémy (vykopnutý datalogger) na lokalitě CS_26. Pro příklad si zobrazíme problematickou lokalitu společně s jinou, bezproblémovou lokalitou. Mimo jiné na liniovém grafu (Obr. 7) vidíme, že se časová řada skládá z různých odečtů (souborů), což nám ukáže parametr color_by_logger = TRUE. Anomálii lze vidět i na rastrovém grafu (Obr. 8). Jedním z dalších kroků přípravy dat pro analýzu bude spojení záznamů do jedné časové řady. Dále vidíme prudký pokles vlhkosti na lokalitě CS_26, který napovídá, že se datalogger dostal mimo kontakt s půdou, jak již víme z terénních zápisků.

## Statické zobrazení průběhu teplot a půdní vlhkosti TOMST TMS-4 loggeru 
## Liniový graf s upravením legendy na 3 řádky 
# Použijeme guides() a guide_legend() pro nastavení počtu řádků v legendě
mc_plot_line(micro.data.clean[c("CS_26","CZ2_LUZNICE")], # Vybereme lokality
             sensors = c("TMS_T3","TMS_moist"), # Použijeme konkrétní senzory
             color_by_logger = TRUE) + ggplot2::guides(color 
                                      = ggplot2::guide_legend(nrow = 3))
Liniový graf teplot a vlhkosti na dvou lokalitách.

Liniový graf teplot a vlhkosti na dvou lokalitách. 

## Rastrový graf
mc_plot_raster(micro.data.clean[c("CS_26","CZ2_LUZNICE")],
               sensors = "TMS_moist")
Rastrový graf vlhkosti na dvou lokalitách.

Rastrový graf vlhkosti na dvou lokalitách.

Přesné trvání období, kdy byl TMS-4 datalogger mimo kontakt s půdou, nám pomůže stanovit funkce mc_prep_TMSoffsoil optimalizovaná pro detekci těchto událostí. Funkce vrací vektor TRUE/FALSE, kdy TRUE znamená, že byl datalogger pravděpodobně venku z půdy (Obr. 9 a 10).

## Detekce vykopnutého TMS4 (bez kontaktu s půdou) 
micro.data.prep <- mc_prep_TMSoffsoil(micro.data.clean)

## Vizualizace detailu lokality CS_26 pomocí mc_plot_line 
# a upravíme na legendu na 3 řádky 
mc_plot_line(
  micro.data.prep[c("CS_26", "CZ2_LUZNICE")],  
  sensors = c("TMS_moist", "off_soil"),        
  color_by_logger = TRUE                       
) + ggplot2::guides(color = ggplot2::guide_legend(nrow = 3))
Vizualizace logické proměnné off_soil pro jednu lokalitu.

Vizualizace logické proměnné off_soil pro jednu lokalitu. 

## Kontrola vykopnutí na všech lokalitách 
mc_plot_line(micro.data.prep, sensors = "off_soil")
Vizualizace logické proměnné off_soil pro více lokalit.

Vizualizace logické proměnné off_soil pro více lokalit.

Z grafu, který zobrazuje půdní vlhkost a výsledek detekce nerelevantních měření funkcí mc_prep_TMSoffsoil, můžeme z kombinace barev pro jednotlivé odečty (stažené soubory) a senzoru TRUE/FALSE = vykopnuto/v půdě vyčíst, že došlo k vytažení dataloggeru z půdy v průběhu června 2022. Následně při podzimní návštěvě lokality v témže roce byl datalogger usazen zpět a v zimě 2022/2023 došlo opětovně k jeho vytažení z půdy. Důležitý je pro nás první výskyt hodnoty TRUE na senzoru off_soil a dále informace o tom, kdy došlo k opětovné instalaci. Tyto údaje nám pomohou označit a smazat kompromitované měření z doby, kdy byl datalogger v nesprávné pozici. Skoky mezi TRUE/FALSE nelze interpretovat jako opakované pokusy o instalaci a následné vytažení, jde o artefakt metody detekce vykopnutí, kdy čidlo registrovalo hodnoty půdní vlhkosti a teplot, které jsou v rozmezí přípustných hodnot pro měření v půdě. Problémy vidíme také na lokalitě CZ2_CEPICKA, což nám potvrzují i terénní poznámky o výpadcích vlhkostního senzoru.

Pro vizuální kontrolu dat můžeme využít interaktivní aplikaci myClimGUI, která umožňuje načítání myClim objektů, zoomování a přepínání mezi lokalitami a senzory. Skrze myClimGUI (Obr. 11) je také možné označovat chybná/nechtěná data ke smazání (nahrazení NA).

## Instalace myClimGUI
remotes::install_github("https://github.com/ibot-geoecology/myClimGui")

## spuštění aplikace, načtení dat do myClimGUI
myClimGui::mcg_run(micro.data.prep)
Náhled aplikace myClimGUI
Náhled aplikace myClimGUI

Označit nechtěná měření je také možné pomocí tabulky, kde označíme začátek a konec nechtěné datové řady. Takový postup má výhodu ve větší reprodukovatelnosti než používání interaktivní aplikace myClimGUI, je však uživatelsky méně přívětivý. V interaktivní aplikaci jsme pomocí zoomu na data z lokality CS_26 zjistili, že datalogger TMS-4 byl vytažen z půdy v čase 2022-06-03 18:30 a opětovně instalován v čase 2022-10-19 12:00. Zadáme tedy tyto hodnoty state do tabulky ke všem senzorům na daném loggeru, kde chceme data označit tagem (nahradit NA), a to sice TMS_T1, TMS_T3, TMS_moist. Liberální přístup velí data ze senzoru TMS_T2 ponechat, protože se jedná o přízemní teplotu a lze předpokládat, že vykopnutý datalogger ležel na zemi. Konzervativní přístup velí smazat i časovou řadu TMS_T2, protože může být narušena integrita měření.

## Získání vzorové tabulky "states/tag"
states <- mc_info_states(micro.data.prep)

## Uložení a editace tabulky v MS excel
write.xlsx(states,"./data/statesR.xlsx")

## Načtení editované tabulky s novými tagy "error"
states.insert <- read.xlsx("./data/states_edit.xlsx")

## Přidaní nové hodnoty tagu do myClim objektu
micro.data.prep.s <- mc_states_insert(micro.data.prep,states_table = states.insert)

## Nahrazení hodnoty tagu "error" za NA (smazat)
micro.data.prep.na <- mc_states_replace(micro.data.prep.s, tags="error", 
                                        replace_value = NA)

## Odstranění pomocného senzoru "off_soil"
micro.data.prep.na <- mc_filter(micro.data.prep.na,sensors = "off_soil", reverse = T)
Označení nepoužitelných měření pomocí tagů (z tabulky states insert)
locality_id logger_name sensor_name tag start end value
CS_26 TMS_2 TMS_T1 error 2022-06-03 18:30:00 2022-10-19 12:00:00 offsoil
CS_26 TMS_2 TMS_T3 error 2022-06-03 18:30:00 2022-10-19 12:00:00 offsoil
CS_26 TMS_2 TMS_moist error 2022-06-03 18:30:00 2022-10-19 12:00:00 offsoil
CS_26 TMS_3 TMS_moist error 2022-12-08 18:15:00 2023-04-06 13:00:00 offsoil
CS_26 TMS_3 TMS_T1 error 2022-12-08 18:15:00 2023-04-06 13:00:00 offsoil
CS_26 TMS_3 TMS_T3 error 2022-12-08 18:15:00 2023-04-06 13:00:00 offsoil
CZ2_CEPICKA TMS_1 TMS_T1 error 2022-01-16 23:45:00 2022-01-17 00:45:00 error
CZ2_CEPICKA TMS_1 TMS_T2 error 2022-01-16 23:45:00 2022-01-17 00:45:00 error
CZ2_CEPICKA TMS_1 TMS_T3 error 2022-01-16 23:45:00 2022-01-17 00:45:00 error
CZ2_CEPICKA TMS_2 TMS_moist error 2022-01-16 23:45:00 2022-01-17 00:45:00 error
CZ2_CEPICKA TMS_3 TMS_moist error 2023-01-16 23:45:00 2023-05-05 22:00:00 error

Výsledek úprav můžeme zkontrolovat pomocí grafu (Obr. 12).

Vizuální kontrola průběhu půdní vlhkosti na vyčištěných časových sériích. Série nyní neobsahují zjevně vadná pozorování.

Vizuální kontrola průběhu půdní vlhkosti na vyčištěných časových sériích. Série nyní neobsahují zjevně vadná pozorování.

Kalibrace teplotních dataloggerů

Kalibrace teploměrů dataloggerů před/po umístění na výzkumné plochy je žádoucím krokem. Ačkoliv výrobce uvádí určitou přesnost měření, nežádoucími vnějšími vlivy či výrobní vadou mohou jednotlivé senzory vykazovat větší než výrobcem udávané odchylky od skutečných hodnot. Kalibrace tak pomůže odhalit případné vadné senzory a zlepšit srovnatelnost měření mezi jednotlivými přístroji. Kalibraci můžeme provádět buď laboratorně, nebo, při nižších požadavcích na přesnost měření, i svépomocí. Například umístěním souboru mikroklimatických dataloggerů na místo se stabilní teplotou, kde je ponecháme po dobu minimálně několika hodin pro ustálení teploty (ideálně alespoň 24 hodin). Z naměřených hodnot pak odvodíme odchylku měření jednotlivých senzorů od kalibrovaného měřidla, případně od společného průměru. Pokud je k dispozici pouze jednobodová kalibrace, aplikuje se pouze korekční faktor (cor_factor) vypočtený jako referenční hodnota - naměřená hodnota. Pokud je k dispozici dvou- či vícebodová kalibrace, potom se aplikuje i sklon závislosti odchylky měření (cor_slope). Kalibrovaná hodnota je vypočtena jako:

kalibrovaná hodnota = naměřená hodnota x (cor_slope + 1) + cor_factor

Kalibrační hodnoty pro jednotlivé senzory, identifikované datem provedení kalibrace, výrobním číslem dataloggeru a názvem senzoru, je možné v knihovně myClim ve dvou krocích přiřadit k odpovídajícím časovým řadám a ty následně očistit od známé chyby měření. Pokud je k dispozici jediný kalibrační údaj, aplikuje se na celou časovou řadu (viz příklad). V případě, že jsou k dispozici opakovaně stanovené kalibrační údaje, aplikují se vždy od datumu provedení kalibrace s výjimkou první kalibrace, která se aplikuje rovněž zpětně od prvního měření.

Kalibraci vlhkostního senzoru se věnujeme v kapitole 2.3.

## Načtení tabulky s hodnotami kalibrace (offset)
calib <- read.xlsx("./data/calibrationR.xlsx", detectDates = FALSE)
calib$datetime <- as.POSIXct(convertToDate(calib$datetime))

## Načtení kalibračních dat do myClim
micro.data.prep.c <-
  mc_prep_calib_load(micro.data.prep.na, calib_table = calib)

## Aplikace kalibračního faktoru na časové řady
micro.data.prp.cal <- mc_prep_calib(micro.data.prep.c)
Kalibrační údaje pro jednotlivé senzory
serial_number datetime sensor_id cor_factor
94171044 2021-10-21 TMS_T1 0.0000
94171044 2021-10-21 TMS_T2 0.0625
94171044 2021-10-21 TMS_T3 0.0625
91201145 2021-10-21 Thermo_T 0.0625
21067844 2021-10-21 HOBO_T 0.0060
21061889 2021-10-28 HOBO_T 0.0110
94181410 2021-10-28 TMS_T1 0.5000
94181410 2021-10-28 TMS_T2 0.1875
94181410 2021-10-28 TMS_T3 0.1250
91201144 2021-10-28 Thermo_T 0.0625

Spojování časových řad čistých mikroklimatických dat

Máme připravená vyčištěná a zkalibrovaná data a přejdeme ke spojení datových řad z jednotlivých odečtů, abychom získali spojitou časovou řadu. Spojování dat pomocí knihovny myClim je interaktivní proces, protože v datech často dochází k překryvům, kdy jsou v tentýž okamžik naměřené 2 různé hodnoty, například pokud se na lokalitě měnily dataloggery a z paměti nového dataloggeru byla vyčtena data z doby před zahájením měření na lokalitě. Vzácně se také objevují jednotlivá duplikovaná měření vzniklá v důsledku záporné korekce času vnitřních hodin dataloggeru při odečtu a následným opětovným měřením v měřícím intervalu. Rovněž se může jednat o případy opakovaného vyčtení celé paměti dataloggerů, kdy jsou tak některé úseky vyčteny vícekrát. V takovém případě myClim automaticky zachová pouze jeden údaj k času měření. V případech, kdy se časové řady překrývají, ale hodnoty měření nejsou shodné, myClim vypíše údaje o tom, na které lokalitě, ve kterém dataloggeru a senzoru k překryvu dochází a jaký je časový rozsah překryvu a zobrazí interaktivní graf přiblížený na místo překryvu, kde si může uživatel danou situaci prohlédnout. Následně myClim čeká na akci uživatele, který musí zadat do konzole jednu z voleb (číslici 1 - 6) podle toho, zda má myClim upřednostnit hodnoty ze starší časové řady (dataloggeru) (1), z mladší časové řady (2), přeskočit tento konkrétní překryv a zachovat hodnoty z obou dataloggerů (3), u tohoto a všech dalších překryvů v celém objektu upřednostnit vždy starší datalogger (4), vždy mladší datalogger (5), nebo proces spojování ukončit (6).

# Spojení časových řad
micro.data.join <- mc_join(micro.data.prp.cal)
## vizualizace detailu lokality CS_26 
## Vzhledem k interaktivní povaze zde není výstup
mc_plot_line(micro.data.join[c("CS_26","CZ2_LUZNICE")],
             sensors = c("TMS_moist","TMS_T1"),
             color_by_logger = TRUE)

Výstup spojování zkontrolujeme vizuálně v myClimGUI, pro příklad je zobrazen znovu graf lokalit ‘CS_26’ a ‘CZ2_LUZNICE’ s parametrem color_by_logger = TRUE podobně jako výše. Na grafu je patrné, že jsme odstranili problematická měření a spojili časovou řadu do jednoho celku.

Doplnění chybějících měření

Posledním krokem v přípravě mikroklimatických dat pro analýzu je nahrazení chybějících hodnot. Pomocí myClim lze bezpečně doplnit zejména krátké výpadky v řádu 1 - 5 chybějících pozorování v závislosti na časovém kroku daného dataloggeru, např. když přesně nenavazují data pro jednotlivé odečty. Takové hodnoty je možné doplnit jednoduchou lineární interpolací dvou sousedních měření implementovanou v myClim. Doplnění delších úseků pomocí myClim nedoporučujeme provádět.

## Doplnění chybějících hodnot v maximální délce 10 za sebou jdoucích měření
micro.data.fill <- mc_prep_fillNA(micro.data.join,maxgap = 10)

Uložení dat připravených pro analýzu

Jestliže jsme při načítání mikroklimatických dat ze souborů nevyužili možnost nahrát doplňující informace o lokalitách, je možné je doplnit dodatečně. V tomto příkladě doplníme informace o geografických souřadnicích lokalit. Ty se budou při dalším zpracování hodit, například pro výpočet posunu místního času vůči koordinovanému světovému času (UTC) pro synchronizaci fotoperiody napříč lokalitami. Na závěr uložíme objekt obsahující veškeré časové řady měření a metadata o lokalitách do souboru .rds. Objekt myClim uložený v souboru .rds je datově úsporný (v tomto případě zabírá soubor ‘prep_CZ2_10_join.rds’ pouze 11 MB ve srovnání s 129 MB zdrojových .csv souborů). Pro uložení myClim objektu doporučujeme vždy používat funkci myClim::mc_save(), která zajistí zpětnou kompatibilitu objektu s novými verzemi knihovny myClim, namísto základní funkce saveRDS().

## Přidat zeměpisnou délku v souřadnicovém systému WGS84 (lon_wgs84) 
# v desetinných stupních
micro.data.fill <- mc_prep_meta_locality(micro.data.fill,
                                         list(CZ2_KRKLAB = 15.5640182, 
                                              CZ2_KLINOV = 12.9735495, 
                                              CZ2_BARINY = 17.9556943, 
                                              CZ2_HRUSUT = 17.4393209, 
                                              CZ2_MASTALE = 16.1369899, 
                                              CZ2_HODDUB = 17.0748401, 
                                              CZ2_CEPICKA = 16.3635405, 
                                              CZ2_LUZNICE = 14.8428597, 
                                              CS_26 = 13.9910361, 
                                              CZ2_BUKOVEC = 15.3606035),
                                         param_name = "lon_wgs84")
## Přidat zeměpisnou šířku v souřadnicovém systému WGS84 (lat_wgs84) 
# v desetinných stupních
micro.data.fill <- mc_prep_meta_locality(micro.data.fill,
                                         list(CZ2_KRKLAB = 50.7512076, 
                                              CZ2_KLINOV = 50.3935486, 
                                              CZ2_BARINY = 49.6310008,
                                              CZ2_HRUSUT = 49.6932542,
                                              CZ2_MASTALE = 49.8115909, 
                                              CZ2_HODDUB = 48.8803506, 
                                              CZ2_CEPICKA = 49.4982277,
                                              CZ2_LUZNICE = 48.9940496, 
                                              CS_26 = 50.5953875,
                                              CZ2_BUKOVEC = 50.8143946),
                                         param_name = "lat_wgs84")
## Uložení finálního myClim objektu
mc_save(micro.data.fill, file = "./data/prep_CZ2_10_join.rds")

2. Mikroklimatické proměnné

V předchozích kapitolách jsme se věnovali tomu, jak vytvořit co nejrelevantnější mikroklimatická měření z dataloggerů, jak je z dataloggerů stáhnout, vyčistit, validovat, kalibrovat a případně napojit na sebe časové řady z jedné lokality s více odečty.

V této sekci si ukážeme, jak z časových řad měření mikroklimatu odvodit biologicky relevantní mikroklimatické proměnné. Ty charakterizují mikroklima na lokalitách za dané časové období (např. roční/sezónní/měsíční/denní průměry, minima, maxima či variabilita).

2.1. Čas / období

Mikroklimatické dataloggery zpravidla měří v poměrně jemném časovém rozlišení (desítky vteřin až jednotky hodin). Pro biologické analýzy je však většinou žádoucí data v čase agregovat pomocí matematických funkcí a, tím získat popisné charakteristiky pro dané období (např. denní, měsíční, sezónní, roční…). Volba délky periody pro výpočet mikroklimatických proměnných závisí na otázce, na niž chceme s pomocí mikroklimatických dat získat odpověď. Pro některé systémy a proměnné je vhodné agregovat časové řady přes kalendářní rok. Tam, kde do úvahy vstupuje například půdní vlhkost ovlivněná táním sněhu, je vhodnější volba hydrologický rok (pro ČR období 1. listopadu - 31. října). Pokud není pro zkoumaný systém některá část roku relevantní, pak může být žádoucí používat mikroklimatická data ze specifického období, např. vegetační sezóna. Volba vhodného časového úseku pro agregaci dat hraje významnou roli, protože ovlivňuje, zda zachytíme v datech biologicky relevantní signál, nebo ho naopak překryje šum z měření v období, které není pro danou otázku relevantní.

Pro zpracování mikroklimatických časových řad je také důležité časové pásmo, ve kterém loggery měří, a to zejména pokud analyzujeme data z většího území např. globální. Časový posun může hrát roli zejména při agregacích dat a výpočtech biologicky relevantních proměnných. Mnohé loggery měří v UTC (TOMST TMS-4) bez ohledu na přání uživatele či časovou zónu, kde jsou instalovány, jiné loggery je možné nastavit na konkrétní časové pásmo. Zde používaná knihovna myClim striktně předpokládá, že časová řada jakýchkoli vstupních dat je v UTC. Pokud uživatel používá loggery, které neměří v UTC, je žádoucí při jejich načítání do myClim zadat parametr tz_offset, tedy počet minut o  kolik je daná časová řada posunutá vůči UTC. myClim pak tento offset při následných agregacích zohlední.

Pokud analyzujeme globální dataset s loggery měřícími v UTC a agregujeme data například na denní průměry, je dobré synchronizovat na všech lokalitách fotoperiodu, tedy dodat informaci o tom, kdy na dané lokalitě nastává poledne. To můžeme udělat buď ručním nastavením parametru tz_offset nebo vypočítat ze zeměpisné délky, kterou musíme mít zadanou v metadatech lokality. Posun vůči UTC pro synchronizaci fotoperiody provedeme pomocí funkce mc_prep_solar_tz(). Tato funkce neovlivní data ani časovou řadu samotnou, pouze přidá do metadat lokality informaci o  posunu. Tento posun se pak aplikuje při agregacích na dny a delší periody. Na kratší periody než dny se pro agregaci posun neaplikuje.

Pro příklady používáme data z deseti lokalit. Na všech lokalitách jsou TOMST TMS-4 loggery měřící vzdušnou teplotu v 15 cm, 2 cm nad povrchem půdy, v 8 cm pod povrchem půdy, vlhkost půdy ve svrchní vrstvě 0 - 15 cm, na devíti lokalitách Thermologger měřící vzdušnou teplotu ve 200 cm nad zemí, a HOBO U23 logger měřící vzdušnou teplotu a vlhkost ve 150 cm nad zemí. Časová řada přibližně od dubna 2021 do října 2023, ale loggery neměřily stejnou dobu, v datech jsou chybějící úseky.

## Načtení připraveného myClim objektu z kapitoly 1 - Příprava dat ze souborů
# Instalace a načtení knihoven
suppressWarnings(lapply(c("myClim", "ggplot2"), function(pkg) {
  if (!require(pkg, character.only = TRUE))
    install.packages(pkg, repos = "https://mirrors.nic.cz/R/")
  library(pkg, character.only = TRUE)
}))

micro.data <- mc_load("./data/prep_CZ2_10_join.rds")
micro.data
## vykreslení grafu průběhu teplot a relativní vlhkosti
p <- mc_plot_line(micro.data[5:7], sensors = c("TMS_T3", "HOBO_RH"))

## manuální uprava barev
p <-
  p + ggplot2::scale_color_manual(values = scales::alpha(c("darkblue", "hotpink"), 0.7),
                                  name = NULL)
print(p)
Vykreslení grafu průběhu teplot a relativní vlhkosti.

Vykreslení grafu průběhu teplot a relativní vlhkosti. 

## výpočet offsetu vůči UTC, synchronizace fotoperiody
micro.data<-mc_prep_solar_tz(micro.data)
Výpis informací o objektu micro.data
locality_id lon_wgs84 lat_wgs84 elevation tz_offset
CS_26 13.99104 50.59539 NA 56
CZ2_BARINY 17.95569 49.63100 NA 72
CZ2_BUKOVEC 15.36060 50.81439 NA 61
CZ2_CEPICKA 16.36354 49.49823 NA 65
CZ2_HODDUB 17.07484 48.88035 NA 68
CZ2_HRUSUT 17.43932 49.69325 NA 70
CZ2_KLINOV 12.97355 50.39355 NA 52
CZ2_KRKLAB 15.56402 50.75121 NA 62
CZ2_LUZNICE 14.84286 48.99405 NA 59
CZ2_MASTALE 16.13699 49.81159 NA 65

2.2. Teplota

Teplota a z ní přímo odvozené proměnné jsou jedním z nejpoužívanějších (mikro)klimatických parametrů. Je to jednak kvůli skutečné biologické relevanci, ale také kvůli nízké technické náročnosti měření teploty. Mezi základní teplotní proměnné bezesporu řadíme průměrnou teplotu za dané období a teplotní extrémy. Průměr je metodicky snadno uchopitelná proměnná, nicméně i zde se setkáme s různými způsoby výpočtu. Aritmetický průměr ze všech pravidelných pozorování v průběhu 24 hodin označuje meteorologický slovník (https://slovnik.cmes.cz/) jako “průměrná denní teplota pravá” a tento způsob výpočtu je v mikrometeorologii preferovaný. V meteorologii na území ČR se z historických důvodů používá průměrná denní teplota vzduchu jako vážený průměr teplot v termínech měření 7h, 14h a 21h místního času, přičemž měření v 21h mělo dvojnásobnou váhu. Americká NOAA zase definuje průměrnou denní teplotu jako průměr z minimální a maximální denní teploty, což může vést k jistým odchylkám takto stanovených teplot a komplikovat srovnatelnost uváděných hodnot z různých zdrojů.

Rovněž způsob stanovení minimálních a maximálních teplot se v praxi liší. Ve zprávách SYNOP je uváděno minimum z teplot za noční období (od 18 do 06 UTC), pro klimatologické účely je referována minimální teplota za období 24 hodin před večerním termínem (21h). Pro odstranění vlivu odlehlých pozorování (v důsledku nahodilých klimatických událostí či ovlivněných měření) se pro charakteristiku teplotních extrémů v mikroklimatologii používají percentily (např. 5 percentil pro minima a 95 percentil pro maxima) z pozorovaných hodnot.

Důležitou odvozenou proměnnou jsou pak sumy teplot nad či pod určitou prahovou teplotou (suma efektivních teplot, Growing degree days, GDD / Freezing degree days, FDD). Ty určují například délku sezóny pro vývoj plodin či škůdců, nebo naopak délku zamrznutí relevantní pro přežití některých organismů (škůdců / parazitů). Volba prahové teploty závisí na výzkumné otázce, v agronomii jsou stanovené rozdílné prahové hodnoty pro jednotlivé plodiny podle jejich fyziologických nároků. V ekologii se pro výpočet zpravidla používá prahové hodnoty 5°C, která je všeobecně přijímána jako limitní teplota pro růst rostlin.

Výpočet virtuálních senzorů z teploty

Pro snazší práci s daty je vhodné některé proměnné, odvozené z jednoho nebo více senzorů, uložit přímo jako tzv. virtuální senzor daného loggeru. Tyto senzory mají stejnou periodu měření jako ty senzory, ze kterých byly odvozeny. Například z teplot dopočteme virtuální senzor Growing degree days (GDD) nebo Freezing degree days (FDD), tedy sumu efektivních/mrazových teplot nad/pod zvolenou prahovou hodnotou - v senzoru jsou uloženy hodnoty, kterými každá perioda přispívá dané proměnné. Teprve po agregaci na dny nebo delší časový úsek můžeme hovořit přímo o FDD nebo GDD.

# Instalace a načtení knihoven
suppressWarnings(lapply(c("myClim", "dplyr", "tidyr"), function(pkg) {
  if (!require(pkg, character.only = TRUE))
    install.packages(pkg, repos = "https://mirrors.nic.cz/R/")
  library(pkg, character.only = TRUE)
}))

# Pro příklad použijeme TMS_T3 senzor 15 cm nad zemí, z TMS-4 TOMST loggeru 
micro.data <- mc_calc_gdd(micro.data, sensor = "TMS_T3", t_base = 5)
micro.data <- mc_calc_fdd(micro.data, sensor = "TMS_T3", t_base = 5)

Výpočet teplotních proměnných

Pro příklad definujme vegetační sezónu jako období 1.5. - 31.8. Pomocí parametru period="custom" můžeme zavést libovolné období, pro které se budou proměnné počítat. Pokud časová řada obsahuje měření z více let, proměnné se vypočítají pro všechny roky, kde je dostatečné pokrytí daty definované parametrem min_coverege=0.9,tedy pokud máme pro dané období více než 90% dat, proměnná je vypočtena, pokud méně než 90% myClim vrátí NA, tedy prázdnou hodnotu. Agregační funkce jsou v příkladu níže definované jako list (seznam) funkcí, kdy specificky vybíráme, pro které senzory se aplikuje jaká funkce. Senzory, pro které není definovaná funkce, do výpočtu nevstupují a ve výsledném objektu nejsou zahrnuty. Parametrem percentiles = c(5,95) definujeme, které percentily se mají vypočítat ve funkci “percentile”.

## výpis senzorů, které jsou v myClim objektu k dispozici pro výpočty
levels(factor(mc_info(micro.data)$sensor_name))
## [1] "FDD5"      "GDD5"      "HOBO_RH"   "HOBO_T"    "Thermo_T"  "TMS_moist"
## [7] "TMS_T1"    "TMS_T2"    "TMS_T3"
## vybereme teplotní senzory a požadované agregační funkce
micro.veget <- mc_agg(micro.data,
                    period = "custom",
                    custom_start = "05-01",
                    custom_end = "08-31",
                    percentiles = c(5,95),
                    min_coverage = 0.9,
                    fun=list(
                      FDD0="sum",
                      GDD5="sum",
                      HOBO_T=c("mean","percentile","range"),
                      Thermo_T=c("mean","percentile","range"),
                      TMS_T1=c("mean","percentile","range"),
                      TMS_T2=c("mean","percentile","range"),
                      TMS_T3=c("mean","percentile","range")))

Výsledný myClim objekt obsahuje všechny vstupní senzory s příponou dané funkce. Každý senzor má počet měření podle toho, kolik úseků časové řady splnilo podmínku “min_coverage”. V našem příkladu jsou to pro některé lokality 3 hodnoty (3 vegetační sezony), v některých případech 2 hodnoty, někdy pouze jedna, protože vstupní časová řada byla příliš krátká.

## výpis výsledných senzorů
levels(factor(mc_info(micro.veget)$sensor_name))
##  [1] "GDD5_sum"              "HOBO_T_mean"           "HOBO_T_percentile5"   
##  [4] "HOBO_T_percentile95"   "HOBO_T_range"          "Thermo_T_mean"        
##  [7] "Thermo_T_percentile5"  "Thermo_T_percentile95" "Thermo_T_range"       
## [10] "TMS_T1_mean"           "TMS_T1_percentile5"    "TMS_T1_percentile95"  
## [13] "TMS_T1_range"          "TMS_T2_mean"           "TMS_T2_percentile5"   
## [16] "TMS_T2_percentile95"   "TMS_T2_range"          "TMS_T3_mean"          
## [19] "TMS_T3_percentile5"    "TMS_T3_percentile95"   "TMS_T3_range"
## převod z myClim objektu na prostou tabulku, dlouhý formát
df.veget <- mc_reshape_long(micro.veget)

## vizualizace
## průměrná, minimální a maximální teplota vzduchu 200 cm nad zemí (Thermo_T)
df.veget$datetime <- as.character(df.veget$datetime)
viz <- df.veget %>%
  filter (
    sensor_name %in% c(
      "Thermo_T_mean",
      "Thermo_T_percentile5",
      "Thermo_T_percentile95"
    )
  ) %>%
  filter (!is.na(value))

viz <- pivot_wider(viz, 
                   names_from = "sensor_name",
                   values_from = "value")

p <-
  ggplot(viz, aes(x = locality_id, y = Thermo_T_mean, col = datetime)) +
  geom_point(position = position_dodge(width = 0.5)) +
  geom_errorbar(
    aes(ymin = Thermo_T_percentile5,
        ymax = Thermo_T_percentile95,
        width = 0.15),
    position = position_dodge(width = 0.5)
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(
    angle = 60,
    vjust = 0.7,
    hjust = 0.5
  )) +
  theme(legend.position = "bottom") +
  labs(title = "Teplota ve 2 m", 
       subtitle = "průměr, minimum (P5) a maximum (P95)")

print(p)
Průměrná, minimální a maximální teplota ve 2 m na několika lokalitách.

Průměrná, minimální a maximální teplota ve 2 m na několika lokalitách.

Výpočet myClim teplotních proměnných

Pro výpočet statistik agregujících měření za delší časové úseky slouží v knihovně myClim funkce mc_agg(). Ta, která aplikuje funkce (předdefinované, nebo uživatelem definované) na série měření v zadaném časovém období.

Pro výpočet standardizované sady teplotních ukazatelů slouží v knihovně myClim funkce mc_env_temp(), která pro zadané časové období vypočítá následující souhrnné ukazatele:

  • min5p: Minimální teplota (5. percentil denních minimálních teplot)
  • mean: Průměrná teplota (průměr denních teplot)
  • max95p: Maximální teplota (95. percentil denních maxim)
  • drange: průměrná denní amplituda teplot (průměrný rozdíl denních maxim a minim)
  • GDD5: suma efektivních teplot (suma teplot nad prahovou hranicí, výchozí hodnota 5°C)
  • FDD0: suma mrazových teplot (suma záporných teplotních odchylek od prahové hodnoty, výchozí hodnota 0°C)
  • frostdays: počet mrazových dnů (dny s minimální teplotou <0°C)

Při ručním výpočtu pomocí funkce mc_agg() se přímo agregují hodnoty dle zadaných parametrů uživatelem. Na rozdíl od toho, při použití funkce mc_env_temp() se teplotní proměnné agregují nejprve na denní krok a odtud pak na periodu zadanou uživatelem. Tento postup je bližší meteorologickým standardům. Dalším rozdílem je, že funkce mc_env_temp() vrací přímo plochou tabulku. S tou nejde již dále pracovat pomocí knihovny myClim, ale je přímo připravena pro následné analýzy.

Funkce mc_env_temp() vyžaduje, aby všechny senzory na lokalitách měly shodný časový krok. V našem souboru dat, který používáme jako příklad, tomu tak není. Na lokalitách jsou 2 typy loggerů: TOMST TMS měří v intervalu 15 minut a HOBO v intervalu 30 minut. Nejjednodušší tedy je sjednotit časový krok na 30 minut pomocí agregace s použitím např. průměru mc_agg(period="30 min", fun="mean").

Pro příklad zde vypočteme proměnné pro kalendářní rok. I zde definujeme minimální pokrytí daty, respektive maximální podíl povolených chybějících hodnot parametrem min_coverage. Ve výsledné tabulce je velké množství řádků s chybějící hodnotou, což indikuje to, že pro dané období nebyla splněna podmínka min_coverage. Můžeme zvážit jiné nastavení prahové hodnoty, chybějící (NA) hodnoty můžeme následně z výsledné tabulky odstranit.

## sjednocení periody TMS a HOBO loggerů na 30 min
micro.data30 <- mc_agg(micro.data, period = "30 min", fun = "mean")

## výpočet myClim teplotních proměnných pro vegetační sezónu
micro.temp <-
  mc_env_temp(
    micro.data30,
    period = "year",
    min_coverage = 0.9,
    gdd_t_base = 5,
    fdd_t_base = 0
  )

## odstranění chybějících hodnot
micro.temp <- filter(micro.temp, !is.na(value))

## vizualizace
micro.temp$fun <- micro.temp$sensor_name %>%
  strsplit(".", fixed = T) %>%
  lapply("[", 3) %>%
  unlist()

viz <- micro.temp %>%
  filter (fun %in% c("mean", "min5p", "max95p")) %>%
  filter (!is.na(value))
viz$sensor_name <- NULL
viz <- pivot_wider(viz, names_from = "fun",
                   values_from = "value")

p <- ggplot(viz, aes(x = locality_id, y = mean, col = height)) +
  geom_point(position = position_dodge(width = 0.5)) +
  geom_errorbar(aes(ymin = min5p,
                    ymax = max95p,
                    width = 0.15),
                position = position_dodge(width = 0.5)) +
  theme_minimal() +
  theme(axis.text.x = element_text(
    angle = 60,
    vjust = 0.7,
    hjust = 0.5
  )) +
  theme(legend.position = "bottom") +
  ylab("Teplota [°C]") +
  xlab(NULL) +
  labs(title = "Teplota v různých výškách /hloubkách",
       subtitle = "průměr, minimum (P5) a maximum (P95) v roce 2022")

print(p)
Průměrná, minimální a maximální teplota v různých výškách na několika lokalitách.

Průměrná, minimální a maximální teplota v různých výškách na několika lokalitách.

2.3. Vlhkost půdy

Numerickým výstupem senzorů půdní vlhkosti u mikroklimatických loggerů často nejsou přímo hodnoty objemové půdní vlhkosti či jiné standardní proměnné. Surové hodnoty vlhkostního signálu z loggerů je možné přímo použít pro relativní porovnání či pro získání hrubé informace o vlhkostních poměrech v daném systému. Vhodnější však je pokusit se surová měření převést na standardní veličiny, pokud k tomu výrobce loggeru poskytuje potřebné informace a kalibrační protokoly.

Výpočet objemové půdní vlhkosti

Senzory půdní vlhkosti (TOMST TMS4) poskytují po zpracování surového signálu informaci o objemové vlhkosti půdy. Tato veličina dává základní představu o dostupnosti vláhy pro rostliny a další procesy závislé na obsahu vody v půdě. Nicméně hydrolimity pro jednotlivé typy půd se liší, a proto je třeba naměřené hodnoty interpretovat ve vztahu k vlastnostem půdy na konkrétním stanovišti. Například bod vadnutí se pohybuje mezi 5% objemové vlhkosti u písčitých půd po 20% u jílovitých půd.

Převod TMS4 surového vlhkostního signálu na objemovou vlhkost je stále předmětem diskusí a optimalizací. Existuje značná škála možností jaký použít půdní typ ať už výběrem z před-definovaných či zadáním kalibračních parametrů vlastního půdního typu. Dále hraje roli kalibrace jednotlivých loggerů (hodnoty, kterých senzor dosahuje na vzduchu a po zanoření do vody). Pro více informací doporučujeme následující zdroje:

Pro zjednodušení zde použijeme univerzální kalibrační křivku odvozenou v práci (Kopecký et al., 2021). Kromě převodu na objemovou vlhkost funkce mc_calc_vwc() v základním nastavení také smaže hodnoty vlhkosti v období, kdy byla půda zamrzlá frozen2NA = TRUE, protože taková měření jsou biologicky málo relevantní.

## výpočet virtuálního senzoru objemové vlhkosti
## převod TMS vlhkostního signálu na objemovou vlhkost
moist.data <- mc_calc_vwc(micro.data, soiltype = "universal",
                          frozen2NA = TRUE)

## vizualizace
p <- mc_plot_line(moist.data,sensors = c("TMS_moist", "VWC_moisture"))
p <- p + ggplot2::scale_color_manual(values = c("darkblue","steelblue"), 
                                   name = NULL)
print(p)
Vizualizace průběhu surových signálů vlhkostí a volumetrických vlhkostí pro několik lokalit.

Vizualizace průběhu surových signálů vlhkostí a volumetrických vlhkostí pro několik lokalit.

Proměnné vlhkosti půdy

Pro výpočet standardizovaných ukazatelů půdní vlhkosti slouží v myClim knihovně funkce mc_env_moist(), která pro časové série vlhkostních měření vypočítá následující parametry pro zadané časové intervaly:

  • VWC.5p: minimální objemová vlhkost půdy (5. percentil objemové vlhkosti půdy)
  • VWC.mean: průměrná objemová vlhkost půdy
  • VWC.95p: maximální objemová vlhkost půdy (95. percentil objemové vlhkosti půdy)
  • VWC.sd: vyrovnanost vodního režimu (směrodatná odchylka objemové vlhkosti půdy)

Podobně jako funkce pro výpočet myClim teplotních proměnných i zde se vypočtené proměnné vrací v podobě ploché tabulky, kterou již není možné používat v knihovně myClim, ale je připravena pro následné analýzy. Před použitím funkce mc_env_moist() je nutné nejprve spočítat virtuální senzory objemové vlhkosti. Funkce nepracuje s TMS jednotkami ani s jinými přímými měřeními, které nejsou objemovou vlhkostí.

I zde funkce vyžaduje, aby všechny senzory na lokalitách měly shodný časový krok. V našem souboru dat, který používáme jako příklad, tomu tak není. Na lokalitách jsou 2 typy loggerů: TOMST TMS měří v intervalu 15 minut a HOBO v intervalu 30 minut. Můžeme podobně jako v příkladu teploty výše senzory převést na shodný krok. V tomto případě však z myClim objektu HOBO loggery odstraníme, protože je pro tento výpočet nepotřebujeme. I zde nám funkce vrátí velké množství chybějících hodnot kvůli parametru min_coverage.

## využijeme myClim objekt s virtuálním senzorem objemové vlhkosti
## odstraníme z myClim objektu HOBO loggery, mají jiný časový krok
moist.data.TMS <-
  mc_filter(moist.data, logger_types = "HOBO_U23-001A", reverse = T)

## výpočet proměnných pouze pro TOMST TMS loggery
micro.moist <-
  mc_env_moist(moist.data.TMS, period = "year", min_coverage = 0.9)

## odstranění chybějících hodnot
micro.moist <- filter(micro.moist, !is.na(value))

## vizualizace
micro.moist$fun <- micro.moist$sensor_name %>%
  strsplit(".", fixed = T) %>%
  lapply("[", 3) %>%
  unlist()
micro.moist$sensor_name <- NULL
micro.moist.w <- pivot_wider(micro.moist, names_from = "fun",
                             values_from = "value") %>%
  filter(datetime == as.Date("2022-01-01"))

ggplot(micro.moist.w, aes(x = locality_id, y = mean)) +
  geom_point() +
  geom_errorbar(aes(ymin = `5p`,
                    ymax = `95p`,
                    width = 0.15)) +
  theme_minimal() +
  theme(axis.text.x = element_text(
    angle = 60,
    vjust = 0.7,
    hjust = 0.5
  )) +
  theme(legend.position = "bottom") +
  ylab("Objemová vlhkost půdy [%]") +
  xlab(NULL) +
  labs(title = "Objemová vlhkost půdy",
       subtitle = "průměr, minimum (P5) a maximum (P95) pro rok 2022")
Objemová vlhkost půdy na několika lokalitách.

Objemová vlhkost půdy na několika lokalitách.

2.4. Vlhkost vzduchu

Z hlediska biologické relevantnosti je důležitým ukazatelem vzdušné vlhkosti sytostní doplněk (vapor pressure deficit, VPD). Ten vyjadřuje rozdíl mezi maximálním tlakem vodní páry při dané teplotě a skutečným tlakem vodních par a odráží tak, lépe než relativní vlhkost vzduchu, např. ztrátu vody transpirací a riziko vodního stresu rostlin či požární riziko. Při znalosti teploty a relativní vlhkosti vzduchu (např. z dat čidel HOBO U23 Pro v2) lze sytostní doplněk vypočítat funkcí mc_calc_vpd().

Funkce také vyžaduje, aby všechny senzory na lokalitách měly shodný časový krok. Opět tedy z myClim objektu odstraníme loggery, které nebudeme pro tento výpočet potřebovat, ponecháme pouze HOBO loggery. I zde nám funkce vrátí velké množství chybějících hodnot kvůli parametru min_coverage.

## ponecháme pouze HOBO loggery
micro.data.HOBO <- mc_filter(micro.data,logger_types = "HOBO_U23-001A",reverse = F)

## výpočet virtuálního senzoru sytostního doplňku (VPD) z loggerů HOBO
vpd.data <- mc_calc_vpd(micro.data.HOBO)

## vizualizace
p <- mc_plot_line(vpd.data,sensors = c("HOBO_RH", "VPD"),scale_coeff = 20)
p <- p + ggplot2::scale_color_manual(values = c("red","black"), 
                                   name = NULL)
print(p)
Průběhy relativní vlhkosti vzduchu a sytostní doplněk pro několik lokalit.

Průběhy relativní vlhkosti vzduchu a sytostní doplněk pro několik lokalit.

myClim proměnné vlhkosti vzduchu

Pro výpočet standardizovaných proměnných vzdušné vlhkosti (sytostního doplňku VPD) slouží v myClim knihovně funkce mc_env_vpd(), která pro časové série virtuálního senzoru VPD vypočítá následující parametry:

  • VPD.mean: průměrný sytostní doplněk (průměr z denních průměrů VPD)
  • VPD.max95p: maximální sytostní doplněk (95. percentil denních maxim VPD)
## použijeme myClim objekt s virtuálním senzorem VPD vypočteným výše
micro.vpd <- mc_env_vpd(vpd.data, period = "year", min_coverage = 0.9)

## odstranění chybějících hodnot
micro.vpd <- filter(micro.vpd, !is.na(value))

## vizualizace
micro.vpd$proměnná <- micro.vpd$sensor_name %>%
  strsplit(".", fixed = T) %>%
  lapply("[", 3) %>%
  unlist()

ggplot(micro.vpd, aes(x = locality_id, y = value, col = proměnná)) +
  geom_point() +
  theme_minimal() +
  theme(legend.position = "bottom") +
  ylab("Sytostní doplněk [kPa]") +
  xlab(NULL) +
  labs(title = "Sytostní doplněk",
       subtitle = "průměr a maximum (P95) pro rok 2022")
Průměrný a maximální (95p) sytostní doplněk na několika lokalitách.

Průměrný a maximální (95p) sytostní doplněk na několika lokalitách.

2.5. Odhad přítomnosti sněhové pokrývky

Sněhová pokrývka zásadním způsobem ovlivňuje teploty při povrchu půdy. Izolační vlastnosti sněhové pokrývky a skupenské teplo tání způsobují, že teplotní výkyvy pod sněhovou pokrývkou jsou výrazně utlumené a v případě částečně mokrého sněhu (což je v podmínkách ČR převládající stav) se teplota na rozhraní půda/sníh drží v blízkosti 0°C. Těchto principů lze využít pro detekci sněhové pokrývky z mikroklimatických měření. K tomu slouží funkce mc_calc_snow(), která detekuje přítomnost sněhu u těch měření, která splňují podmínky pro maximální teplotu (výchozí hodnota 1,25°C) a maximální rozsah teplot (výchozí hodnota 1°C) v plovoucím časovém okně (výchozí hodnota 3 dny) na zvoleném teplotním senzoru (výchozí volba přízemní senzor čidla Tomst TMS4). Spolehlivost detekce sněhu s těmito výchozími hodnotami se pohybuje okolo 95% (srovnáním s každodenními fotkami lokality fotopastí). Pro odvození základních statistik (doba trvání sněhové pokrývky, první/poslední den sněhové pokrývky a první/poslední den souvislé sněhové pokrývky) slouží funkce mc_calc_snow_agg(). Souvislou sněhovou pokrývku lze uživatelsky definovat (výchozí hodnota 3 dny).

## Detekce sněhu na základě přízemních teplot (TMS_T2)
## Výstup uložen jako virtuální senzor "snih"
micro.snow <- mc_calc_snow(micro.data, 
                           sensor = "TMS_T2", 
                           output_sensor = "snih", 
                           range = 1,
                           tmax = 1.25, 
                           days = 3)
## Vizualizace
mc_plot_line(micro.snow[c("CZ2_BUKOVEC",
                           "CZ2_KRKLAB",
                           "CS_26")], sensors = c("TMS_T2","snih"))
Vizualizace detekovaného sněhu na časových řadách několika stanic.

Vizualizace detekovaného sněhu na časových řadách několika stanic. 

## Výpočet základní statistiky sněhové pokrývky pro zimní sezónu 2021/2022
## Výběr dat pro zvolené období
micro.snow.2021 <- mc_prep_crop(micro.snow, 
                                start = as.POSIXct("2021-10-01", tz="UTC"), 
                                end = as.POSIXct("2022-06-01", tz="UTC"))

## Výpočet doby trvání sněhové pokrývky, 
## souvislá sněhová pokrývka definována jako minimálně 7 po sobě jdoucích dní
tabulka.snih.2021 <- mc_calc_snow_agg(micro.snow.2021, 
                                      snow_sensor = "snih", 
                                      period = 7)
Trvání sněhové pokrývky automaticky detekované z průběhu teplot
locality_id snow_days first_day last_day first_day_period last_day_period
CS_26 0 NA NA NA NA
CZ2_BARINY 0 NA NA NA NA
CZ2_BUKOVEC 153 2021-11-24 2022-04-25 2021-11-24 2022-04-25
CZ2_CEPICKA 40 2021-12-04 2022-02-06 2021-12-04 2022-02-06
CZ2_HODDUB 0 NA NA NA NA
CZ2_HRUSUT 5 2022-01-25 2022-01-29 NA NA
CZ2_KLINOV 119 2021-12-04 2022-05-01 2021-12-04 2022-05-01
CZ2_KRKLAB 165 2021-11-30 2022-05-13 2021-11-30 2022-05-13
CZ2_LUZNICE 0 NA NA NA NA
CZ2_MASTALE 0 NA NA NA NA

3. Vyhodnocení a srovnání

Mikroklima v terénu měříme většinou s cílem kvantifikovat mikroklimatické rozdíly mezi jednotlivými stanovišti nebo pro kvantifikaci rozdílu mezi mikroklimatem zájmového stanoviště a referenční hodnotou pro makroklima, například ze staničních měření sítě ČHMÚ. Následující příklady názorně ukazují, jak s daty pracovat a jak ověřit, zda se mikroklima liší.

3.1 Vliv druhu stromu na mikroklima podrostu

Úvod

Často využívaným postupem je hledání signifikantního rozdílu mikroklimatu mezi kategoriemi. Mohou to být například různé druhy stromů, typy stanovišť, různé disturbance aj. Typicky se používá alespoň 6–10 opakování na kategorii, pokud očekáváte střední efekt a máte relativně nízkou variabilitu. Pokud je variabilita vyšší nebo očekáváte menší efekt, může být zapotřebí více nezávislých opakování (10–20). Kromě rozdílů mezi kategoriemi je užitečný i rozdíl vůči kontrole. V případě disturbancí je dobrým typem kontrolní kategorie vzrostlý les. Pokud nás zajímá vliv jednotlivých dřevin, může kontrolní kategorie být na bezlesí.

V ideálním případě by jednotlivé mikroklimatické stanice měly být na sobě nezávislé, například by neměly být příliš blízko u sebe. Pokud sledujeme vliv jednotlivých stromů, odpovídá minimální vzdálenost nejdelšímu stínu daného stromu. Kromě sledovaných kategorií by se jednotlivé lokality neměly lišit v ostatních charakteristikách prostředí - například nadmořské výšce, expozici, zástinu atd.

Následující příklad čerpá z dat získaných v roce 2024 v Průhonickém parku. Cílem bylo srovnat vliv konifer na mikroklima podrostu. Abychom eliminovali vliv topografie, která v rámci Průhonického parku vytváří rozdíly denních průměrů až 12 °C (Brůna et al., 2023), založili jsme pokus v rámci jedné topoklimatické jednotky „mírné severní svahy“ (Wild et al., 2014).

Následující druhy jsme našli v dostačujícím počtu a kvalitě, abychom mohli vyhodnotit jejich vliv na mikroklima: Abies alba, Abies grandis, Picea abies a Pseudotsuga menziesii. Vybrali jsme jednotlivé stromy 5-9 metrů vysoké, pro každou dřevinu 8 opakování. Ze severní strany ke kmeni každého stromu jsme umístili mikroklimatickou stanici TOMST TMS-4. K ověření vlivu porostu konifer na mikroklima obecně bylo 8 mikroklimatických stanic umístěno v rámci blízké mladé výsadby na otevřené nestíněné ploše.

Metodika

Začneme načtením informací o souborech z tabulky files_table.csv, kterou je třeba si připravit dle informací v kapitole 1.4

## Načtení tabulky 'files_table'
ft <- read.table("./data/dreviny/files_table.csv", sep=";", header = TRUE)
# Pro načtení z podsložky je třeba upravit sloupec path (relativní cesta vůči aktuální složce)
ft$path<-paste0("./data/dreviny/data/",ft$path)
# Načtení dat loggerů s metadaty
tms <- mc_read_data(files_table = ft, silent = TRUE)

Měření probíhalo od 15.12.2023 do 22.7.2024. Pro příklad jsme vybrali jaro 21.3.-20.6.2024. Pro tento úsek jsme spočítali průměrnou, minimální a maximální teplotu a objemovou vlhkost půdy. Pro výpočty využívající vlhkost půdy je důležité správně nastavit druh půdy na každém stanovišti, případně lze využít přímo surová data ze senzoru, pokud je půda všude stejná.

## Oříznutí časové řady na vybrané období - jaro
start <- as.POSIXct("2024-03-21", tz = "UTC")
end <- as.POSIXct("2024-06-21", tz = "UTC")
tms.jaro <- mc_prep_crop(tms, start, end, end_included=FALSE)
filename_prefix<-"jaro_"

## Výpočet virtuálního senzoru VWC ze surových dat vlhkosti TMS
tms.jaro <- mc_calc_vwc(tms.jaro, soiltype = "loamy sand A")

## Výpočet vegetačních a mrazových stupňů (GDD a FDD)
tms.jaro <- mc_calc_gdd(tms.jaro, sensor = "TMS_T3")
tms.jaro <- mc_calc_fdd(tms.jaro, sensor = "TMS_T3")

## Agregace všech časových řad, výpočet průměrů, rozptylů 
tms.jaro.vse <- mc_agg(tms.jaro, fun = c("mean", "range", "coverage", 
                                         "percentile", "min", "max", "sum"),
                         percentiles = 95, period = "all", min_coverage = 0.95)

## Převod na dlouhý formát (pro lepší manipulaci s daty)
tms.jaro.vse.long <- mc_reshape_long(tms.jaro.vse)

## Převod na široký formát pro snadnější analýzu
df_wide.jaro <- tms.jaro.vse.long %>%
  select(-height) %>%
  pivot_wider(names_from = sensor_name, values_from = value)

## Načtení dat o lokalitách a odstranění lokality "Plaste_40"
lokality_df <- read.csv("./data/dreviny/lokality.csv", sep=";")
lokality_df <- lokality_df[lokality_df$locality_id != "Plaste_40",]
# Spojení s daty
df_wide.jaro <- lokality_df %>%
  left_join(df_wide.jaro, by = "locality_id")

V létě 2024 byly pořízeny hemisférické fotografie všech lokalit ve výšce 80 cm nad zemí a odvozen korunový zápoj pomocí knihovny HemisphereR (Chianucci a Macek 2023). Zápoj byl hodnocen pro celou hemisféru i pro jednotlivé výseče po 10°. Díky tomu můžeme ověřit rozdíly mikroklimatu způsobené korunovým zápojem. Z terénních zápisů víme, že lokalita Plaste_40 byla zničena a data z ní nejsou dostupná, proto ji vyřadíme i z tabulky lokalit.

V tabulce lokalit máme též sloupec “drevina”, kde je uvedeno, pod kterým stromem byla která stanice TMS-4. Stejně tak zde mohou být uvedeny další relevantní proměnné pro analýzu, například nadmořská výška, souřadnice, svažitost, přítomnost druhu aj.

Náhled tabulky lokalit včetně metadat.
locality_id lat_wgs84 lon_wgs84 drevina GapFraction Openness
Plaste_01 49.99036 14.56065 Picea abies 4.97 4.97
Plaste_02 49.99040 14.56075 Abies alba 5.08 5.05

Pro vizualizaci vytvoříme krabicové diagramy (boxplot), které znázorňují medián, kvartilové rozpětí a odlehlé hodnoty naměřených veličin pro jednotlivé druhy a kontrolu.

Statistickou průkaznost na 95 % hladině významnosti zjistíme analýzou variance s post-hoc testem (Tukey HSD) pro srovnání jednotlivých druhů a kontroly (více v Crawley 2013). Signifikance rozdílů je vyznačena přímo v grafu pomocí písmen (CLD, Compact letter display). Boxploty se stejným písmenem nejsou signifikantně rozdílné, naopak ty, které nemají žádné shodné písmeno, jsou signifikantně odlišné. Rozdíly, které nejsou signifikantní, mohou být způsobeny variabilitou v datech. Pro zvýšení šance na odhalení signifikatních rozdílů je třeba zvětšit počet nezávislých měření.

Výsledky

Nejprve si vygenerujeme graf pro jednu proměnnou - průměrnou teplotu pod zemí (TMS_T1_mean).

# Proměnná pro analýzu
promenna <- "TMS_T1_mean"

# ANOVA
model <- aov(as.formula(paste(promenna, "~ drevina")), data = df_wide.jaro)
tukey <- TukeyHSD(model)

# Kompaktní zobrazení písmen (CLD)
cld <- multcompLetters4(model, tukey)$drevina
cld_df <- data.frame(drevina = names(cld$Letters), Letters = cld$Letters)

# Úprava názvů druhů pro dvouřádkové zobrazení
df_wide.jaro$drevina <- gsub(" ", "\n", df_wide.jaro$drevina)

# Spojení s maximální hodnotou pro správné umístění písmen v grafu
max_values <- df_wide.jaro %>%
  group_by(drevina) %>%
  summarize(max_value = max(.data[[promenna]], na.rm = TRUE))
cld_df <- merge(cld_df, max_values, by = "drevina")

# Vytvoření boxplotu
ggplot(df_wide.jaro, aes(x = drevina, y = .data[[promenna]], fill = drevina)) +
  geom_boxplot() +
  geom_text(
    data = cld_df,
    aes(label = Letters, y = max_value + 0.1 * max_value),
    position = position_dodge(width = 0.75),
    vjust = -0.5
  ) +
  labs(
    x = "Druh",
    y = promenna
  ) +
  theme_minimal() +
  theme(legend.position = "none", axis.text.x = element_text(face = "italic"))
Boxplot pro jednotlivé druhy s průměrnou teplotou pod zemí na jaře.

Boxplot pro jednotlivé druhy s průměrnou teplotou pod zemí na jaře.

Pro průzkum rozdílů v jednotlivých klimatických proměnných mezi jednotlivými druhy se hodí si vygenerovat grafy pro jednotlivé proměnné. Následující kód generuje soubory do složky data/dreviny/vysledky.

# Seznam proměnných k analýze - pro každou vznikne i graf
variables <-
  c("TMS_T1_mean","TMS_T3_mean","TMS_T3_min", 
    "TMS_T3_max", "TMS_T3_range","TMS_T3_percentile95",
    "VWC_moisture_mean","GDD5_sum","FDD0_sum","Openness")

# Funkce pro provedení ANOVA, Tukey HSD testu a vytvoření boxplotů
analyzuj_data <- function(df, promenna, prefix) {

# ANOVA
  model <- aov(as.formula(paste(promenna, "~ drevina")), data = df)
  tukey <- TukeyHSD(model)
  
# Kompaktní zobrazení písmen (CLD)
  cld <- multcompLetters4(model, tukey)$drevina
  cld_df <-data.frame(drevina = names(cld$Letters),
           Letters = cld$Letters)
  
# Spojení s maximální hodnotou pro správné umístění písmen v grafu
  max_values <- df %>%
    group_by(drevina) %>%
    summarize(max_value = max(.data[[promenna]], na.rm = TRUE))
  cld_df <- merge(cld_df, max_values, by = "drevina")
  
# Vytvoření boxplotu
  p <-
    ggplot(df, aes(x = drevina, y = .data[[promenna]], fill = drevina)) +
    geom_boxplot() +
    geom_text(
      data = cld_df,
      aes(label = Letters, y = max_value + 0.1 * max_value),
      position = position_dodge(width = 0.75),
      vjust = -0.5
    ) +
    labs(
      title = paste("Boxplot", prefix, promenna),
      x = "Druh",
      y = promenna
    ) +
    theme_minimal() +
    theme(legend.position = "none",
          axis.text.x = element_text(face = "italic"))
  
# Export boxplotu ve formátu .png
  ggsave(paste0( "./data/dreviny/vysledky/",
      prefix, promenna,"_boxplot.png"),
    plot = p, width = 8, height = 6 , bg = "white")
}

for (promenna in variables) {
  analyzuj_data(df_wide.jaro, promenna, filename_prefix)
}

V tomto příkladě byly grafy vygenerovány do souborů, které najdete v části soubory ve složce data/dreviny/vysledky. Níže jsou vložené přímo jako obrázky.

Korunový zápoj jednotlivých druhů se v žádném ze sledovaných parametrů signifikantně nelišil, kontrola byla vždy odlišná. Rozdíly v mikroklimatu tak lze přisoudit především odlišnostem jednotlivých druhů.

Otevřenost lesního porostu spočtená z hemisférických fotografií slouží pro ověření, že lokality se liší jen druhem a že mezi zástinem není rozdíl.
Otevřenost lesního porostu spočtená z hemisférických fotografií slouží pro ověření, že lokality se liší jen druhem a že mezi zástinem není rozdíl.

Většina klimatických parametrů ve výšce 15 cm nad zemí byla na kontrolních plochách signifikantně odlišná, vliv porostu konifer na mikroklima je nepopiratelný.

Jednotlivé stromy dokáží pufrovat extrémy teploty 15 cm nad zemí. V rámci hodnot pro jaro jde o snižování maximálních teplot až o 5 °C a zvyšování minimálních teplot o 2 °C.

Teplotní parametry v půdě se signifikantně nelišily mezi druhy a ve většině případů ani mezi porostem a kontrolou.

Boxploty průměrných půdních teplot dle jednotlivých dřevin.
Boxploty průměrných půdních teplot dle jednotlivých dřevin.

Zjištěné rozdíly v teplotních parametrech ve výšce 15 cm nad zemí mezi druhy nejsou velké a pouze v několika parametrech byly signifikantně odlišné. Největší rozdíly byly v maximálních teplotách.

Maximální teploty v jarním období pod stromy Abies alba se signifikantně liší od stromů Pseudotsuga menziesii a od kontroly. Teploty pod všemi druhy se na jaře signifikantně liší od kontroly.

Boxploty maximálních jarních teplot v 15 cm nad zemí pro lokality rozdělené podle jednotlivých druhů.
Boxploty maximálních jarních teplot v 15 cm nad zemí pro lokality rozdělené podle jednotlivých druhů.

V ostatních parametrech nejsou signifikantní rozdíly, ale i tak je viditelný vyšší pufrační vliv především Abies alba a Picea abies. Pod stromy Abies grandis a Pseudotsuga menziesii jsou často vyšší teploty.

Objemová vlhkost půdy je signifikantně odlišná jen mezi kontrolou a druhy, mezi jednotlivými druhy není v tomto parametru signifikantní rozdíl. Na kontrole je během jara vlhkost půdy o 15 procentních bodů vyšší. To je pravděpodobně dáno především intercepcí srážek na větvích stromů a využitím vody kořeny stromů.

Boxploty průměrné objemové půdní vlhkosti pro lokality rozdělené podle dřevin během jara.
Boxploty průměrné objemové půdní vlhkosti pro lokality rozdělené podle dřevin během jara.

3.2. Srovnání s daty z meteorologických stanic ČHMÚ

Úvod

Často nás zajímá, zda se mikroklima naší lokality liší od nejbližší meteorologické stanice. V červnu 2024 publikovalo ČHMÚ denní data svých stanic ve strojově čitelném formátu JSON. V datech jsou obsažena i  historická data ve formě denních průměrů. Pro teploty jde o SYNOP průměr, který je počítaný z hodnot v 6:00, 13:00 a 20:00 UTC. Pro jednoduchost budeme tyto hodnoty porovnávat s denními průměry všech měření ze stanic TMS.

Pro příklady v části 3.2 a 3.3 načítáme data z připraveného souboru, jelikož jinak je třeba stáhnout několik GB dat a ty filtrovat. Na konci části 3.3 přikládáme skript, kterým jsme tato data získali a upravili. Při hledání relevantní stanice ČHMÚ je třeba zohlednit nejen její vzdálenost, ale i podobnou nadmořskou výšku a zda měří požadovanou proměnnou. Pokud by nadmořská výška byla výrazně odlišná, lze předpokládat, že velká část rozdílu teplot je způsobená především tímto rozdílem.

Metodika

Pro dvě lokality na Šumavě najdeme nejbližší meteorologickou stanici. Použijeme připravený soubor, kde jsou pouze stanice, které měřily teplotu v roce 2022 a navíc mají vypočtenou vzdálenost od naší stanice. Úpravou skriptu na konci části 3.3 můžete stáhnout data i pro jiné stanice dle potřebné proměnné. Nejbližší stanice je na Horské Kvildě s  WSI: 0-203-0-11494 a na Churáňově, s WSI kódem: 0-20000-0-11457. Churáňov má podobnější nadmořskou a pro další zpracování vybereme pouze tuto stanici.

# Najděte nejbližší stanice v zadaném poloměru - souřadnice WGS84 a poloměr v kilometrech
stanice_s_informacemi<-readRDS("./data/stanice_s_informacemi.rds")
max_vzdalenost<-5 # 5 km
stanice_v_blizkosti<-stanice_s_informacemi[stanice_s_informacemi$vzdalenost_km
                                           <max_vzdalenost,]
print(stanice_v_blizkosti$WSI)
## [1] "0-20000-0-11457" "0-203-0-11494"
# Vybereme jen Churáňov, jelikož má podobnější nadmořskou výšku
stanice_v_blizkosti<-stanice_s_informacemi[stanice_s_informacemi$WSI=="0-20000-0-11457",]

# Načtení připravených dat
climate_data<-readRDS("./data/chmu_data_2022.rds")

# Filtrování climate_data na základě odpovídajících hodnot WSI 
climate_data <- climate_data[climate_data$STATION %in% stanice_v_blizkosti$WSI, ]

Nyní načteme naše 2 lokality - tentokrát ve formátu TMS_join - v této podobě často dodáváme zprocesovaná data, která mohou pocházet z kombinace čidel, která byla na lokalitě v průběhu několika let. V tomto příkladu využíváme data z Thermologgeru umístěného ve výšce 2 m nad zemí na severní straně stromu.

# Načtení souborů stanic TOMST Thermologger pomocí files_table
ft <- read.table("./data/stanice/files_table.csv", sep=";", header = TRUE)
ft$path<-paste0("./data/stanice/",ft$path)

# Načtení dat loggerů s metadaty
tms <- mc_read_data(files_table = ft, silent = TRUE)

# Oříznutí časové řady na vybrané období
start <- as.POSIXct("2022-01-01", tz = "UTC")
end <- as.POSIXct("2023-01-01", tz = "UTC")
tms <- mc_prep_crop(tms, start, end, end_included=FALSE)

# Agregace všech časových řad, výpočet průměrů, rozptylů 
tms.agg <- mc_agg(tms, fun = c("mean"), period = "day", min_coverage = 0.95)

# Převod na dlouhý formát (pro lepší manipulaci s daty)
tms.agg.long <- mc_reshape_long(tms.agg)

## Převod na široký formát pro snadnější analýzu
df_wide.tms <- tms.agg.long %>%
  select(-height) %>%
  pivot_wider(names_from = sensor_name, values_from = value)

# Uložení pro pozdější využití (zde už ukládá tabulku, proto stačí 
# saveRDS a nemusím volat myClim::mc_save)
saveRDS(df_wide.tms, file = "./data/dve_tms_sumava.rds")

Výsledky

Nyní vykreslíme jednoduchý graf s daty ze stanice Churáňov a z našich mikroklimatických stanic.

# Kombinace datasetů do jednoho s přidaným sloupcem 'locality'
combined_data <- rbind(
  data.frame(Date = climate_data$DT, 
             Temp = as.numeric(climate_data$VAL), Locality = "CHMU"),
  data.frame(Date = as.Date(df_wide.tms[df_wide.tms$locality_id==
                                          "NPS_4210_D_T1",]$datetime), 
             Temp = df_wide.tms[df_wide.tms$locality_id=="NPS_4210_D_T1",]$Thermo_T_mean, 
             Locality = "NPS_4210_D_T1"),
  data.frame(Date = as.Date(df_wide.tms[df_wide.tms$locality_id=="4220_TMS",]$datetime), 
             Temp = df_wide.tms[df_wide.tms$locality_id=="4220_TMS",]$Thermo_T_mean, 
             Locality = "4220_TMS")
)

# Vykreslení pomocí palety snadno rozpoznatelné i pro barvoslepé (paleta Okabe-Ito)
ggplot(combined_data, aes(x = Date, y = Temp, color = Locality)) +
  geom_line() +
  scale_color_manual(values = c("#E69F00", "#000000", "#009E73")) +  
  labs(title = "Průměrná denní teplota ze tří zdrojů",
       x = "Datum", y = "Teplota (°C)",
       color = "Lokalita") +  
  theme_minimal()
Na první pohled je vidět, že průměrné teploty mají podobné průběhy a že na stanicích TMS jsou většinou nižší.

Na první pohled je vidět, že průměrné teploty mají podobné průběhy a že na stanicích TMS jsou většinou nižší.

Zajímavý je velký rozdíl například v březnu, na který se podíváme dále.

# Definujte rozsah dat pro výběr
start_date <- as.Date("2022-03-01")
end_date <- as.Date("2022-03-31")

filtered_data <- combined_data[combined_data$Date >= start_date
                               & combined_data$Date <= end_date,]

# Vytvoření grafu pro výběr
ggplot(filtered_data, aes(x = Date, y = Temp, color = Locality)) +
  geom_line() +
  scale_color_manual(values = c("#E69F00", "#000000", "#009E73")) +
  labs(title = "Průměrná denní teplota ze tří zdrojů",
       x = "Datum", y = "Teplota (°C)",
       color = "Lokalita"
      ) +
  theme_minimal()
V některých dnech se průměrné teploty lišily až o 10°C. Během teplých dní byly průměrné teploty na stanicích TMS nižší, naopak během chladnějších byly průměrné teploty na stanicích TMS vyšší.

V některých dnech se průměrné teploty lišily až o 10°C. Během teplých dní byly průměrné teploty na stanicích TMS nižší, naopak během chladnějších byly průměrné teploty na stanicích TMS vyšší.

Pro snazší vyhodnocení se můžeme podívat na histogram rozdílu teplot mezi každou mikroklimatickou stanicí a stanicí na Churáňově.

# Spojení datasetů GSOD, NPS_4210_D_T1 a 4220_TMS podle data
merged_data <- combined_data %>%
  spread(Locality, Temp)

# Výpočet denních rozdílů
merged_data <- merged_data %>%
  mutate(Difference_NPS_CHMU = NPS_4210_D_T1 - CHMU,
         Difference_TMS_CHMU = `4220_TMS` - CHMU)

# Získání rozsahu osy x a y pro oba grafy
x_limits <- range(c(merged_data$Difference_NPS_CHMU,
                    merged_data$Difference_TMS_CHMU),
            na.rm = TRUE)

y_limits <- c(0, max(table( cut(merged_data$Difference_NPS_CHMU,
            breaks = seq(min(x_limits), max(x_limits), by = 0.5))),
            table(cut(merged_data$Difference_TMS_CHMU,
            breaks = seq(min(x_limits), max(x_limits), by = 0.5)))))

# Histogram pro rozdíl mezi NPS_4210_D_T1 a GSOD
p1 <- ggplot(merged_data, aes(x = Difference_NPS_CHMU)) +
  geom_histogram(binwidth = 0.5, fill = "#E69F00",
                 color = "black") +
  geom_vline(xintercept = 0, linetype = "dashed",
             color = "red") +  # Přidání čáry na x=0
  scale_x_continuous(limits = x_limits) +  # Stejná osa x pro oba grafy
  scale_y_continuous(limits = y_limits) +  # Stejná osa y pro oba grafy
  labs(title = "NPS_4210_D_T1 - ČHMÚ",
       x = "Rozdíl teplot (°C)", y = "Počet dnů") +
  theme_minimal()

# Histogram pro rozdíl mezi 4220_TMS a CHMU
p2 <- ggplot(merged_data, aes(x = Difference_TMS_CHMU)) +
  geom_histogram(binwidth = 0.5, fill = "#56B4E9",
                 color = "black") +
  geom_vline(xintercept = 0, linetype = "dashed",
             color = "red") +  # Přidání čáry na x=0
  scale_x_continuous(limits = x_limits) +  # Stejná osa x pro oba grafy
  scale_y_continuous(limits = y_limits) +  # Stejná osa y pro oba grafy
  labs(title = "4220_TMS - ČHMÚ",
       x = "Rozdíl teplot (°C)", y = "Počet dnů") +
  theme_minimal()

# Zobrazení grafů vedle sebe
grid.arrange(p1, p2, ncol = 2)
Histogram rozdílů průměrných denních teplot jednotlivých stanic TMS a průměrných denních teplot ze stanice ČHMÚ na Churáňově.

Histogram rozdílů průměrných denních teplot jednotlivých stanic TMS a průměrných denních teplot ze stanice ČHMÚ na Churáňově.

Dobrou pomůckou pro vizualizaci může být i scatterplot, kdy na ose X vidíme staniční data a na ose Y data z mikroklimatických stanic.

# Získání minimální a maximální hodnoty pro osy X a Y (zahrnutí 0 a stejné limity pro obě osy)
min_limit <- min(c(merged_data$CHMU, 
                   merged_data$NPS_4210_D_T1,
                   merged_data$`4220_TMS`), 
                 na.rm = TRUE)
max_limit <- max(c(merged_data$CHMU, 
                   merged_data$NPS_4210_D_T1,
                   merged_data$`4220_TMS`),
                 na.rm = TRUE)
axis_limits <- c(min(0, min_limit), max_limit)

# Vytvoření scatterplotu
ggplot(merged_data, aes(x = CHMU)) +
  # Body pro NPS_4210_D_T1 s průhledností
  geom_point(aes(y = NPS_4210_D_T1, color = "NPS_4210_D_T1"), alpha = 0.6) + 
  # Body pro 4220_TMS s průhledností
  geom_point(aes(y = `4220_TMS`, color = "4220_TMS"), alpha = 0.6) +           
  # Čára 1:1
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "red") +  
  # Černá čára na Y = 0
  geom_hline(yintercept = 0, color = "black", size = 0.5) +                    
  # Černá čára na X = 0
  geom_vline(xintercept = 0, color = "black", size = 0.5) +       
  # Stejné limity osy X a Y
  scale_x_continuous(limits = axis_limits) +  
  scale_y_continuous(limits = axis_limits) +  
  labs(title = "Porovnání stanic ČHMÚ s NPS_4210_D_T1 a 4220_TMS",
       x = "Teplota CHMU (°C)", y = "Teplota (°C)",
       color = "Lokalita") +  # Popisek legendy
  theme_minimal() +
  scale_color_manual(values = c("NPS_4210_D_T1" = "#E69F00", 
                                "4220_TMS" = "#56B4E9")) +  
  theme(axis.line = element_blank(),)  # Odstranění osových čar
Při použití scatterplotu je dobré sledovat pozici oproti linii 1:1. Díky tomu můžeme vidět, že během chladných i teplých dní byly průměry ze stanic TMS spíše nižší, ale výjimkou nejsou i vyšší hodnoty, nejen v chladné části. Na stanici NPS_4210_D_T1 jsou průměry většinou nižší než data z Churáňova.

Při použití scatterplotu je dobré sledovat pozici oproti linii 1:1. Díky tomu můžeme vidět, že během chladných i teplých dní byly průměry ze stanic TMS spíše nižší, ale výjimkou nejsou i vyšší hodnoty, nejen v chladné části. Na stanici NPS_4210_D_T1 jsou průměry většinou nižší než data z Churáňova.

Pro snazší předání výsledku se hodí jednoduchý výpočet průměrného rozdílu denních teplot. Při studiu mikroklimatu se často používá “pufrovací efekt” (“buffering effect”) lokality.

# Výpočet průměrného rozdílu (pufrovací efekt) a směrodatné odchylky
buffering_effect <- merged_data %>%
  summarise(
    Mean_Buffering_NPS_4210_D_T1 = mean(Difference_NPS_CHMU, na.rm = TRUE),
    SD_Buffering_NPS_4210_D_T1 = sd(Difference_NPS_CHMU, na.rm = TRUE),
    Mean_Buffering_4220_TMS = mean(Difference_TMS_CHMU, na.rm = TRUE),
    SD_Buffering_4220_TMS = sd(Difference_TMS_CHMU, na.rm = TRUE)
  )

# Vytvoření textového výstupu pro obě lokality
cat(sprintf("Pufrovací efekt vegetace na lokalitě NPS_4210_D_T1 je %.2f ± %.2f °C\n", 
            buffering_effect$Mean_Buffering_NPS_4210_D_T1,
            buffering_effect$SD_Buffering_NPS_4210_D_T1))
## Pufrovací efekt vegetace na lokalitě NPS_4210_D_T1 je -1.13 ± 1.33 °C
cat(sprintf("Pufrovací efekt vegetace na lokalitě 4220_TMS je %.2f ± %.2f °C\n", 
            buffering_effect$Mean_Buffering_4220_TMS,
            buffering_effect$SD_Buffering_4220_TMS))
## Pufrovací efekt vegetace na lokalitě 4220_TMS je -0.57 ± 0.81 °C

Je vidět, že efekt není příliš velký a má velkou směrodatnou odchylku. Testování signifikance takového rozdílu je složitější. Běžně používaný párový T-test předpokládá, že rozdíly mezi páry jsou nezávislé a mají normální rozdělení. Pokud jednotlivé dny nejsou nezávislé (například teploty mezi po sobě jdoucími dny mohou být silně korelované kvůli sezónním a denním trendům), může dojít k porušení předpokladu nezávislosti, což ovlivní výsledky testu. To může vést k:

  • podhodnocení variability: Pokud jsou hodnoty závislé (např. teploty v po sobě jdoucích dnech), může to vést k podhodnocení variability, což může způsobit, že test bude falešně detekovat statisticky významné rozdíly.

  • větší pravděpodobnost chyby I. druhu: Může dojít k vyšší pravděpodobnosti chyby I. druhu, tedy že zamítneme nulovou hypotézu (že rozdíl neexistuje), i když je pravdivá.

Jak dál?

Správným postupem je pak například využití ARIMA modelů. (AutoRegressive Integrated Moving Average) (Box et Jenkins, 1970). Tyto modely jsou speciálně navrženy k zachycení vzorů v časových řadách, kde jsou hodnoty mezi jednotlivými dny závislé na předchozích hodnotách.

Jak ARIMA funguje:

- AR (AutoRegressive): Část AR modelu předpokládá, že aktuální hodnota je lineární kombinací předchozích hodnot (lagů).

- I (Integrated): Pomáhá odstranit trend nebo sezónní komponentu pomocí diferenciace (odečítání sousedních hodnot).

- MA (Moving Average): MA model zohledňuje závislost na předchozích chybách modelu.

To je však již nad rámec této metodiky.

3.3. Srovnání se stanicemi v celé ČR

Úvod

Zde využijeme opět knihovnu GSODR, ale podíváme se, jak si stojí naše stanice oproti teplotám v celé České republice na gradientu nadmořské výšky.

Metodika

Vybereme si jeden den a pro ten stáhneme data ze stanic v ČR (stahují se po celých letech, následně se filtrují)

# Zadejte datum zájmu
date_of_interest <- as.Date("2022-09-01")

# Načtení připravených dat
stanice_s_informacemi<-readRDS("./data/stanice_s_informacemi.rds")
climate_data<-readRDS("./data/chmu_data_2022.rds")

# Spojení informací o stanicích s daty o teplotě
climate_data <- climate_data %>%
  left_join(stanice_s_informacemi, by = c("STATION" = "WSI"), 
            relationship = "many-to-many")

# Úprava dat
climate_data$VAL<-as.numeric(climate_data$VAL)
climate_data$ELEVATION<-as.numeric(climate_data$ELEVATION)

# Filtrování dat pro konkrétní datum
filtered_data <- subset(climate_data, DT == date_of_interest)

A nyní načteme data ze stejných mikroklimatických stanic jako v minulé ukázce. Navíc doplníme jejich nadmořskou výšku, abychom je mohli v grafu srovnat se stanicemi v ČR a vyfiltrujeme pro vybraný den.

# Načtení dat
df_wide.tms<-readRDS("./data/dve_tms_sumava.rds")
# Přidání nadmořské výšky
df_wide.tms$elevation <- 1120
df_wide.tms$mean_temp <- df_wide.tms$Thermo_T_mean
# Filtrování na vybraný den
filtered_tms <- subset(df_wide.tms, as.Date(datetime) == date_of_interest)

Výsledky

A teď již vytvoříme graf pro vybraný den (1.9.2022), kam přidáme i jednoduchou regresi teploty a nadmořské výšky s intervaly spolehlivosti.

# Kontrola, zda jsou pro dané datum dostupná data
if (nrow(filtered_data) > 0) {
  # Vytvořit scatter plot teplota vs. nadmořská výška s trendovou čarou
  p <- ggplot(filtered_data, aes(x = ELEVATION, y = VAL)) +
    geom_point(color = "blue", alpha = 0.7) +
    # Přidání regresní čáry a intervalů spolehlivosti
    geom_smooth(method = "lm", se = TRUE, color = "black", linetype = "dashed",
                fill = "lightblue", alpha = 0.3) +  
    labs(title = paste("Teplota vs. nadmořská výška dne", date_of_interest, 
                       "zdroj: ČHMÚ"),
         x = "Nadmořská výška (m n. m.)", y = "Průměrná denní teplota (°C)") +
    theme_minimal()
  
  # Přidání dat ze stanic loggerů jako červené body
  p <- p + 
    geom_point(data = filtered_tms, aes(x = elevation, y = mean_temp), 
               color = "red", size = 3)
  # Zobrazení grafu
  print(p)
} else {
  cat("Pro vybrané datum nejsou dostupná data:", date_of_interest, "\n")
}
Je vidět, že 1.9.2022 průměrná teplota na stanici nevybočovala z hodnot v rámci ČR. Z grafu je ale také patrné, že v horských polohách je stanic málo.

Je vidět, že 1.9.2022 průměrná teplota na stanici nevybočovala z hodnot v rámci ČR. Z grafu je ale také patrné, že v horských polohách je stanic málo.

A jak to bylo v jiných dnech? Využijeme toho, že data už máme načtená, pouze vybereme (filtrujeme) jiný den.

# Zadejte datum zájmu
date_of_interest <- as.Date("2022-01-12")

# Filtrovat data pro konkrétní datum
filtered_data <- subset(climate_data, DT == date_of_interest)
filtered_tms <- subset(df_wide.tms, as.Date(datetime) == date_of_interest)

# Zkontrolovat, zda jsou pro dané datum dostupná data
if (nrow(filtered_data) > 0) {
  # Vytvořit scatterplot teplota vs. nadmořská výška s trendovou čarou
  p <- ggplot(filtered_data, aes(x = ELEVATION, y = VAL)) +
    geom_point(color = "blue", alpha = 0.7) +
    # Přidání regresní čáry a intervalů spolehlivosti
    geom_smooth(method = "lm", se = TRUE, color = "black", linetype = "dashed",
                fill = "lightblue", alpha = 0.3) +  
    labs(title = paste("Teplota vs. nadmořská výška dne", date_of_interest, 
                       "zdroj: ČHMÚ"),
         x = "Nadmořská výška (m n. m.)", y = "Průměrná denní teplota (°C)") +
    theme_minimal()
  
  # Přidání dat ze stanic loggerů jako červené body
  p <- p + 
    geom_point(data = filtered_tms, aes(x = elevation, y = mean_temp), 
               color = "red", size = 3)
  # Zobrazení grafu
  print(p)
} else {
  cat("Pro vybrané datum nejsou dostupná data:", date_of_interest, "\n")
}
12.1.2022 byla teplota na našich stanicích nižší, než bychom čekali dle stanic v ČR v této nadmořské výšce, tomu odpovídá i hodnota na Churáňově. Můžeme spekulovat, zda šlo o vliv sněhu v lese nebo šlo jen o prostý zástin během slunečného zimního dne.

12.1.2022 byla teplota na našich stanicích nižší, než bychom čekali dle stanic v ČR v této nadmořské výšce, tomu odpovídá i hodnota na Churáňově. Můžeme spekulovat, zda šlo o vliv sněhu v lese nebo šlo jen o prostý zástin během slunečného zimního dne.

A jak to vypadá v případě teplotní inverze? k té došlo například 31.10.2022.

# Datum zájmu
date_of_interest <- as.Date("2022-10-31")

# Filtrovat data pro konkrétní datum
filtered_data <- subset(climate_data, DT == date_of_interest)
filtered_tms <- subset(df_wide.tms, as.Date(datetime) == date_of_interest)

# Zkontrolovat, zda jsou pro dané datum dostupná data
if (nrow(filtered_data) > 0) {
  # Vytvořit scatter plot teplota vs. nadmořská výška s trendovou čarou
  p <- ggplot(filtered_data, aes(x = ELEVATION, y = VAL)) +
    geom_point(color = "blue", alpha = 0.7) +
    # Přidání regresní čáry a    intervalů spolehlivosti
    geom_smooth(method = "lm", se = TRUE, color = "black", linetype = "dashed",
                fill = "lightblue", alpha = 0.3) +  
    labs(title = paste("Teplota vs. nadmořská výška dne", date_of_interest, 
                       "zdroj: ČHMÚ"),
         x = "Nadmořská výška (m n. m.)",y = "Průměrná denní teplota (°C)") +
    theme_minimal()
  # Přidání dat ze stanic loggerů jako červené body
  p <- p + 
    geom_point(data = filtered_tms, aes(x = elevation, y = mean_temp), 
               color = "red", size = 3)
  # Zobrazení grafu
  print(p)
} else {
  cat("Pro vybrané datum nejsou dostupná data:", date_of_interest, "\n")
}
I v tomto případě byla teplota na našich stanicích nižší, než by odpovídalo meteorologickým stanicím na bezlesí.

I v tomto případě byla teplota na našich stanicích nižší, než by odpovídalo meteorologickým stanicím na bezlesí.

Skript pro stažení a filtrování dat ČHMÚ

Následuje skript, kterým jsme připravili podklady pro kapitolu 3.2 a 3.3. Modifikací lze stáhnout data pro jiné stanice, proměnné nebo časový úsek.

# Seznam potřebných knihoven
required_packages <- c("jsonlite", "lubridate", "sf", "dplyr", "stringr")

# Instalace a načtení knihoven
suppressWarnings(
  lapply(required_packages, function(pkg) {
    if (!require(pkg, character.only = TRUE)) install.packages(pkg, 
                               repos = "https://mirrors.nic.cz/R/")
    library(pkg, character.only = TRUE)
  })
)
rm(required_packages)

# Definice funkcí
# Funkce pro stažení metadat a filtrování stanic na základě dostupné proměnné a roku
ziskat_stanice_s_dostupnou_promennou <- function(rok, element) {
  # Stažení metadatového souboru meta2.json pomocí curl
  metadata_url <- 
    "https://opendata.chmi.cz/meteorology/climate/historical/metadata/meta2.json"
  metadata_file <- tempfile(fileext = ".json")
  system(paste("curl -o", metadata_file, metadata_url))
  # Načtení metadat ze souboru
  metadata <- fromJSON(metadata_file)   
  # Zpracování metadat, extrahování hodnot a filtrování podle proměnné a roku
  metadata_df <- data.frame(metadata$data$data$values)
  colnames(metadata_df) <- unlist(strsplit(metadata$data$data$header, ","))
  # Filtrování stanic, které mají dostupná data pro zvolený rok a proměnnou (element)
  stanice_s_promennou <- metadata_df %>%
    filter(EG_EL_ABBREVIATION == element) %>%
    mutate(BEGIN_DATE = as.Date(BEGIN_DATE),
           END_DATE = as.Date(END_DATE),
           HEIGHT=as.numeric(HEIGHT)) %>%
    filter(year(BEGIN_DATE) <= rok & year(END_DATE) >= rok & HEIGHT < 2.5 & HEIGHT > 1.9)
  return(stanice_s_promennou)
}

# Funkce pro načtení meta1.json a připojení informací o stanicích s filtrací dle roku
pripojit_meta1 <- function(stanice_s_promennou, rok) {
  # Stažení meta1.json pomocí download.file
  metadata1_url <- 
    "https://opendata.chmi.cz/meteorology/climate/historical/metadata/meta1.json"
  metadata1_file <- tempfile(fileext = ".json")
  download.file(metadata1_url, metadata1_file)
  metadata1 <- fromJSON(metadata1_file)
  # Extrakce informací z meta1.json a vytvoření datového rámce
  metadata1_df <- data.frame(metadata1$data$data$values)
  colnames(metadata1_df) <- c("WSI", "GH_ID", "BEGIN_DATE", "END_DATE", "FULL_NAME", 
                              "GEOGR1", "GEOGR2", "ELEVATION")
  # Převod BEGIN_DATE a END_DATE na Date formát
  metadata1_df <- metadata1_df %>%
    mutate(BEGIN_DATE = as.Date(BEGIN_DATE),
           END_DATE = as.Date(END_DATE))
  # Filtrování na základě zadaného roku a zajištění, že GEOGR1 a GEOGR2 nejsou NA
  metadata1_df_filtered <- metadata1_df %>%
    filter(year(BEGIN_DATE) <= rok & year(END_DATE) >= rok) %>%
    filter(!is.na(GEOGR1) & !is.na(GEOGR2))
  # Připojení informací z meta1 k tabulce stanice_s_promennou pomocí left_join
  stanice_s_informacemi <- stanice_s_promennou %>%
    left_join(metadata1_df_filtered, by = "WSI")
  stanice_s_informacemi <- stanice_s_informacemi %>%
    filter(!is.na(GEOGR1) & !is.na(GEOGR2))
  return(stanice_s_informacemi)
}

# Funkce pro výpočet vzdálenosti v kilometrech v UTM33N
vypocet_vzdalenosti_utm33n <- function(lat1, lon1, lat2, lon2) {
  # Vytvoření sf objektů pro body
  bod1 <- st_as_sf(data.frame(lon = lon1, lat = lat1), 
                   coords = c("lon", "lat"), crs = 4326)
  bod2 <- st_as_sf(data.frame(lon = lon2, lat = lat2), 
                   coords = c("lon", "lat"), crs = 4326)
  # Transformace souřadnic do UTM33N (EPSG:32633)
  bod1_utm <- st_transform(bod1, 32633)
  bod2_utm <- st_transform(bod2, 32633)
  # Výpočet vzdálenosti v metrech a převedení na kilometry
  vzdalenost_m <- st_distance(bod1_utm, bod2_utm)
  vzdalenost_km <- as.numeric(vzdalenost_m) / 1000
  return(vzdalenost_km)
}

# Funkce pro filtrování stanic na základě vzdálenosti od zadaného bodu
filtrovat_stanice_dle_vzdalenosti <- function(stanice_df, lat, lon, max_vzdalenost_km) {
  # Výpočet vzdálenosti pro každou stanici
  stanice_df <- stanice_df %>%
    rowwise() %>%
    mutate(vzdalenost_km = vypocet_vzdalenosti_utm33n(lat, lon, as.numeric(GEOGR2), 
                                                      as.numeric(GEOGR1))) %>%
    ungroup() %>%
    # Filtrování podle vzdálenosti
    filter(vzdalenost_km <= max_vzdalenost_km)  
  return(stanice_df)
}

# Funkce pro načítání a čištění JSON souboru s jsonlite
ziskat_data_stanice_pro_rok <- function(stanice_id, element, vtype, rok) {
  # Složka pro ukládání dat
  dir.create("./data/stanice", recursive = TRUE, showWarnings = FALSE)
  filename <- paste0("./data/stanice/dly-", stanice_id, ".json")
  
  # Pokud soubor existuje, načti ho, jinak ho stáhni
  if (!file.exists(filename)) {
    # Upravený formát URL pro stahování bez uvedení roku
    data_url <- 
      paste0("https://opendata.chmi.cz/meteorology/climate/historical/data/daily/dly-", stanice_id, ".json")
    # Stáhnout soubor pomocí curl
    download_command <- paste("curl -w '%{http_code}' -o", shQuote(filename), data_url)
    status_code <- system(download_command, intern = TRUE)
    
    # Zajisti, že kontrolujeme jen první hodnotu status_code
    if (length(status_code) == 0 || status_code[1] != "200") {
      message("Soubor není dostupný nebo nelze stáhnout (HTTP status: ", status_code, "): ", data_url)
      return(NULL)
    }
  }
  
  # Pokus o načtení JSON dat pomocí jsonlite s chybovým ošetřením
  data <- tryCatch({
    fromJSON(filename, simplifyVector = TRUE)  # Načítání JSON souboru pomocí jsonlite
  }, error = function(e) {
    # Pokud dojde k chybě, vypíše chybu a přeskočí soubor
    message("Chyba při načítání JSON souboru pomocí jsonlite: ", filename, " - ", e$message)
    return(NULL)  # Přeskočí tento soubor a vrátí NULL
  })
  
  if (is.null(data)) {
    return(NULL)  # Pokud JSON soubor nebyl načten, přeskoč soubor
    }
  
  # Kontrola, zda JSON obsahuje data
  if (!"data" %in% names(data)) {
    message("Chybí data ve staženém souboru: ", filename)
    return(NULL)
  }
  
  # Pokus o převod dat na datový rámec
  data_filtered <- tryCatch({
    as.data.frame(data$data$data$values)
  }, error = function(e) {
    message("Chyba při převodu dat na data frame: ", e$message)
    return(NULL)
  })
  
  if (is.null(data_filtered)) {
    return(NULL)  # Přeskočí soubor, pokud nelze převést data na dataframe
  }

  # Přiřazení názvů sloupců
  colnames(data_filtered) <- unlist(strsplit(data$data$data$header, ","))
  
  # Filtrování podle ELEMENT a VTYPE
  data_filtered <- data_filtered %>%
    filter(ELEMENT == element, VTYPE == vtype)
  
  # Převod DT na datum a filtrování pro celý rok 2022
  data_filtered$DT <- as.Date(data_filtered$DT)
  data_filtered <- data_filtered %>%
    filter(year(DT) == rok)
  
  # Přidání identifikátoru stanice
  data_filtered$stanice_id <- stanice_id
  
  return(data_filtered)
}

# Funkce pro získání dat pro všechny stanice za rok 2022
ziskat_data_ze_vsech_stanic <- function(stanice_s_informacemi, element, vtype, rok) {
  # Inicializace prázdného dataframe
  vysledna_data <- data.frame()  
  for (i in 1:nrow(stanice_s_informacemi)) {
    stanice_id <- stanice_s_informacemi$WSI[i]
    message("Zpracovávám stanici s WSI: ", stanice_id)
    data <- ziskat_data_stanice_pro_rok(stanice_id, element, vtype, rok)
    if (!is.null(data)) {
      # Přidání dat pro každou stanici
      vysledna_data <- rbind(vysledna_data, data)  
    }
  }
  return(vysledna_data)
}

Nyní nastavujeme proměnné: souřadnice, sledovanou proměnnou a její typ, maximální vzdálenost a rok. V tomto případě jsme nastavili vzdálenost 750 km, abychom dostali všechny stanice, které v roce 2022 měřily teplotu.

# Zadání souřadnic, roku zájmu, ELEMENT, VTYPE, maximální vzdálenosti a konkrétního datumu
lat <- 49.0679850   # Zeměpisná šířka
lon <- 13.5678442   # Zeměpisná délka
element <- "T"      # Parametr ELEMENT, např. "T" (teplota)
vtype <- "AVG"      # Parametr VTYPE, např. "AVG" (průměrná hodnota)
max_vzdalenost_km <- 750  # Maximální vzdálenost v kilometrech
rok <- 2022

# Získání seznamu stanic, které mají dostupná data pro daný rok a proměnnou
stanice_s_promennou <- ziskat_stanice_s_dostupnou_promennou(rok, element)

# Připojení informací z meta1 s filtrováním na základě roku a dostupných souřadnic
stanice_s_informacemi <- pripojit_meta1(stanice_s_promennou, rok)
stanice_s_informacemi <- stanice_s_informacemi[!is.na(stanice_s_informacemi$GEOGR1) & 
                                                 !is.na(stanice_s_informacemi$GEOGR2), ]

# Filtrování stanic na základě vzdálenosti od zadaného bodu
stanice_s_informacemi <- filtrovat_stanice_dle_vzdalenosti(stanice_s_informacemi, 
                                                           lat, lon, max_vzdalenost_km)

# Uložit seznam stanic, které vyhovují podmínkám včetně jejich metadat
saveRDS(stanice_s_informacemi, file = "./data/stanice_s_informacemi.rds")

Nyní už stačí stáhnout data pro vybrané stanice a uložit je do tabulky chmu_data_2022.

# Získání dat o teplotě pro všechny stanice za celý rok 2022
chmu_data_2022 <- ziskat_data_ze_vsech_stanic(stanice_s_informacemi, element, vtype, rok)

# Uložení dat do souboru ve formátu RDS pro další zpracování
saveRDS(chmu_data_2022, file = "./data/chmu_data_2022.rds")

3.4. Srovnání s daty ve světě

Denní data z různých stanic celého světa lze získat pomocí knihovny GSODR. Data pochází z GSOD neboli Global Surface Summary of the Day. Data jsou poskytovaná Národními centry pro environmentální informace USA (NCEI) a jsou cenným zdrojem meteorologických dat s celosvětovým pokrytím. V tuto chvíli je zahrnuto 40 stanic v České republice. Kromě dostupnosti zahraničních dat je výhodou i rychlost zpracování bez nutnosti stahovat velké množství JSON souborů. Tímto způsobem jsou dostupná data z cca 9000 stanic, z nichž některá mají data od roku 1929. Globální pokrytí stanic je vidět na Obrázku 32.

Mapa stanic v databázi GSOD.
Mapa stanic v databázi GSOD.

Metodika

Opět zkusíme stáhnout data ze stanice na Churáňově, jejíž ID v databázi GSOD je “114570-99999”. To použijeme pro získání denních dat z roku 2022. Opět načteme data ze dvou našich stanic a připojíme k datům z GSOD.

# Najděte nejbližší stanice v zadaném poloměru - souřadnice WGS84 a poloměr v kilometrech
GSODR::nearest_stations(LAT = 49.0679850, LON = 13.5678442, distance = 20)

# Stáhněte data pro zadanou lokalitu a období
station_id <- "114570-99999"  # Churáňov
start_year <- 2022
end_year <- 2022

# Ošetření pokusu o stažení dat - GSOD bývá občas nedostupný
tryCatch({
  # Pokusíme se stáhnout data z GSODR
  climate_data_gsod <- GSODR::get_GSOD(start_year:end_year, station = station_id)
}, error = function(e) {
  # V případě chyby načteme připravená data ze souboru
  message("Nepodařilo se stáhnout data, načítám uložená data ze souboru.")
  climate_data_gsod <- readRDS("./data/gsod_Churanov_2022.rds")
})

# Načtení uložených TMS stanic (viz kapitola 3.2)
df_wide.tms<-readRDS("./data/dve_tms_sumava.rds")

# Kombinace datasetů do jednoho s přidaným sloupcem 'locality'
combined_data <- rbind(
  data.frame(Date = as.Date(climate_data_gsod$YEARMODA), 
             Temp = climate_data_gsod$TEMP, Locality = "GSOD Churáňov"),
  data.frame(Date = as.Date(df_wide.tms[df_wide.tms$locality_id==
                                          "NPS_4210_D_T1",]$datetime), 
             Temp = df_wide.tms[df_wide.tms$locality_id=="NPS_4210_D_T1",]$Thermo_T_mean, 
             Locality = "NPS_4210_D_T1"),
  data.frame(Date = as.Date(df_wide.tms[df_wide.tms$locality_id==
                                          "4220_TMS",]$datetime), 
             Temp = df_wide.tms[df_wide.tms$locality_id==
                                  "4220_TMS",]$Thermo_T_mean, 
             Locality = "4220_TMS")
)

Výsledky

Nyní se podíváme na průběh průměrné teploty v březnu 2022.

# Definujte rozsah dat pro výběr
start_date <- as.Date("2022-03-01")
end_date <- as.Date("2022-03-31")
filtered_data <- combined_data[combined_data$Date >= start_date
                               & combined_data$Date <= end_date,]

# Vytvoření grafu pro výběr
ggplot(filtered_data, aes(x = Date, y = Temp, color = Locality)) +
  geom_line() +
  scale_color_manual(values = c("#E69F00", "#000000", "#009E73")) +
  labs(
    title = "Průměrná denní teplota ze tří zdrojů",
    x = "Datum",
    y = "Teplota (°C)",
    color = "Lokalita"
  ) +  # Popisek legendy
  theme_minimal()
Zde je vidět, že data v databázi GSOD nejsou shodná s daty ČHMÚ. Částečně to je dáno jiným způsobem výpočtu průměru, částečně to může být další úpravou dat v rámci GSOD.

Zde je vidět, že data v databázi GSOD nejsou shodná s daty ČHMÚ. Částečně to je dáno jiným způsobem výpočtu průměru, částečně to může být další úpravou dat v rámci GSOD.

Případně můžeme stáhnout data pro celý stát a filtrovat pro vybraný den jako v části 3.3.

# Zadejte datum zájmu
date_of_interest <- as.Date("2022-09-01")

# Stáhnout data pro rok vybraný podle data zájmu pro Českou republiku
year_of_interest <- format(date_of_interest, "%Y")
climate_data_cz <- GSODR::get_GSOD(as.integer(year_of_interest),
                            station = NULL,
                            country = "CZE",
                            max_missing = NULL,
                            agroclimatology = FALSE)

# Filtrovat data pro konkrétní datum
filtered_data <- subset(climate_data_cz, YEARMODA == date_of_interest)

Zajímavější však bude srovnat data ze stanic TMS s daty stanic v okruhu 250 km, namísto celé ČR, jelikož lokalita leží blízko hranic.

# Období
start_year <- 2022
end_year <- 2022

# Najděte nejbližší stanice v zadaném poloměru (v km)
stations_nearby <- GSODR::nearest_stations(LAT = 49.0679850, LON = 13.5678442, 
                                           distance = 250)

# Je vidět, že v údajích o stanicích jsou chyby, ty očividné odstraníme
stations_nearby<-stations_nearby[!stations_nearby$COUNTRY_NAME=="GUAM",]
stations_nearby<-stations_nearby[!stations_nearby$COUNTRY_NAME=="IRAN",]
stations_nearby<-stations_nearby[!stations_nearby$COUNTRY_NAME=="IRAQ",]
stations_nearby<-stations_nearby[!is.na(stations_nearby$`ELEV(M)`),]

# Stáhnutí dat GSOD pro všechny stanice v okolí 250 km
# Použijeme funkci lapply pro stažení dat z více stanic
climate_data_nearby <- lapply(stations_nearby$STNID, function(station_id) {
  get_GSOD(years = start_year:end_year, station = station_id)
})

# Sloučení všech stažených dat do jednoho datového rámce
climate_data_nearby <- do.call(rbind, climate_data_nearby)

# Pro zrychlení práce uložíme data, tato část kódu se nespouští automaticky, 
# aby nedocházelo k přetěžování úložiště, pokud změníte poloměr, spusťte znovu. 
# Aktuální nastavení stahuje 291 stanic.

saveRDS(climate_data_nearby, file = "./data/gsod_nearby_2022.rds")

A nyní vytvoříme opět graf pro vybraný den (1.9.2022), kam přidáme i jednoduchou regresi teploty a nadmořské výšky s intervaly spolehlivosti. Tentokrát jsou však využita data z okolních stanic v rámci GSOD.

# Načtení uložených TMS stanic (viz kapitola 3.2)
df_wide.tms<-readRDS("./data/dve_tms_sumava.rds")
climate_data_nearby<-readRDS("./data/gsod_nearby_2022.rds")
# Přidání nadmořské výšky
df_wide.tms$elevation <- 1120

# Zadejte datum zájmu
date_of_interest <- as.Date("2022-09-01")

# Filtrování dat pro konkrétní datum
filtered_data <- subset(climate_data_nearby, as.Date(YEARMODA) == date_of_interest)
filtered_data_tms <- subset(df_wide.tms, as.Date(datetime) == date_of_interest)

# Kontrola, zda jsou pro dané datum dostupná data
if (nrow(filtered_data) > 0) {
  # Vytvoření scatter plotu teplota vs. nadmořská výška s trendovou čarou
  p <- ggplot(filtered_data, aes(x = ELEVATION, y = TEMP)) +
    geom_point(color = "blue", alpha = 0.7) +
    # Přidání regresní čáry a intervalů spolehlivosti
    geom_smooth(method = "lm", se = TRUE, color = "black", linetype = "dashed",
                fill = "lightblue", alpha = 0.3) +  
    labs(title = paste("Teplota vs. nadmořská výška dne", date_of_interest, "zdroj: GSOD"),
         x = "Nadmořská výška (m n. m.)", y = "Průměrná denní teplota (°C)") +
    theme_minimal()
  # Přidání dat ze stanic loggerů jako červené body
  p <- p + geom_point(data = filtered_data_tms, aes(x = elevation, y = Thermo_T_mean), 
           color = "red", size = 3)
  print(p)
} else {
  cat("Pro vybrané datum nejsou dostupná data:", date_of_interest, "\n")
}
Přestože graf vypadá podobně, z rozložení stanic podél gradientu nadmořské výšky je vidět, že jde o jiné stanice.

Přestože graf vypadá podobně, z rozložení stanic podél gradientu nadmořské výšky je vidět, že jde o jiné stanice.

4. Práce s rastry

V současnosti je mnoho klimatických dat dostupných ve formě rastrů, od globálních až po ryze lokální. Liší se mezi sebou nejen rozlišením (pixel může mít hranu v řádu metrů i kilometrů), ale také časovým rozsahem dat a typem proměnných. Pro mnoho otázek tak můžeme získat relevantní data, aniž bychom museli instalovat klimatické stanice. Vždy je však třeba posoudit, zda jsou daná data vhodná k řešení dané otázky, především s ohledem na rozlišení, ale i zdroj dat. V této kapitole ukážeme několik příkladů rastrů a jak s nimi pracovat v prostředí R. Pro rastry je často jednodušší využít geografický informační systém (GIS), například QGIS.

4.1. Rastry mikroklimatu pro NP Šumava zpracované BÚ

Pro NP Šumava jsme v roce 2021 připravili soubor mikroklimatických map pokrývajících území NP Šumava a NP Bavorský les. Data jsou veřejně publikovaná (Brůna et al., 2023) a lze s nimi pracovat přímo v prostředí R. Mapy vznikly s využitím dat ze sítě mikroklimatických stanic z období říjen 2019 - říjen 2020. Efekt vegetačního pokryvu a topografie byl zahrnut s využitím dat laserového skenování (LiDAR). Modely pro interpolaci zohledňují nejen topografické vlivy na mikroklima, ale i vliv hustoty a typu porostu. Mapy se hodí pro řešení detailních otázek týkajících se klimatu v lesních porostech a na přilehlém bezlesí ve 3 různých výškách - pod zemí, 15 cm nad zemí a 200 cm nad zemí.

Příprava potřebných knihoven prostředí R.

V úložišti Zenodo jsou dostupné (název souboru můžeme využít v následujícím kódu):

Teplota půdy (-8 cm)

  • Průměrná teplota = průměrná teplota T.soil_8_cm.mean.tif

Teplota vzduchu u země (15 cm)

  • Průměrná teplota T.air_15_cm.mean.tif

  • Maximální teplota = 95. percentil denních maximálních teplot T.air_15_cm.max.95p.tif

  • Minimální teplota = 5. percentil denních minimálních teplot T.air_15_cm.min.5p.tif

Teplota vzduchu (200 cm)

  • Průměrná teplota T.air_200_cm.mean.tif

  • Maximální teplota = 95. percentil denních maximálních teplot T.air_200_cm.max.95p.tif

  • Minimální teplota = 5. percentil denních minimálních teplot T.air_200_cm.min.5p.tif

  • GDD = součet stupňů růstu nad základní teplotou (základ 5°C) T.air_200_cm.GDD5.tif

Pomocí následujícího kódu stáhneme mapu průměrné teploty ve výšce 15 cm nad zemí a načteme ji do objektu zeta.

# URL ke stažení souboru (soubor T.air_15_cm.mean.tif z archivu Zenodo, velikost 97,0 MB)
soubor <- c("T.air_15_cm.mean.tif")

url <- paste0("https://zenodo.org/record/6352641/files/",soubor)

# Cesta, kam uložit soubor
destfile <- paste0("./data/rastry/", soubor)

# Vytvoření složky, pokud neexistuje
dir.create("./data/rastry", recursive = TRUE, showWarnings = FALSE)

# download rasteru pomocí curl
curl::multi_download(url, destfile, resume = T)

# Načtení rastru pomocí knihovny terra
zeta <- terra::rast(destfile)

Mapa

Po stažení vizualizujeme celou mapu pomocí základního grafu.

# Nastavení prostoru pro legendu
par(mar = c(5, 4, 4, 5))  # Přidání prostoru na pravé straně pro legendu
par(oma = c(2, 1, 1, 1))  # Zvýšení vnějšího spodního okraje

# Základní graf rastru pomocí funkce plot()
plot(zeta, col = viridis::viridis(100), 
     xlab = "Zeměpisná délka", ylab = "Zeměpisná šířka")

# Přidání popisku k legendě s jednotkou (°C)
mtext("Teplota (°C)", side = 4, line = 3)  # Přidání popisku na pravou stranu
Mapa roční průměrné teploty vzduchu 15 cm nad zemí na území NP Šumava (10-2019 - 10-2020). Mapa ukazuje výrazný vliv nadmořské výšky na teplotu. Patrný je též vliv topografie.

Mapa roční průměrné teploty vzduchu 15 cm nad zemí na území NP Šumava (10-2019 - 10-2020). Mapa ukazuje výrazný vliv nadmořské výšky na teplotu. Patrný je též vliv topografie.

Histogram

Nyní se podíváme na rozložení teplot pomocí histogramu, který zobrazí počet výskytů různých hodnot.

# Výpočet průměrné hodnoty rastru
values_zeta <- values(zeta)
mean_zeta <- mean(values_zeta, na.rm = TRUE)

# Vytvoření histogramu pro hodnoty rastru
hist_zeta <- hist(values_zeta, main = "Histogram teplot", 
                  xlab = "Průměrná roční teplota (°C)", ylab = "Četnost výskytů",
                  col = "lightblue", breaks = 20, 
                  probability = TRUE, border = "gray")

# Přidání vertikální čáry označující průměr
abline(v = mean_zeta, col = "red", lwd = 2)

# Přidání popisku s hodnotou průměru
text(mean_zeta, par("usr")[4] * 0.92,
     labels = paste0("Průměr: ", round(mean_zeta, 2)," °C"),
     col = "red", pos = 4)
Na histogramu vidíme, že roční průměr je 6.09 °C, ale v různých oblastech se může pohybovat od 4.4°C do 8.6°C.

Na histogramu vidíme, že roční průměr je 6.09 °C, ale v různých oblastech se může pohybovat od 4.4°C do 8.6°C.

Výřez

Nyní se podíváme na menší oblast, abychom využili vysokého rozlišení map. Pro příklad ořízneme mapu na čtverec 1000 x 1000 m v okolí Březnické hájenky.

# Definování souřadnic středu ve WGS84( EPSG:4326)
lat <- 48.9712522 # Březnická hájenka
lon <- 13.4825903
coords <- data.frame(lon, lat)

# Vytvoření sf objektu se souřadnicemi v EPSG:4326
point_wgs84 <- st_as_sf(coords, coords = c("lon", "lat"), crs = 4326)

# Transformace bodu do souřadného systému odpovídajícímu rastrovým podkladům
# (3-degree Gauss-Kruger zone 4, EPSG:31468)
point_epsg31468 <- st_transform(point_wgs84, 31468)

# Získání souřadnic pro vytvoření obdélníku (bounding box) kolem bodu
point_coords <- st_coordinates(point_epsg31468)
x_center <- point_coords[1, "X"]
y_center <- point_coords[1, "Y"]

# Definování velikosti výřezu (např. 500x500 metrů kolem bodu)
buffer_size <- 500  # polovina velikosti (v metrech)
xmin <- x_center - buffer_size
xmax <- x_center + buffer_size
ymin <- y_center - buffer_size
ymax <- y_center + buffer_size

# Vytvoření obdélníku (bbox)
bbox <- terra::ext(xmin, xmax, ymin, ymax)

# Oříznutí rastru pomocí obdélníku
zeta_cropped <- crop(zeta, bbox)

# Nastavení prostoru na okrajích mapy pro legendu
par(mar = c(4, 4, 4, 5))  # Přidání prostoru na pravé straně pro legendu
par(oma = c(2, 1, 1, 1))  # Zvýšení vnějšího spodního okraje

# Základní graf rastru pomocí funkce plot()
plot(zeta_cropped, col = viridis::viridis(100), 
     xlab = "Zeměpisná délka (m), EPSG:31468", 
     ylab = "Zeměpisná šířka (m)")

# Přidání popisku k legendě s jednotkou (°C)
mtext("Teplota (°C)", side = 4, line = 3)  # Přidání popisku na pravou stranu

# Přidání bodu se souřadnicemi na mapu (symbol 'x')
points(x_center, y_center, pch = 4, col = "red", cex = 1.5)  # pch=4 je symbol 'x'
Výřez: na této mapě jsou patrné teplejší, k jihu orientované luční porosty pod Březnickou hájenkou (v mapě označená x) a chladnější údolní polohy. Opět jde o roční průměrnou teplotu (10-2019 - 10-2020).

Výřez: na této mapě jsou patrné teplejší, k jihu orientované luční porosty pod Březnickou hájenkou (v mapě označená x) a chladnější údolní polohy. Opět jde o roční průměrnou teplotu (10-2019 - 10-2020).

Získání hodnoty z rastrů pro body

Hodnoty ve formě rastru jsou dobrým podkladem pro tvorbu map, častěji však potřebujeme hodnoty pro zájmové body jako podklad pro analýzu. Následující ukázka nejprve převede body do souřadnicového systému mapy, následně pro body extrahuje hodnoty a přidá je do tabulky.

# Načtení souřadnic lat/lon z CSV, které obsahuje sloupce 'lat' a 'lon'
coords_data <- read.csv("data/rastry/souradnice.csv", sep=";")

# Vytvoření sf objektu ze souřadnic v EPSG:4326 (WGS84)
points_wgs84 <- st_as_sf(coords_data, coords = c("lon", "lat"), crs = 4326)

# Transformace bodů do EPSG:31468
points_epsg31468 <- st_transform(points_wgs84, 31468)

# Převod sf objektu na formát vhodný pro terra::extract
coords_matrix <- st_coordinates(points_epsg31468)

# Extrakce hodnot z rastru pro souřadnice pomocí bilineární interpolace
values <- terra::extract(zeta, coords_matrix, method = "bilinear")

# Spojení extrahovaných hodnot se souřadnicemi a uložení do nové tabulky
final_data <- cbind(coords_data, values)
default_kable(final_data)
Poslední bod má ve sloupci s teplotou hodnotu NA, patrně byl mimo oblast mapy. To můžeme snadno ověřit vykreslením bodů do mapy.
id lat lon altitude T.air_15_cm.mean
1 49.11978 13.30771 1140.50 5.333770
2 49.07077 13.38393 1132.40 5.249397
3 48.98027 13.63936 1015.00 5.057910
4 49.05457 13.49517 907.00 6.836538
5 48.98433 13.47457 1207.10 4.736646
6 48.99739 13.52126 1136.30 4.898316
7 48.92075 13.66820 902.20 4.993112
8 48.87504 13.76450 1055.60 5.354843
9 48.96517 13.49173 1176.43 5.175743
10 48.98573 13.79490 1247.00 NaN

Body též můžeme vizualizovat přímo v mapě spolu s rastrem. Tím si navíc ověříme správnost bodů.

# Nastavení prostoru pro legendu a graf
par(mar = c(4, 4, 4, 5))  # Přidání prostoru na pravé straně pro legendu
par(oma = c(2, 1, 1, 1))  # Zvýšení vnějšího spodního okraje

# Vykreslení rastru - pro kontrolu, že všechny body jsou uvnitř mapy
plot(zeta, main = "Vizualizace rastru se souřadnicemi", 
     col = viridis(100),  # Použití barevné škály viridis
     xlab = "Zeměpisná délka (m), EPSG:31468", 
     ylab = "Zeměpisná šířka (m)")

# Přidání legendy s jednotkou (°C)
mtext("Teplota (°C)", side = 4, line = 3)

# Přidání všech bodů ze souřadnic jako křížky (pch=4 je symbol 'x')
points(coords_matrix[,1], coords_matrix[,2], pch = 4, col = "red", cex = 1.5)
Zde je vidět, že jeden z bodů byl opravdu mimo oblast mapy.

Zde je vidět, že jeden z bodů byl opravdu mimo oblast mapy.

4.2. Externí zdroje

ForestClim

Nyní se podíváme na externí zdroj ForestClim (Haesen et al., 2023). Jedná se o soubor map bioklimatologických proměnných pokrývajících lesy Evropy s rozlišením 25 m. Soubor map je založen na měření mikroklimatickými dataloggery, převážně TOMST TMS-4. Využití tedy opět cílí na rozdíly teploty v porostech. Vysoké rozlišení umožňuje porovnávat i velmi blízké lokality. V tomto případě jsou dostupné bioklimatologické proměnné, které sumarizují mikroklima pomocí 11 odvozených proměnných:

ForestClim_01 = průměrná roční teplota

ForestClim_02 = průměrný denní rozsah teplot

ForestClim_03 = izotermalita

ForestClim_04 = sezónnost teplot

ForestClim_05 = maximální teplota nejteplejšího měsíce

ForestClim_06 = minimální teplota nejchladnějšího měsíce

ForestClim_07 = roční rozsah teplot

ForestClim_08 = průměrná teplota v nejvlhčím čtvrtletí

ForestClim_09 = průměrná teplota v nejsušším čtvrtletí

ForestClim_10 = průměrná teplota v nejteplejším čtvrtletí

ForestClim_11 = průměrná teplota v nejchladnějším čtvrtletí

Jednotlivé soubory mají velikost okolo 6 GB, a tak je snazší si data stáhnout předem zde: https://figshare.com/articles/dataset/ForestClim_Bioclimatic_variables_for_microclimate_temperatures_of_European_forests/22059125?file=39164684

Pro náš příklad využijeme průměrnou roční teplotu, tedy vrstvu ForestClim_01. Do složky s daty jsme uložili výřez pro oblast NP Šumava. Následující ukázka opět vygeneruje jednoduchou mapu.

# Načtení rastru ForestClim pomocí terra
forest_raster <- terra::rast("./data/rastry/ForestClim_01_CZE_NPS.tif")

# Zobrazení nového rastru
par(mar = c(5, 4, 4, 5))  # Nastavení okrajů
plot(forest_raster, main = paste0("Vizualizace rastru ", "ForestClim_01"), 
     col = viridis::viridis(100), 
     xlab = "Zeměpisná délka", ylab = "Zeměpisná šířka")

# Přidání popisku k legendě
mtext("Teplota (°C)", side = 4, line = 3)  # Přidání popisku na pravou stranu
Vizualizace výřezu rastru ForestClim_01 pro oblast NP Šumava

Vizualizace výřezu rastru ForestClim_01 pro oblast NP Šumava

Opět se podíváme na histogram průměrných ročních teplot.

# Histogram pro ForestClim raster
forest_mean_value <- mean(values(forest_raster), na.rm = TRUE)

# Histogram pro ForestClim
hist_forest <- hist(values(forest_raster), main = "Histogram ForestClim", 
                    xlab = "Průměrná roční teplota (°C)", ylab = "Počet výskytů",
                    col = "lightgreen", breaks = 20, 
                    probability = TRUE, border = "gray")

# Přidání vertikální čáry označující průměr pro ForestClim raster
abline(v = forest_mean_value, col = "darkred", lwd = 2)

# Přidání popisku pro průměrnou hodnotu
text(forest_mean_value, par("usr")[4] * 0.9, labels = 
     paste0("Průměr: ", round(forest_mean_value, 2),"°C"), 
     col = "darkred", pos = 4)
Přestože průměr je hodně podobný, počty výskytů hodnot v tomto rastru jsou jiné.

Přestože průměr je hodně podobný, počty výskytů hodnot v tomto rastru jsou jiné.

A opět se podíváme na detail mapy okolo Březnické hájenky

# Využíváme stejný point_wgs84, definovaný v předchozí ukázce (Březnická hájenka)
# Transformace bodu do EPSG:3035 (European Lambert Conformal Conic)
point_epsg3035 <- st_transform(point_wgs84, 3035)

# Získání souřadnic pro vytvoření obdélníku (bounding box) kolem bodu
point_coords <- st_coordinates(point_epsg3035)
x_center <- point_coords[1, "X"]
y_center <- point_coords[1, "Y"]

# Definování velikosti výřezu (např. 500 metrů kolem bodu)
buffer_size <- 500  
xmin <- x_center - buffer_size
xmax <- x_center + buffer_size
ymin <- y_center - buffer_size
ymax <- y_center + buffer_size

# Vytvoření obdélníku (bbox)
bbox <- terra::ext(xmin, xmax, ymin, ymax)

# Oříznutí rastru pomocí obdélníku (pro ForestClim raster)
forest_raster_cropped <- crop(forest_raster, bbox)

# Nastavení prostoru pro legendu
par(mar = c(4, 4, 4, 5))  # Přidání prostoru na pravé straně pro legendu
par(oma = c(2, 1, 1, 1))  # Zvýšení vnějšího spodního okraje

# Základní graf rastru pomocí funkce plot()
plot(forest_raster_cropped, main = paste0("Vizualizace výřezu rastru ForestClim"),
     col = viridis::viridis(100),
     xlab = "Zeměpisná délka (m), EPSG:3035", 
     ylab = "Zeměpisná šířka (m)")

# Přidání popisku k legendě
mtext("Hodnoty rastru", side = 4, line = 3)  # Popisek na pravé straně

# Přidání bodu se souřadnicemi na mapu (symbol 'x')
points(x_center, y_center, pch = 4, col = "red", cex = 1.5)  # pch=4 je symbol 'x'
Zde je vidět slabina tohoto datasetu. Kromě nižšího prostorového rozlišení (větších pixelů) je to pokrytí pouze lesních porostů, které jsou definovány vrstvami Copernicus (více v Haesen et al., 2023), a tudíž disturbované plochy, zmlazení a nižší porosty nejsou zahrnuty.

Zde je vidět slabina tohoto datasetu. Kromě nižšího prostorového rozlišení (větších pixelů) je to pokrytí pouze lesních porostů, které jsou definovány vrstvami Copernicus (více v Haesen et al., 2023), a tudíž disturbované plochy, zmlazení a nižší porosty nejsou zahrnuty.

Srovnání histogramů rastrů Forestclim a BÚ

Nabízí se srovnání s rastrem průměrné teploty, který vytvořil Botanický ústav. Tato část spoléhá na proměnné, které jste vytvořili při spuštění kódu v části 4.1.

# Hodnoty pro oba rastry
values_forest <- values(forest_raster)

# Výpočet společných breaků histogramu
breaks <- seq(min(c(values_zeta, values_forest), na.rm = TRUE), 
              max(c(values_zeta, values_forest), na.rm = TRUE), length.out = 30)

# Vykreslení prvního histogramu pro T.air_15_cm.mean.tif (v paměti z předchozích výpočtů)
plot(hist_zeta, col = rgb(0, 0, 1, 0.5), xlim = range(breaks), 
     ylim = c(0, max(hist_zeta$counts)), 
     main = "Kombinovaný histogram rastrů s dvojí osou Y", 
     xlab = "Hodnota", ylab = soubor, 
     cex.main = 1.2)  # Zmenšení velikosti titulku

# Přidání čáry a popisku průměru pro T.air_15_cm.mean.tif (průměr nad čárou)
abline(v = mean_zeta, col = "blue", lwd = 2)
text(mean_zeta, max(hist_zeta$counts) * 0.75, labels = 
       paste0("Průměr: ", round(mean_zeta, 2),"°C"), col = "blue", pos = 4)

# Vykreslení druhé osy y na pravé straně
par(new = TRUE)
plot(hist_forest, col = rgb(0, 1, 0, 0.5), xlim = range(breaks), 
     ylim = c(0, max(hist_forest$counts)), 
     axes = FALSE, xlab = "", ylab = "", main="")
axis(side = 4)
mtext("Frekvence ForestClim", side = 4, line = 3)

# Přidání čáry a popisku průměru pro ForestClim (průměr pod čárou)
abline(v = forest_mean_value, col = "darkgreen", lwd = 2)
text(forest_mean_value, max(hist_forest$counts) * 0.5, labels = 
       paste0("Průměr: ", round(forest_mean_value, 2),"°C"), col = "darkgreen", pos = 4)

# Přidání legendy
legend("topright", legend = c("T.air_15_cm.mean", "ForestClim 01"), 
       fill = c(rgb(0, 0, 1, 0.5), rgb(0, 1, 0, 0.5)))
Přestože průměr obou datasetů se příliš neliší, histogram ukazuje, že mikroklimatická mapa zpracovaná Botanickým ústavem má více chladnějších hodnot, než je průměr a zároveň místa s vyššími teplotami. Oboje patrně souvisí s ořezem ForestClim na les, nicméně může to ukazovat i větší vliv porostů na mikroklima v datasetu vytvořeném Botanickým ústavem.

Přestože průměr obou datasetů se příliš neliší, histogram ukazuje, že mikroklimatická mapa zpracovaná Botanickým ústavem má více chladnějších hodnot, než je průměr a zároveň místa s vyššími teplotami. Oboje patrně souvisí s ořezem ForestClim na les, nicméně může to ukazovat i větší vliv porostů na mikroklima v datasetu vytvořeném Botanickým ústavem.

Downscalovaná data E-OBS pomocí knihovny easyclimate

Knihovna easyclimate zpřístupňuje denní klimatická data s rozlišením 1 km pro Evropu (srážky, minimální a maximální teploty) z evropské klimatické databáze hostované na Univerzitě přírodních zdrojů a biologických věd, Vídeň, Rakousko. Data jsou v současné době k dispozici od roku 1950 do roku 2022.

Vstupní klimatický dataset je založen na downscalovaných rastrech E-OBS (European Observations). Původně byl vytvořen A. Moreno a H. Hasenauerem - https://doi.org/10.1002/joc.4436 a dále vyvíjen W. Rammerem, C. Pucherem a M. Neumannem. Aktuální verze je v4 - detailně popsaná v tomto dokumentu - https://doi.org/10.6084/m9.figshare.22962671.v1). Pro podrobný popis knihovny easyclimate si přečtěte tento článek - https://doi.org/10.1016/j.envsoft.2023.105627 (verze s otevřeným přístupem https://doi.org/10.32942/osf.io/mc8uj) nebo navštivte webové stránky knihovny - https://verughub.github.io/easyclimate/.

Vstupní rastry E-OBS jsou založené na interpolaci dat z 7852 klimatických stanic v Evropě. Díky tomu jsou vhodné pro charakterizaci makroklimatu bez ohledu na vliv lesa, případně jiného krajinného pokryvu. Díky downscalování bylo zvýšeno rozlišení na 1 km, které ale stále nedokáže podchytit mikroklimatické rozdíly.

Příklad stažení rastru pro vybranou oblast. Výhodou je, že můžeme stahovat i více dní najednou. Načteme hranici NP Šumava v souřadnicovém systému UTM 33N, převedeme do souřadnicového systému zdroje dat (WGS 84) a stáhneme data pro 1.-3.5.2020.

## Načtení shapefilu z určené cesty (EPSG:32633)
nps_region <- vect("./data/rastry/NPS_UTM33N.shp")

## Převod shapefilu na požadovaný souřadnicový systém (EPSG:4326 pro souřadnice lon lat)
nps_region <- project(nps_region, "EPSG:4326")

## Stažení hodnot Tmax pro oblast NP Šumava mezi 1. a 3. květnem 2020
nps_temp <- get_daily_climate(
  coords = nps_region,
  climatic_var = "Tmax", # lze změnit na "Prcp" pro srážky nebo "Tmin" pro min. teploty
  # průměrná teplota zde není, ale běžně se počítá jako (Tmax+Tmin)/2
  period = "2020-05-01:2020-05-03",
  output = "raster"
)

Výstup nps_temp je objekt SpatRaster s 3 vrstvami (pro každý ze 3 dnů). Nyní si zobrazíme mapy pro jednotlivé dny.

# Vykreslení jednotlivých vrstev a přidání hranic polygonu nps_region
ggplot() +
  geom_spatraster(data = nps_temp) +
  facet_wrap(~lyr, ncol = 3) +
  scale_fill_distiller(palette = "RdYlBu", na.value = "transparent") +
  geom_spatvector(data = nps_region, fill = NA) +
  labs(fill = "Maximum\ntemperature \n(ºC)") +
  scale_x_continuous(breaks = seq(13.3, 13.9, by = 0.3)) +
  scale_y_continuous(breaks = seq(48.8, 49.2, by = 0.2)) +
  theme_minimal()
Na třech mapách je vidět maximální teploty pro jednotlivé dny 1.-3.5.2020.

Na třech mapách je vidět maximální teploty pro jednotlivé dny 1.-3.5.2020.

Díky tomu, že data jsou dostupná od roku 1950, můžeme je dobře využít pro srovnání.

## Stažení hodnot Tmax pro oblast NPS mezi 1. a 3. květnem 1950
nps_temp_1950 <- get_daily_climate(
  coords = nps_region,
  climatic_var = "Tmax",
  period = "1950-05-01:1950-05-03",
  output = "raster"
)

# Vykreslení jednotlivých vrstev (dnů) a přidání hranic polygonu nps_region
ggplot() +
  geom_spatraster(data = nps_temp_1950) +
  facet_wrap(~lyr, ncol = 3) +
  scale_fill_distiller(palette = "RdYlBu", na.value = "transparent") +
  geom_spatvector(data = nps_region, fill = NA) +
  labs(fill = "Maximum\ntemperature \n(ºC)") +
  scale_x_continuous(breaks = seq(13.3, 13.9, by = 0.3)) +
  scale_y_continuous(breaks = seq(48.8, 49.2, by = 0.2)) +
  theme_minimal()
Na třech mapách je vidět maximální teploty pro jednotlivé dny 1.-3.5.1950. Legenda má jiný rozsah hodnot.

Na třech mapách je vidět maximální teploty pro jednotlivé dny 1.-3.5.1950. Legenda má jiný rozsah hodnot.

Nyní můžeme spočítat rozdíl pro jednotlivé dny mezi roky 2020 a 1950 a opět zobrazit v mapách. Heterogenita denních teplot je vysoká, v analýzách je tak většinou lepší spočítat agregovat data na za delší časová období - měsíce, sezóny nebo roky a ty porovnávat mezi sebou.

# Vypočet rozdílu mezi teplotami (nps_temp - nps_temp_1950)
nps_temp_diff <- nps_temp - nps_temp_1950

ggplot() +
  geom_spatraster(data = nps_temp_diff) +
  facet_wrap(~lyr, ncol = 3) +
  scale_fill_distiller(palette = "RdYlBu", na.value = "transparent") +
  geom_spatvector(data = nps_region, fill = NA) +
  labs(fill = "Rozdíl \nmaximální \nteploty \n2020-1950 \n(ºC)") +
  scale_x_continuous(breaks = seq(13.3, 13.9, by = 0.3)) +
  scale_y_continuous(breaks = seq(48.8, 49.2, by = 0.2)) +
  theme_minimal()
V jednotlivých dnech byly rozdíly mezi rastry mezi lety 2020 a 1950 velmi odlišné. Při srovnání jednotlivých dnů jde především o specifické počasí a těžko lze usuzovat na vliv změny klimatu.

V jednotlivých dnech byly rozdíly mezi rastry mezi lety 2020 a 1950 velmi odlišné. Při srovnání jednotlivých dnů jde především o specifické počasí a těžko lze usuzovat na vliv změny klimatu.

Interaktivní mapový výstup

Následující ukázka představuje postup, jak vizualizovat mapu v interaktivním prostředí s podkladovou mapou. To se může hodit, pokud chcete výsledky ukázat někomu, kdo nepoužívá R ani GIS, ale chce data prozkoumat více, než jen z obrázků.

if (is_html_output || is_interactive) { # Kód pro HTML verzi nebo interaktivní režim
  nps_temp_diff_1 <- nps_temp_diff[[1]]  # Vyberte vrstvu dle potřeby
  nps_temp_diff_1_raster <- as(nps_temp_diff_1, "Raster")
  raster_values <- raster::values(nps_temp_diff_1_raster)
  
  # Barevná paleta pro vykreslení
  pal <- colorNumeric(palette = "RdYlBu", domain = raster_values, 
                      na.color = "transparent")
  
  # Vytvoření interaktivní mapy
  leaflet() %>%
    addTiles() %>%
    addRasterImage(nps_temp_diff_1_raster, colors = pal, opacity = 0.8) %>%
    addPolygons(data = as(nps_region, "Spatial"), fill = NA, color = "black") %>%
    addLegend(pal = pal, values = raster_values, title = "Rozdíl Tmax 2020-1950 (ºC)")
  
} else { # Kód pro ne-HTML verzi (PDF/Word - statický obrázek)
  nps_temp_diff_1 <- nps_temp_diff[[1]]  # Vyberte vrstvu dle potřeby
  nps_temp_diff_1_raster <- as(nps_temp_diff_1, "Raster")
  raster_values <- raster::values(nps_temp_diff_1_raster)
  
  # Barevná paleta
  pal <- colorNumeric(palette = "RdYlBu", domain = raster_values, 
                      na.color = "transparent")
  
  # Vytvoření interaktivní mapy, kterou uložíme jako HTML
  map <- leaflet() %>%
    addTiles() %>%
    addRasterImage(nps_temp_diff_1_raster, colors = pal, opacity = 0.8) %>%
    addPolygons(data = as(nps_region, "Spatial"), fill = NA, color = "black") %>%
    addLegend(pal = pal, values = raster_values, title = "Rozdíl Tmax 2020-1950 (ºC)")
  
  # Uložení widgetu jako HTML souboru
  htmlwidgets::saveWidget(map, "./images/map.html", selfcontained = TRUE)
  
  # Vytvoření statického snímku mapy pomocí webshot
  webshot2::webshot("./images/map.html", file = "./images/map.png", cliprect = "viewport",
                    vwidth = 1000, vheight = 800, delay = 5)
  
  # Vložení obrázku do dokumentu
  knitr::include_graphics("./images/map.png", dpi=150)
}

Interaktivní mapa využívá online podkladovou mapu. Hlavní výhodou je však možnost snadné manipulace s mapou, což však není možné ve formátu PDF ukázat.

Pokud používáte knihovnu easyclimate, uveďte prosím citaci jak pro příslušný zdroj dat, tak pro knihovnu: Cruz-Alonso V, Pucher C, Ratcliffe S, Ruiz-Benito P, Astigarraga J, Neumann M, Hasenauer H, Rodríguez-Sánchez F (2023). “The easyclimate R package: Easy access to high-resolution daily climate data for Europe.” Environmental Modelling & Software, 105627. https://doi.org/10.1016/j.envsoft.2023.105627.

Pro další tutoriály navštivte články na webu knihovny: https://verughub.github.io/easyclimate/.

Data z rastrů E-OBS pro body - knihovna easyclimate

Spojité rastry lze též použít pro získání klimatických dat pro lokality. Je třeba mít na paměti, že rozlišení rastru je 1 km, takže se hodí spíše pro srovnání s lokálním klimatem. V následujícím příkladu stáhneme data o maximální teplotě a srovnáme s hodnotami na naší stanici.

Příprava knihoven prostředí R.

Zde definujeme body a stahujeme pro ně data pro studované období 1.5. - 30.6.2022.

# Definujte souřadnice lokalit (longitude a latitude) a období
locations <- data.frame(
  lon = c(13.5678442),  # příklad souřadnic, lze přidat další např. c(13.3, 13.8)
  lat = c(49.0679850)   # příklad souřadnic, lze přidat další např. c(49.0, 49.5)
)

start_date <- "2022-05-01"
end_date <- "2022-06-30"
period <- paste(start_date, end_date, sep = ":")

# Funkce pro stažení dat pro konkrétní lokalitu
get_climate_data <- function(lon, lat, period) {
  climate_data <- get_daily_climate(
    coords = data.frame(lon = lon, lat = lat),
    period = period,
    climatic_var = "Tmax",  # lze změnit na "Prcp" pro srážky nebo "Tmin" pro min. teploty
    # průměrná teplota zde není, ale běžně se počítá jako (Tmax+Tmin)/2
    version = 4  # standardní verze datasetu
  )
  return(climate_data)
}

# Stáhněte data pro všechny lokality
climate_data_list <- lapply(1:nrow(locations), function(i) {
  get_climate_data(locations$lon[i], locations$lat[i], period)
})

# Spojte data do jednoho dataframe
climate_data <- bind_rows(climate_data_list)

# Převeďte sloupec 'date' z typu character na Date
climate_data <- climate_data %>%
  mutate(date = as.Date(date, format = "%Y-%m-%d"))

rm(climate_data_list, locations, get_climate_data)

Nyní načteme naše data ze stanice NPS_4210_D_T1, ořízneme na stejné časové období a agregujeme na denní maximální teplotu.

## Načtení souborů stanic TOMST Thermologger pomocí files_table
ft <- read.table("./data/stanice/files_table.csv", sep=";", header = TRUE)
ft$path<-paste0("./data/stanice/",ft$path)
ft<-ft[1,] # jen jedna stanice

# Načtení dat loggerů s metadaty
tms <- mc_read_data(files_table = ft, silent = TRUE)

## Oříznutí časové řady na vybrané období
start <- as.POSIXct(start_date, tz = "UTC")
end <- as.POSIXct(end_date, tz = "UTC")
tms <- mc_prep_crop(tms, start, end, end_included=TRUE)

## Agregace všech časových řad, výpočet průměrů, rozptylů 
tms.agg <- mc_agg(tms, fun = c("max"), period = "day", min_coverage = 0.95)

## Převod na dlouhý formát (pro lepší manipulaci s daty)
tms.agg.long <- mc_reshape_long(tms.agg)

## Převod na široký formát pro snadnější analýzu
df_wide.tms <- tms.agg.long %>%
  dplyr::select(-height) %>%
  pivot_wider(names_from = sensor_name, values_from = value)
df_wide.tms$date<-as.Date(df_wide.tms$datetime)
df_selected <- df_wide.tms[, c("date", "Thermo_T_max")]

Nyní můžeme data z naší stanice srovnat s hodnotami z rastru E-OBS.

# Srovnání dat
comparison_data <- merge(climate_data, df_selected, 
                         by.x = "date", by.y = "date", all = TRUE)

# Vizualizace srovnání
ggplot(comparison_data, aes(x = date)) +
  geom_line(aes(y = Tmax, color = "E-OBS Data")) +
  geom_line(aes(y = Thermo_T_max, color = "NPS_4210_D_T1")) +
  labs(title = paste0("Období ", period),
       x = "Datum",
       y = "Maximální teplota (°C)") +
  theme_minimal()
Srovnání maximálních teplot pro období 1.5.2022 - 30.6.2022. Přestože je E-OBS založen na datech meteorologických stanic a má rozlišení 1 km, při práci s denními maximálními teplotami jsou data srovnatelná.

Srovnání maximálních teplot pro období 1.5.2022 - 30.6.2022. Přestože je E-OBS založen na datech meteorologických stanic a má rozlišení 1 km, při práci s denními maximálními teplotami jsou data srovnatelná.

Nyní můžeme stejné hodnoty srovnat pomocí scatterplotu a podívat se, zda není mezi hodnotami nějaký trend. Pokud by byly rozdíly náhodné, byl by trend rovnoběžný s linií 1:1.

# Zjistěte minimální a maximální hodnoty pro nastavení rozsahu obou os
min_val <- min(c(comparison_data$Tmax, comparison_data$Thermo_T_max), na.rm = TRUE)
max_val <- max(c(comparison_data$Tmax, comparison_data$Thermo_T_max), na.rm = TRUE)

# Vytvořte scatter plot s trendovou čarou, osami stejného rozsahu a čarou 1:1
ggplot(comparison_data, aes(x = Tmax, y = Thermo_T_max)) +
  geom_point(color = "blue", alpha = 0.6) +  # scatter plot bodů
  geom_abline(slope = 1, intercept = 0, color = "red", linetype = "dashed") +  # 1:1 linie
  geom_smooth(method = "lm", color = "green", se = FALSE) +  # regresní linie
  geom_hline(yintercept = 0, color = "black") +  # přidání horizontální osy 0
  geom_vline(xintercept = 0, color = "black") +  # přidání vertikální osy 0
  coord_equal(xlim = c(min_val, max_val), ylim = c(min_val, max_val)) + # stejné limity os
  labs(title = paste0("Období ", period),
       x = "E-OBS Tmax (°C)",
       y = "Tmax na NPS_4210_D_T1 (°C)") +
  theme_minimal()
Srovnání NPS_4210_D_T1 a E-OBS včetně 1:1 linie a regresní linie (červená tečkovaná). Regresní přímka je strmější než 1:1 linie, takže při nižších maximálních teplotách E-OBS jsou teploty na mikroklimatické stanici nižší, ale při vyšších teplotách jsou naopak vyšší, než uvádí E-OBS. Rozdíly mohou být dány vlivem porostu, pro přesnější analýzu by bylo třeba získat více dat a také zohlednit jejich závislost v čase.

Srovnání NPS_4210_D_T1 a E-OBS včetně 1:1 linie a regresní linie (červená tečkovaná). Regresní přímka je strmější než 1:1 linie, takže při nižších maximálních teplotách E-OBS jsou teploty na mikroklimatické stanici nižší, ale při vyšších teplotách jsou naopak vyšší, než uvádí E-OBS. Rozdíly mohou být dány vlivem porostu, pro přesnější analýzu by bylo třeba získat více dat a také zohlednit jejich závislost v čase.

4.3. Rastry mikroklimatu pro NP České Švýcarsko zpracované BÚ

Podobně jako pro NP Šumava vznikl mikroklimatický atlas i pro NP České Švýcarsko a NP Saské Švýcarsko. Pomocí následujícího kódu si opět můžete stáhnout vybranou mapu a vizualizovat ji pomocí mapy v prostředí R.

https://gitlab.ibot.cas.cz/matej.man/microclimate-atlas-public/

V tomto případě jsou dostupné vrstvy pro roky 2018-2020, rok 2018 odpovídá období listopad 2018 - říjen 2019. Odlišně od mikroklimatických rastrů NP Šumava byly nelesní oblasti v NP České Švýcarsko vyříznuty.

Dostupné jsou tyto proměnné:

  • 2018_FDD_baseT0_air_200cm_epsg32633.tif

  • 2018_GDD_baseT5_air_200cm_epsg32633.tif

  • 2018_T_max95p_air_15cm_epsg32633.tif

  • 2018_T_max95p_air_200cm_epsg32633.tif

  • 2018_T_mean_air_15cm_epsg32633.tif

  • 2018_T_mean_air_200cm_epsg32633.tif

  • 2018_T_mean_soil_8cm_epsg32633.tif

  • 2018_T_min5p_air_200cm_epsg32633.tif

  • 2018_snow_sum_epsg32633.tif

(pro další roky je v názvu pouze jiný začátek)

# URL ke stažení souboru (soubor T.air_200_cm.mean.tif z Zenodo)
soubor<-"2020_T_mean_air_15cm_epsg32633.tif"
url <- paste0("https://gitlab.ibot.cas.cz/",
              "matej.man/microclimate-atlas-public/-/",
              "raw/main/geodata/EPSG32633/",soubor)

# Cesta, kam uložit soubor
destfile <- paste0("./data/rastry/", soubor)

# Vytvoření složky, pokud neexistuje
dir.create("./data/rastry", recursive = TRUE, showWarnings = FALSE)

# Stažení souboru, pokud neexistuje
if (!file.exists(destfile)) {
  download.file(url, destfile, mode = "wb")
}

# Načtení rastru pomocí knihovny terra
npcs <- terra::rast(destfile)

Mapa

Opět se podíváme na mapu jedné vrstvy. Jakmile máme mapu načtenou jako rastr, jsou už všechny další kroky stejné, jako v kapitole 4.1.

# Nastavení prostoru pro legendu
par(mar = c(5, 4, 4, 5))  # Přidání prostoru na pravé straně pro legendu
par(oma = c(2, 1, 1, 1))  # Zvýšení vnějšího spodního okraje

# Základní graf rastru pomocí funkce plot()
plot(npcs, col = viridis::viridis(100), 
     xlab = "Zeměpisná délka", ylab = "Zeměpisná šířka")

# Přidání popisku k legendě s jednotkou (°C)
mtext("Teplota (°C)", side = 4, line = 3)  # Přidání popisku na pravou stranu
Jednoduchá mapa ukazuje průměrnou teplotu pro období listopad 2020 - říjen 2021.

Jednoduchá mapa ukazuje průměrnou teplotu pro období listopad 2020 - říjen 2021.

5. Závěr

Cílem metodiky bylo usnadnit využití mikroklimatických dat v praxi ochrany přírody a zajistit reprodukovatelnost výsledků. V rámci projektu a jeho implementace budeme metodiku dále udržovat, tak aby byla funkční a sloužila svému účelu. Pokud se najdou vhodná témata, plánujeme ji rozšířit o další příklady, především na základě požadavků z praxe ochrany přírody.

Literatura

BOX, G. E. P., and JENKINS, G. M. (1970) Time Series Analysis: Forecasting and Control. San Francisco: Holden-Day.

BRŮNA, J., KLINEROVÁ, T., KONOPOVÁ, Z., KALČÍK, V. a KIRSCHNER, J. (2023) Využití mikroklimatických dat v památkách zahradního umění. Zahradnictví. vol. 11.

BRŮNA, J., MACEK, M., MAN, M., HEDEROVÁ, L., KLINEROVÁ, T., MOUDRÝ, V., HEURICH, M., ČERVENKA, J., WILD, J. a KOPECKÝ, M. (2023) High-resolution microclimatic grids for the Bohemian Forest Ecosystem (1.0) [Data set]. Zenodo. doi: https://doi.org/10.5281/zenodo.6352641

BRŮNA J., WILD J., HEDEROVÁ L., KLINEROVÁ T., MACEK M. a URBANOVÁ R. (2021) Metodika měření mikroklimatu pomocí mikroklimatických stanic systému TMS. Botanický ústav AVĆR http://hdl.handle.net/11104/0317943

CHIANUCCI F. a MACEK M. (2023) hemispheR: an R package for fisheye canopy image analysis. Agricultural and Forest Meteorology. vol.  336: 109470. doi: https://doi.org/10.1016/j.agrformet.2023.109470

CRAWLEY, M. J. (2013). The R Book. John Wiley & Sons.

CRUZ-ALONZO, V., PUCHER, C., RARCLIFFE, S., RUIZ-BENITO, P., ASTIGARRAGA, J., NEUMANN, M., HASENAUER, H., RODRÍGUEZ-SÁNCHEZ, F. (2023) “The easyclimate R package: Easy access to high-resolution daily climate data for Europe.” Environmental Modelling & Software, 105627. doi: https://doi.org/10.1016/j.envsoft.2023.105627.

HAESEN, S., et al. (2021) ForestTemp – Sub-canopy microclimate temperatures of European forests. Global Change Biology. vol. 27(23), s. 6307–6319. doi: https://doi.org/10.1111/gcb.15892

HAESEN, S., et al. (2023) ForestClim – Bioclimatic variables for microclimate temperatures of European forests. Global Change Biology. vol. 29: s. 2886-2892 doi: https://doi.org/10.1111/gcb.16678

LEMBRECHTS J.J., et al. (2022) Global maps of soil temperature. Global Change Biology vol. 28: s. 3110–3144. doi: https://doi.org/10.1111/gcb.16060

KOPECKÝ, M., MACEK, M. a WILD, J. (2021) Topographic Wetness Index calculation guidelines based on measured soil moisture and plant species composition. Science of the Total Environment. vol. 757, 143785. doi: https://doi.org/10.1016/j.scitotenv.2020.143785

MACEK, M., KOPECKÝ, M., a WILD, J. (2019) Maximum air temperature controlled by landscape topography affects plant species composition in temperate forests. Landscape Ecology. vol. 34: s. 2541–2556. doi: https://doi.org/10.1007/s10980-019-00903-x

MAN, M., KALČÍK, V., MACEK, M., BRŮNA, J., HEDEROVÁ, L., WILD, J. a KOPECKÝ, M. (2023) myClim: Microclimate data handling and standardised analyses in R. Methods in Ecology and Evolution. vol. 14: s. 2308-2320. doi: https://doi.org/10.1111/2041-210X.14192

MORENO, A., HASENAUEER, H. (2016) Spatial downscaling of European climate data. International Journal of Climatology, 1444–1458. doi: https://doi.org/10.1002/joc.4436

NOAA National Centers of Environmental Information (1999) Global Surface Summary of the Day - GSOD. 1.0. NOAA National Centers for Environmental Information. id: gov.noaa.ncdc:C00516

PUCHER, C., NEUMANN, M. (2022). Description and Evaluation of Downscaled Daily Climate Data Version 3. figshare doi: https://doi.org/10.6084/m9.figshare.19070162.v1

PUCHER, C. (2023). Description and Evaluation of Downscaled Daily Climate Data Version 4. figshare doi: https://doi.org/10.6084/m9.figshare.22962671.v1

SPARKS, A. H., HENGL, T. a NELSON, A. (2017) GSODR: Global Summary Daily Weather Data in R. The Journal of Open Source Software, 2(10). doi: https://doi.org/10.21105/joss.00177

WILD, J., KIRSCHNER, J., MORAVEC, D. a KOHLOVÁ, J. (2014) Microclimate measurement as one of the prerequisites for succesful introduction of ornamental trees. Acta Pruhoniciana. vol. 108: s. 5–13.

WILD, J., KOPECKÝ, M., MACEK, M., ŠANDA, M., JANKOVEC, J. a T. HAASE. (2019) Climate at ecologically relevant scales: A new temperature and soil moisture logger for long-term microclimate measurement. Agricultural and Forest Meteorology. vol. 268: s. 40–47. doi: https://doi.org/10.1016/j.agrformet.2018.12.018