153ZODH / 4. cvičení: Porovnání verzí

← Porovnání se starší verzí Porovnání s novější verzí →

VizuálníWikitext

Verze z 2. 10. 2011, 21:42

< Stránky předmětu • Předchozí cvičení • Další 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 GRASSu na starost modul r.mapcalc (více tutoriál).

Modul r.mapcalc lze ovládat ve dvou módech - v interaktivním, kdy spustíte modul bez parametrů a zadáváte jednotlivé výrazy, pro ukončení slouží klíčové slovo "end". Poznámka: V GRASS 7 podporuje r.mapcalc pouze neinteraktivní mód.

 # interaktivní mód modulu r.mapcalc
 #
 r.mapcalc

Enter expressions, "end" when done.
mapcalc> p12 = tm1 / tm2
mapcalc> p23 = tm2 / tm3
mapcalc> end

Výhodnější je však předat dané výrazy modulu jako parametr (nutné uzavřít do uvozovek), nebo je uložit do textového souboru (při opětovném použití) a přesměrovat do modulu přes standardní vstup.

 # zadání výrazu jako parametr
 #
 r.mapcalc 'p12 = tm1 / tm2'
 r.mapcalc 'p23 = tm2 / tm3'
 #
 # výrazy uložené v textovém souboru
 #
 r.mapcalc < 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,

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.

Tab. č.1: Operátory
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ý	-

Tab č.2: Funkce
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	`e^x`	F
`exp(x,y)`	mocnina y z x	`x^y`	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	`log_b(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.

Tab č.3: Speciální proměnné
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.

Obr č.2: (a) tm1 + tm2; (b) tm1 - tm2; (c) tm1 * tm2; (d) tm1 / tm2; (e) 1.0 * tm1 / tm2

Několik drobných poznámek:

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

Obr č.3: Podíl obrazu (tm7/tm4) - tabulka barev rainbow (výchozí)

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, současné zobrazení legendy a rastrové vrstvy nabízí (za cenu omezené možnosti ovlivnit podobu legendy) modul d.rast.leg.

 # 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

Obr č.9: Reklasifikovaný NDVI s legendou

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

< Stránky předmětu • Předchozí cvičení • Další cvičení

@@ Řádek 52: / Řádek 52: @@
 ''Poznámka:'' Modul {{GrassPrikaz|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,
+Ve wxGUI je mapový kalkulátor dostupný z menu
+ Raster -> Raster map calculator
 [[Image:grass-mapcalc.png|thumb|640px|center|Grafické rozhraní mapového kalkulátoru]]