Přeskočit obsah

07. Zpracování časoprostorových dat v GRASS GIS

GRASS GIS nabízí specializované nástroje pro zpracování časoprostorových dat, viz GRASS dokumentace.

GRASS zavádí tři datové typy, které jsou určeny pro zpracování časovéprostorové řady dat:

  • časoprostorové rastrové datové sady (strds)
  • časoprostorové 3D rastrové datové sady (str3ds)
  • časoprostorové vektorové datové sady (stvds)

Úvod (simulace povodně)

V této části navážeme na simulaci v povodně. Simulaci povodně tentokrát provedeme ve více časových krocím pomocí nástroje r.lake.series.

Note

Nástroj r.lake.series není výchozí součástí instalace GRASS GIS, jde o rozšíření - tzv. GRASS Addons. Tato rozšíření lze snadno nainstalovat pomocí g.extension nebo z hlavní nabídky Settings > Addons extensions > Install extension from addons.

Nástroj r.lake.series vytváří pro každý časový okamžik rastrovou vrstvu s výškou aktuální hladiny a tyto vrstvy přidá do výstupního časoprostorového datasetu.

Před výpočtem vytvoříme nový mapset zaplava a nastavíme výpočetní region:

g.region vector=praha_7 align=DTM grow=1000

Vytvořme časoprostorovou řadu s vodními hladinami mezi 180.5m a 185.5m s krokem 0,1m:

r.lake.series elevation=DTM out=zaplava start_water_level=180.5 end_water_level=185.5 \
water_level_step=0.1 coordinates=-741903.25,-1040411.12 \
time_step=10 nproc=3

Časový krok je 10 minut. Pro zvýšení rychlosti výpočtu použijme více jader procesoru (nproc). Výsledek je uložen v časoprostorové datové sadě s názvem zaplava.

Základní informace o časoprostorové datové sadě získáme pomocí nástroje t.info:

t.info input=zaplava

...
+-------------------- Relative time -----------------------------------------+
 | Start time:................. 1
 | End time:................... 501
 | Relative time unit:......... minutes
 | Granularity:................ 10
 | Temporal type of maps:...... point
...

Informace o časové topologii vypišme pomocí t.topology:

t.topology input=zaplava

 ...
 +-------------------- Temporal topology -------------------------------------+
 ...
 | Number of points: .......... 51
 | Number of gaps: ............ 50
 | Granularity: ............... 10
 ...

Dotazování

Pomocí t.rast.list lze vypsat rastrové vrstvy splňující danou podmínku. Vypišme rastrové vrstvy v rámci první hodiny simulované povodně:

t.rast.list input=zaplava order=start_time where="start_time <= 60"

name|mapset|start_time|end_time
zaplava_180.5|zaplava|1|None
zaplava_180.6|zaplava|11|None
zaplava_180.7|zaplava|21|None
...

Statistiku pro jednotlivé rastrové vsrtvy lze vypsat pomocí nástroje t.rast.univar. Vypočítejme statistiku pouze pro první hodinu povodně:

 t.rast.univar input=zaplava where="start_time < 60"

id|semantic_label|start|end|mean|min|max|...
zaplava_180.5@zaplava||1|None|1.5259037543889|1.52587890625e-05|6.46075439453125|...
zaplava_180.6@zaplava||11|None|1.50206839513946|1.52587890625e-05|6.56076049804688|...
zaplava_180.7@zaplava||21|None|1.55157986180243|1.52587890625e-05|6.66075134277344|...

Agregace

Agregaci dat lze provést pomocí t.rast.aggregate. Agregujme data záplavy po jedné hodině:

t.rast.aggregate input=zaplava output=zaplava_h basename=ag granularity=60 nproc=3

Příkaz vygeneruje novou časoprostorovou datovou sadu, kterou lze použít pro následné analýzy, jako je například jednorozměrná statistika:

t.rast.univar input=zaplava_h

id|semantic_label|start|end|mean|min|max|...
ag_00001@zaplava|61|1|None|1.43664612331051|1.52587890625e-05|6.71075439453125|...
ag_00002@zaplava|121|61|None|1.71285097692945|1.52587890625e-05|7.31075541178385|...
ag_00003@zaplava|181|121|None|2.05479336695301|1.52587890625e-05|7.91075388590495|...

Agregaci časoprostorových dat do jedné rastrové vrstvy lze provést pomocí t.rast.series. Agregujme data v rámci první hodiny povodně s cílem vypočítat maximální výšku hladiny:

