153ZODH / 4. cvičení
Mapová algebra, lokální operace, podíl obrazu, NDVI
Osnova
Úvod do rastrové algebry, lokální operace - podíl obrazu, normalizovaný diferenční vegetační index (NDVI).
Seznam použitých příkazů
Rastrová algebra
Rastrovou (mapovou) algebru (Map algebra) má v systému GRASS na starost modul r.mapcalc (více tutoriál).
Dané výrazy lze zadat jako parametr expr (nutné uzavřít do uvozovek), nebo je uložit do textového souboru (při opětovném použití) a použít parametr file.
# zadání výrazu jako parametr
#
r.mapcalc expr='p12 = tm1 / tm2'
r.mapcalc expr='p23 = tm2 / tm3'
#
# výrazy uložené v textovém souboru
#
r.mapcalc file=podily.txt
Poznámka: Modul r.mapcalc provádí výpočet podobně jako ostatní moduly pro zpracování rastrových dat v aktuální výpočetním regionu, tento fakt je nutné mít na paměti.
Ve wxGUI je mapový kalkulátor 
dostupný z menu
Raster -> Raster map calculator
Obecná forma výrazu:
nova_mapa = vyraz (mapa1, mapa2, ...)
kde vyraz může zahrnovat názvy rastrových vrstev, celočíselné či reálné konstanty a podporované funkce.
| operátor | popis | typ | priorita |
|---|---|---|---|
| ^ | exponent | aritmetický | 5 |
| % | modulo (zbytek po celočíselném dělení) | aritmetický | 4 |
| / | dělení | aritmetický | 4 |
| * | násobení | aritmetický | 4 |
| + | sčítání | aritmetický | 3 |
| - | odečítání | aritmetický | 3 |
| == | rovnost | logický | 2 |
| != | nerovnost | logický | 2 |
| > | větší než | logický | 2 |
| < | menší než | logický | 2 |
| >= | větší rovno než | logický | 2 |
| <= | menší rovno než | logický | 2 |
| && | AND (a) | logický | 1 |
| || | OR (nebo) | logický | 1 |
| # | rozděluje tabulku barev | aritmetický | - |
| funkce | význam | mat. zápis | typ |
|---|---|---|---|
| abs(x) | absolutní hodnota z x | |x| | |
| atan(x) | arkustangens z x (výsledek ve stupních) | arctan(x) | F |
| atan(x,y) | arkustangens z x/y (výsledek ve stupních) | arctan(x/y) | F |
| cos(x) | kosinus z x (x ve stupních) | cos(x) | F |
| double(x) | převede x na číslo s plovoucí desetinnou čárkou (výstupní mapa DCELL) | F | |
| eval([x,y,...,]z) | vyhodnotí výrazy x a y a výsledek uloží do z | ||
| exp(x) | exponenciální funkce z x | ex | F |
| exp(x,y) | mocnina y z x | xy | F |
| float(x) | převede x na číslo s plovoucí desetinnou čárkou (výstupní mapa FCELL) | F | |
| if | podmíněný příkaz | ||
| if(x) | 1, když x není 0, jinak 0 | ||
| if(x,a) | a, když x není 0, jinak 0 | ||
| if(x,a,b) | a, když x není 0, jinak b | ||
| if(x,a,b,c) | a, když x > 0, b když x rovno 0, c když x < 0 | ||
| int(x) | převede x na celé číslo, zbytek odřízne | ||
| isnull() | 1, když x rovno NULL (no data) | F | |
| log(x) | přirozený logaritmus z x | ln(x) | F |
| log(x,b) | logaritmus z x při základu b | logb(x) | F |
| max(x,y [,z...]) | vrátí největší ze zadaných hodnot | * | |
| median(x,y[,z...]) | vrátí prostřední hodnotu (medián) ze zadaných hodnot | * | |
| min(x,y[,z...]) | vrátí nejmenší ze zadaných hodnot | * | |
| mode(x,y[,z...]) | vrátí nejčetnější hodnotu (modus) ze zadaných hodnot | * | |
| not(x) | 1 pokud x je nula, jinak 0 | ||
| rand(low,high) | vrátí náhodnou číselnou hodnotu z intervalu [low; high] | ||
| round(x) | zaokrouhlí x na nejbližší celočíselnou hodnotu | I | |
| sin(x) | sinus z x (x ve stupních) | sin(x) | F |
| sqrt(x) | odmocnina z x | sqrt(x) | F |
| tan(x) | tangens z x (x ve stupních) | tan(x) | F |
Vysvětlivky:
| * | ... | Vrací reálné číslo, pokud je x reálné číslo. |
| F | ... | Výsledek je vždy reálné číslo. |
| I | ... | Výsledek je celé číslo. |
| proměnná | význam |
|---|---|
| x() | proměnná pro aktuální souřadnici x |
| y() | proměnná pro aktuální souřadnici y |
| col() | proměnná pro aktuální číslo sloupce |
| row() | proměnná pro aktuální číslo řádku |
| nsres() | proměnná pro aktuální rozlišení sever-jih |
| ewres() | proměnná pro aktuální rozlišení východ-západ |
V následující části se zaměříme na lokalní operace, kdy hodnota výsledné buňky je vypočtena na základě hodnot korespondujících hodnot vstupní vrstvy/vrstev. V dalších cvičeních budeme používat operace fokální, které se provádí v rámci množiny rastrových buněk pokryté tzv. moving window - plovoucím oknem.
Poznámky:
- Datový typ buněk výsledné vrstvy (celočíselná; hodnoty s plovoucí desetinnou čárkou) se odvozuje z datového typu buněk vstupních vrstev a datovém typu konstant uvedených ve výrazu. Výsledkem dělení dvou celočíselných rastrových vrstev je tedy opět vrstva s celočíselnými hodnotami buněk!
- Vstupují-li do výpočtu rastrové buňky s hodnotou NULL ("žádná data") je výsledkem opět hodnota NULL.
Podíl obrazu
Podíl obrazu je základní algebraickou metodou používanou pro detekci příznaků, zvýraznění obrazu, výpočet vegetačních indexů a řadu dalších úloh.
Poznámka: Datový typ buněk vstupních vrstvev je často celočíselný (např. jednotlivé kánaly LandSat-TM). Abychom získali rastrovou vrstvu obsahující hodnoty s plovoucí desetinnou čárkou (FCELL, DCELL) je potřeba alespoň jeden z operandů přetypovat, např. přenásobit konstantou "1.0" nebo použít vestavěnou funkci float().
Jako příklad můžeme vypočítat podíl sedmého a čtvrtého kanálu družicové scény LandSat-TM5:
# podíl obrazových kanálů sedm a čtyři (výsledná mapa DCELL)
#
r.mapcalc 'p74 = 1.0 * tm7 / tm4'
d.rast map=p74
NDVI
Praktickou úlohou podílu obrazu je výpočet normalizovaného diferenčního vegetačního indexu (NDVI).
# výpočet NDVI
#
r.mapcalc 'ndvi = 1.0 * (tm4 - tm3) / (tm4 + tm3)'
Dominuje-li DH v blízkém infračerveném oboru (tm4) nad hodnotou ve viditelném (červeném) spektru (tm3), hodnota NDVI je kladná a ukazuje na oblast s bohatou vegetací. Hodnoty NDVI kolem nuly reprezentují půdu bez výraznější vegetace, nula potom vodní plochy - pouze teoreticky.
Pro lepší vizuální interpretaci obrazu zvolíme vlastní tabulku barev, nejprve zjistíme rozsah dat (r.info -r), pro názornější představu můžeme vytisknout statistiku rastrové vrstvy pomocí modulu r.report:
# rozsah hodnot
#
r.info -r map=ndvi
min=-1.000000 max=0.861314
Rastrové hodnoty se tedy pohybují v intervalu [-1.0; +0.87].
# vytisknout statistiku (počet buněk, procenta, plocha v hektarech); přerozdělit na 25 intervalů
#
r.report -n map=ndvi units=c,p,h nstep=25
+-----------------------------------------------------------------------------+ | RASTER MAP CATEGORY REPORT | |LOCATION: sevcech Sat Oct 29 14:43:14 2005| |-----------------------------------------------------------------------------| | north: -957500 east: -763855.0602047 | |REGION south: -1007318 west: -830529 | | res: 29.99277544 res: 29.99277544 | |-----------------------------------------------------------------------------| |MASK:none | |-----------------------------------------------------------------------------| |MAP: (untitled) (ndvi in student) | |-----------------------------------------------------------------------------| | Category Information | cell| % | | | #|description | count| cover| hectares| |-----------------------------------------------------------------------------| | -1--0.925547|from to . . . . . . . . . . . | 37| 0.00| 3.33| |-0.553285--0.478832|from to . . . . . . . . . . . | 1| 0.00| 0.09| | -0.478832--0.40438|from to . . . . . . . . . . . | 5| 0.00| 0.45| | -0.40438--0.329927|from to . . . . . . . . . . . | 114| 0.00| 10.26| |-0.329927--0.255474|from to . . . . . . . . . . . | 2258| 0.06| 203.12| |-0.255474--0.181022|from to . . . . . . . . . . . | 13612| 0.38| 1224.49| |-0.181022--0.106569|from to . . . . . . . . . . . | 11660| 0.32| 1048.89| |-0.106569--0.032117|from to . . . . . . . . . . . | 70927| 1.96| 6380.36| | -0.032117-0.042336|from to . . . . . . . . . . . | 88807| 2.46| 7988.78| | 0.042336-0.116788|from to . . . . . . . . . . . | 171190| 4.74| 15,399.68| | 0.116788-0.191241|from to . . . . . . . . . . . | 240714| 6.67| 21,653.83| | 0.191241-0.265693|from to . . . . . . . . . . . | 265466| 7.35| 23,880.43| | 0.265693-0.340146|from to . . . . . . . . . . . | 342896| 9.50| 30,845.78| | 0.340146-0.414599|from to . . . . . . . . . . . | 594617| 16.47| 53,489.76| | 0.414599-0.489051|from to . . . . . . . . . . . | 822636| 22.79| 74,001.59| | 0.489051-0.563504|from to . . . . . . . . . . . | 752950| 20.86| 67,732.87| | 0.563504-0.637956|from to . . . . . . . . . . . | 223519| 6.19| 20,107.02| | 0.637956-0.712409|from to . . . . . . . . . . . | 8126| 0.23| 730.99| | 0.712409-0.786861|from to . . . . . . . . . . . | 1| 0.00| 0.09| | 0.786861-0.861314|from to . . . . . . . . . . . | 1| 0.00| 0.09| |-----------------------------------------------------------------------------| |TOTAL |3609537|100.00|324,701.88| +-----------------------------------------------------------------------------+
Tabulka barev:
-1.00 red 0.00 red 0.30 yellow 0.50 green 0.60 0 136 26 0.87 0 87 16
# nastavit tabulku barev ze souboru tb-ndvi.txt
#
r.colors map=ndvi color=rules rules=tb-ndvi.txt
#
# ... a překreslit obsah grafického okna
#
d.redraw
Pro zobrazení samostatné legendy slouží modul d.legend. Zobrazení legendy je ve wxGUI dostupné z nástrojové lišty mapového okna.
Současné zobrazení legendy a rastrové vrstvy nabízí (za cenu omezené možnosti ovlivnit podobu legendy) modul d.rast.leg.
Poznámka: Modul d.rast.leg není dostupný pro wxGUI. Následující ukázka nebude v prostředí wxGUI fungovat.
# současné zobrazení rastrových dat a odpovídající legendy
#
d.rast.leg map=ndvi
Flexibilnější podoba legendy vyžaduje více příkazů:
# rozdělit grafické okno na dva rámce
#
d.frame frame=f1 at=0,100,0,15;d.frame frame=f2 at=0,100,15,100
#
# v jednom zobrazit ndvi a popisek
#
echo "NDVI" | d.text size=5 color=black;d.rast map=ndvi
#
# a v druhém legendu (pouze pro rozsah [0;1])
#
d.frame -s frame=f1;d.legend map=ndvi range=0,1 labelnum=7
Reklasifikace rastrových dat
Na základě NDVI provedeme reklasifikaci na tři základní třídy:
- bez vegetace, vodní plochy
- plochy s minimální vegetací
- plochy pokryté vegetací
Modul určený pro reklasifikaci rastrových dat r.reclass nepodporuje data s plovoucí desetinnou čárkou, musíme si pomoci malým trikem - vytvoříme pomocnou vrstvu, která vznikne přenásobením NDVI konstantou "100". Zachováme si tak dostačující přesnost.
# pomocná vrstva pro r.reclass
#
r.mapcalc 'temp1 = 100 * ndvi'
Pro kontrolu vytiskneme rozsah hodnot rastrové vrstvy 'temp1':
r.info -r map=temp1
min=-100.000000 max=86.131387
Připravíme si tzv. reklasifikační tabulku:
-100 thru 5 = 1 bez vegetace, vodni plochy 5 thru 35 = 2 plochy s minimalni vegetaci 35 thru 87 = 3 plochy pokryte vegetaci
Jednotlivé řádky obsahují reklasifikační pravidla, např. první reklasifikační pravidlo přiřadí hodnotám rastrových buněk -100 až 5 výstupní hodnotu 1 s textovým popiskem "bez vegetace, vodni plochy". Reklasifikovaná rastrová vrstva tak bude obsahovat celkem tři kategorie 1, 2, 3.
# reklasifikace rastrových dat (reklasifikační pravidla uložena v textovém souboru)
#
r.reclass input=temp1 output=r_ndvi rules=reklas-ndvi.txt
Dále nastavíme vhodnou tabulku barev
1 red 2 yellow 3 0 136 26
# nastavení tabulky barev z textového souboru
#
r.colors map=r_ndvi color=rules rules=tb-r_ndvi.txt
#
# zobrazení reklasifikovaných dat
#
d.rast map=r_ndvi
d.rast.leg map=r_ndvi
Reklasifikační modul nevytváří ve skutečnosti plnohodnotnou rastrovou vrstvu, pouze zapíše zvolená reklasifikační pravidla.
Souhrná statistika reklasifikované rastrové vrstvy:
r.report map=r_ndvi units=c,p,h
+-----------------------------------------------------------------------------+ | RASTER MAP CATEGORY REPORT | |LOCATION: sevcech Mon Oct 31 01:15:33 2005| |-----------------------------------------------------------------------------| | north: -957500 east: -763855.0602047 | |REGION south: -1007318 west: -830529 | | res: 29.99277544 res: 29.99277544 | |-----------------------------------------------------------------------------| |MASK:none | |-----------------------------------------------------------------------------| |MAP: Reclass of temp1 in student (r_ndvi in student) | |-----------------------------------------------------------------------------| | Category Information | cell| % | | |#|description | count| cover| hectares| |-----------------------------------------------------------------------------| |1|bez vegetace, vodni plochy . . . . . . . . . . . | 146375| 3.96| 13,167.41| |2|plochy s minimalni vegetaci. . . . . . . . . . . | 820825| 22.23| 73,838.67| |3|plochy pokryte vegetaci. . . . . . . . . . . . . |2642337| 71.56|237,695.81| |*|no data. . . . . . . . . . . . . . . . . . . . . | 82866| 2.24| 7454.35| |-----------------------------------------------------------------------------| |TOTAL |3692403|100.00|332,156.23| +-----------------------------------------------------------------------------+
Praktická úloha: vytvoření rastrové mapy s teplotou povrchu
Nejprve vytvoříme kopii šestého kanálu LandSat-TM5 (pomocí modulu g.copy):
# kopie tm6@PERMANENT -> tm6_null@<mapset>
#
g.copy rast=tm6,tm6_null
Nyní nahradíme hodnotu "0" hodnotou NULL (žádná data) - nulové hodnoty obsahuje levý dolní roh snímku, který neobsahuje žádná data. Nahradíme-li tuto hodnotu hodnotou NULL, nebude tato část snímku vstupovat do výpočtu. Nahrazení hodnotou NULL provedeme příkazem r.null:
g.region rast=tm6
#
# nahrazení hodnoty "0" hodnotou NULL
#
r.null map=tm6_null setnull=0
Výpočet převzat z Open source GIS: A GRASS GIS Approach.
r.mapcalc 'tm6rad = 1.238 + (15.6 - 1.238) * tm6_null / 255'
r.mapcalc 'temp_kelvin = 1260.56 / (log (607.76 / tm6rad + 1))'
r.mapcalc 'temp_celsius = temp_kelvin - 273.15'
Před vizualizací výsledné rastrové vrstvy nastavíme vhodnou tabulku barev:
0 blue 15 green 30 yellow 45 red
r.colors map=temp_celsius col=rules < tb-teplota.txt
d.rast temp_celsius
Na závěr můžeme vzniknuvší rastrové vrstvy odstranit pomocí modulu g.remove.
# odstranění zvolených rastrových vrstev
#
g.remove rast=tm6_null,tm6rad,temp_kelvin