t.rast.series input=zaplava output=zaplava_1_max method=maximum where="start_time < 60"

Statistiku vypišme pomocí r.univar:

r.univar map=zaplava_1_max

...
minimum: 1.52588e-05
maximum: 6.96075
mean: 1.65698
...

Ukázka maximální výšky povodně za první hodinu (na pozadí ČÚZK WMS Ortofoto):

Vizualizace

Vizualizaci časoprostorových dat ve formě animace umožňuje nástroj g.gui.timeline (File -> Animation tool). Vytvořme animaci na základě časoprostorového datasetu:

g.gui.animation strds=zaplava

Ukázka animace:

Tip

Podkladovou vrstvu ortofota lze z kontextového menu nad vrtsvu v záložce Layers stáhnout (Save web service layer) a poté použít v animačním nástroji jako běžnou rastrovou vrstvu.

Analýza teploty povrchu zemského povrchu

K výpočtu teploty povrchu zemského povrchu použijeme data MODIS (Moderate Resolution Imaging Spectroradiometer). MODIS je 36kanálový sensor se spekrálním rozsahem od viditelného po tepelně infračervené záření. Byl vypuštěn jako součást družice Terra v roce 1999 a dále družice Aqua v květnu 2002. Družice Terra (produkty MOD) prolétá dvakrát denně (v 10:30 a 22:30 místního času), družice Aqua (produkty MYD) prolétá také dvakrát denně (v 01:30 a 13:30 místního času). Zdroj: GRASS Wiki

Naše zájmová oblast (Česká republika) je pokryta čtyřmi dlaždicemi MODIS (viz. MODLAND grid):

  • h18v03
  • h18v04
  • h19v03
  • h19v04

Data MODIS jsou poskytována ve třech projekcích (sinusoidální, Lambertova azimutální Equal-Area a geografické). Pro naše účely data transformujeme do ETRS89 / LAEA Europe (EPSG 3035).

Vytvořme novou GRASS lokaci fgis-modis pomocí kódu EPSG (Select CRS from a list by EPSG or description).

Vstupme do nové lokace (mapset PERMANENT) a nainstalujme rozšíření (AddOn) i.modis pro stažení a import dat MODIS. Kromě tohoto rozšíření je nutné také doinstalovat další závislosti jako je např. knihovna pyMODIS. V záložce Console spusťme následující příkazy:

python3 -m pip install setuptools slugify pymodis
g.extension extension=i.modis

Rozšíření pro MODIS se skládá ze dvou nástrojů:

  • i.modis.download a
  • i.modis.import

Stažení dat

Poznámka

Předstažená data najdete na GIS.labu v adresáři Repository/155FGIS/07. Data překopírujte do domovského adresáře (nástroj i.modis.import vyžaduje právo zápisu do adresáře se vstupními daty).

Stáhněme si požadované dlaždice pro rok 2023 pomocí příkazu i.modis.download. Dále zvolme produkt (teplota zemského povrchu agregovaná po osmi dnech s rozlišením 1 km Terra/Aqua) :

i.modis.download settings=settings.txt folder=modis-2023 tiles=h18v03,h18v04,h19v03,h19v04 \
 product=lst_aqua_eight_1000,lst_terra_eight_1000 \
 startday=2023-01-01 endday=2023-12-31

Poznámka

Výstupní složka (modis-2023 v tomto případě) musí existovat, jinak nástroj selže.

Soubor settings.txt obsahuje dva řádky: jméno uživatele a heslo pro přístup ke službě stahování MODIS. Přečtěte si dokumentaci pyModis jak se zaregistrovat a nastavit svůj účet.

Import dat

Vstupní data MODIS lze importovat a transformovat do cílového souřadnicového systému pomocí i.modis.import. Naimportujme postupně data z družice Terra a poté Aqua:

i.modis.import -mwq files=modis-2023/listfileMOD11A2.061.txt \
outfile=tlist-mod.txt

i.modis.import -mwq files=modis-2023/listfileMYD11A2.061.txt \
outfile=tlist-myd.txt

Přepínač -m má za následek vytvoření mozaiky ze vstupních dlaždic. Vynecháme také z důvodu rychlosti importu QA (Quality Assestment) vrstvy (-q). Poslední uvedný přepínač (-w) má za následek zápis souboru outfile, který následně použijeme pro registraci vrstev do časoprostorového datasetu.

Ukázka mozaiky vytvořené ze vstupních dlaždic:

Volba zájmového území

Tip

Pro stažení administrativní hranice ČR lze použít zásuvný modul QGIS RÚIAN s volbou "Výšší územně správní celky (VUSC)".

Vytvořený soubor ve formátu GPKG lze poté importovat do GRASS lokace.

Výpočet teploty povrchu země

Nejprve nastavme pomocí r.mask masku pro výpočet. Nezapomeňme, že před vytvořením masky musí být nastavena výpočetní oblast. Ta bude definována vektorovou vrstvou státní hranice a zarovnána podle vstupní MODIS vrstvy.

g.region vector=staty align=MOD11A2.A2023001_mosaic_LST_Day_1km
r.mask vector=staty

Zkontrolujeme hodnoty rozsahu importovaných dat MODIS pomocí nástroje r.info:

r.info -r map=MOD11A2.A2023001_mosaic_LST_Day_1km

min=0
max=14532

Aby bylo možné určit LST (Land Surface Temperature) ze vstupních dat, musíme převést digitální hodnoty (DN) na stupnici Celsia nebo Kelvina:

C = DN * 0,02 - 273,15

Převod do Celsiovy stupnice lze provést pomocí r.mapcalc. Vhodné je také nahradit nulové hodnoty pomocí hodnotou žádná data (v terminologii GRASS hodnota NULL).

r.mapcalc expression="MOD11A2.A2023001_mosaic_LST_Day_1km_c = \
 if(MOD11A2.A2023001_mosaic_LST_Day_1km != 0, \
 MOD11A2.A2023001_mosaic_LST_Day_1km * 0.02 - 273.15, null())"

Zkontrolujme hodnoty rozsahu nové datové vrstvy LST.

r.info -r map=MOD11A2.A2023001_mosaic_LST_Day_1km_c

min=-18.55
max=13.47

Rekonstrukce LST ve stupnici Celsia (tabulka barev celsius):

Vytvoření časoprostorového datasetu

Nový časoprostorový dataset lze vytvořit pomocí nástroje t.create. Vytvořme dataset s názvem modis:

t.create output=modis title="MODIS LST 2023" desc="MODIS LST 2023"

V dalším kroku se importovaná data MODIS zaregistrují do časoprostorového datasetu pomocí t.register. Příkaz je třeba spustit dvakrát, jednou pro data z družice Terra a podruhé pro družici Aqua:

t.register input=modis file=tlist-mod.txt
t.register input=modis file=tlist-myd.txt

Zkontrolujme základní metadata o vytvořené datové sadě pomocí t.info:

t.info input=modis

...
| Start time:................. 2023-01-01 00:00:00
| End time:................... 2024-01-01 00:00:00
...
| Number of registered maps:.. 184

Úkol

Zkontrolujte granularitu. Produkt je agregován z hodnot v osmi denních intervalech.

MOD11A2.A2023001_mosaic_LST_Day_1km|PERMANENT|2023-01-01 00:00:00|2023-01-09 00:00:00

Digitální hodnoty (DN) je třeba převést na stupnici Celsia. Výraz mapové algebry lze na časoprostorovém datasetu aplikovat pomocí nástroje t.rast.mapcalc. Tento nástroj aplikuje výraz mapové algebry na všechny rastrové vrstvy registrované ve vstupní časoprostorové datové sadě.

Tip

Mnoho nástrojů pro zpracování časoprostorových dat (t.*) podporuje paralelizaci výpočtu (viz volba nproc).

t.rast.mapcalc input=modis output=modis_c nproc=3 basename=c \
 expression="if(modis != 0, modis * 0.02 - 273.15, null())"

Příkaz vytvoří novou časoprostorovou rastrovou datovou sadu a do ní registruje vytvořené rastrové vrstvy s přepočítanými digitálními hodnotami. Název nově vytvořených rastrových vrstev definuje parametr basename.

Poznámka

Všimněte si, že nové rastrové mapy budou vytvořeny v aktuální výpočetní oblasti a aplikovanou maskou.

Nástroj t.rast.colors umožňuje nastavit tabulku barev pro všechny rastrové vrstvy v časoprostorové datové sadě. Nástroj odpovídá funkcionalitou r.colors.

t.rast.colors input=modis_c color=celsius

Úkol

Spočtěte statistiku pomocí nástroje pro vrstvy mezi 1.4.2023 a 30.4.2023:

t.rast.univar modis_c where="start_time >= '2023-04-01' and start_time <= '2023-04-30'"

id|semantic_label|start|end|mean|min|max|...
c_049@spectral||2023-04-07 00:00:00|2023-04-15 00:00:00|12.684121214772|-11.3099999999999|22.99|...
c_050@spectral||2023-04-07 00:00:00|2023-04-15 00:00:00|0.757510639154788|-17.63|8.21000000000004|...
c_051@spectral||2023-04-07 00:00:00|2023-04-15 00:00:00|13.698337137793|-16.3699999999999|23.77|...

Agregace po měsících

Agregaci dat na základě zadané granularity provádí nástroj t.rast.aggregate. Vytvořme novou časoprostorovou řadu s agregovanými hodnotami LST po měsících:

t.rast.aggregate input=modis_c output=modis_cm basename=cm granularity="1 months" nproc=3

Spočtěme statistiku pomocí nástroje t.rast.univar:

t.rast.univar modis_cm

id|semantic_label|start|end|mean|min|max|...
cm_2023_01@spectral||2023-01-01 00:00:00|2023-02-01 00:00:00|1.44984784545898|-7.7133333333333|6.75333333333337|...
cm_2023_02@spectral||2023-02-01 00:00:00|2023-03-01 00:00:00|-1.81418308716353|-10.1566666666666|2.8266666666667|...
cm_2023_03@spectral||2023-03-01 00:00:00|2023-04-01 00:00:00|2.47208359286264|-10.3966666666666|10.2736363636364|...

Úkol

Spočtěte statistiku pro měsíce červen a červenec.

Rozložení dat v časovém ose umožňuje vizualizovat nástroj g.gui.timeline:

g.gui.timeline inputs=modis_cm

Úkol

Vytvořte animaci pomocí nástroje g.gui.animation:

Rozsah tabulky barev nastavte pomocí volby range. Minimální a maximální teplotu určíme pomocí t.info.

Výpočet průměrné teploty pro roční období

Data z časoprostorového datasetu mohou být extrahována pomocí t.rast.extract. Vytvořme nový dataset pokrývající vybrané roční období:

t.rast.extract input=modis_c where="start_time >= '2023-03-20' and start_time < '2023-06-20'" \
 output=modis_c_spring

Nástroj t.rast.series provádí agregregaci dat časoprostorového datasetu do formy rastrové vrstvy. Vypočtěme průměrné teploty pro vybrané roční období:

t.rast.series input=modis_c_spring output=modis_c_spring_avg method=average

Tip

Nástroj t.rast.series má parametr where, který umožňuje aplikovat výběr bez nutnosti vytvářet nový dataset pomocí t.rast.extract.

Průměrné teploty pro jaro 2023 vypočítejme pomocí r.univar:

minimum: 3.6155
maximum: 17.9013
mean: 11.744

Extrakce denních a nočních teplot

Vytvořme dva nové časoprostorové datasety pro denní teploty a zvlášť pro noční teploty, a to pomocí t.rast.extract:

t.rast.extract input=modis where="name LIKE '%Day%'" output=modis_day
t.rast.extract input=modis where="name LIKE '%Night%'" output=modis_night

t.rast.mapcalc input=modis_day output=modis_c_day nproc=3 basename=c \
 expression="if(modis_day != 0, modis_day * 0.02 - 273.15, null())"
t.rast.mapcalc input=modis_night output=modis_c_night nproc=3 basename=c \
 expression="if(modis_night != 0, modis_night * 0.02 - 273.15, null())"

Porovnejme teploty pro dané období v noci a ve dne:

t.rast.univar modis_c_day where="start_time >= '2023-04-01' and start_time < '2023-04-08'"

id|semantic_label|start|end|mean|min|max|...
c_25@spectral||2023-04-07 00:00:00|2023-04-15 00:00:00|12.684121214772|-11.3099999999999|22.99|...
c_26@spectral||2023-04-07 00:00:00|2023-04-15 00:00:00|13.698337137793|-16.3699999999999|23.77|...

t.rast.univar modis_c_night where="start_time >= '2023-04-01' and start_time < '2023-04-08'"

id|semantic_label|start|end|mean|min|max|...
c_25@spectral||2023-04-07 00:00:00|2023-04-15 00:00:00|0.757510639154788|-17.63|8.21000000000004|...
c_26@spectral||2023-04-07 00:00:00|2023-04-15 00:00:00|-0.465633774755256|-14.79|6.95000000000005|...

Ukázka modelu

Následující model provede import MODIS dat a přepočet digitálních hodnot do stupňů Celsia:

Model ke stažení zde.

Další materiály

Automatizaci zpracování časoprostorových dat pomocí specializovaného Python API najdete v materiálech GRASS GIS Jena Workshop.