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

Z GeoWikiCZ
VizuálníWikitext
m (sablona)
m (155zddp)
 
(Není zobrazeno 102 mezilehlých verzí od stejného uživatele.)
Řádek 1: Řádek 1:
<div style="font-size: 13pt; font-weight: bold">Mapová algebra, sčítání, odčítání rastrových dat, podíl obrazu, NDVI</div>
{{Zastaralé|155ZDDP}}
----
{{Cvičení|153ZODH|4|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ů ===
 
* {{GrassPrikaz|r.mapcalc}}
* {{GrassPrikaz|r.info}}
* {{GrassPrikaz|r.report}}
* {{GrassPrikaz|r.colors}}
* {{GrassPrikaz|r.reclass}}
* {{GrassPrikaz|g.copy}}
* {{GrassPrikaz|r.null}}
* {{GrassPrikaz|g.remove}}
 
== Rastrová algebra ==
 
'''Rastrovou (mapovou) algebru''' ({{wikipedia|Map algebra|lang=en}}) má v systému GRASS na starost modul {{GrassPrikaz|r.mapcalc}} (více [http://grass.osgeo.org/gdp/raster/mapcalc-algebra.pdf tutoriál]).
<!--
Modul {{GrassPrikaz|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:'' Interaktivní mód je dostupný pouze z příkazové řádky, ve {{grassPrikaz|wxGUI}} není podporován podporován.
 
<source lang="bash">
# interaktivní mód modulu r.mapcalc
#
r.mapcalc
</source>
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.
-->
 
Dané výrazy lze zadat jako parametr <code>expr</code> (nutné uzavřít do uvozovek), nebo je uložit do textového souboru (při opětovném použití) - viz parametr <code>file</code>.
 
<source lang="bash">
# 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
</source>
 
''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, tento fakt je nutné mít na paměti.
 
Ve {{GrassPrikaz|wxGUI}} je mapový kalkulátor [[Image:Calculator.png]] dostupný z menu {{menu|Raster|Raster map calculator}}
 
{{fig|grass-mapcalc|Grafické rozhraní mapového kalkulátoru|size=640}}
 
Obecná forma výrazu:
 
nova_mapa = vyraz (mapa1, mapa2, ...)
 
kde <tt>vyraz</tt> může zahrnovat názvy rastrových map, celočíselné či reálné konstanty a podporované funkce.
 
<center>
{|class=wikitable
|+Tab. č.1: Operátory
|-
!operátor !! popis !! typ !! priorita
|-
|<tt>^</tt>
|exponent
|aritmetický
|align="right" | 5
|-
|<tt>%</tt>
|modulo (zbytek po celočíselném dělení)
|aritmetický
|align="right" | 4
|-
|<tt>/</tt>
|dělení
|aritmetický
|align="right" | 4
|-
|<tt>*</tt>
|násobení
|aritmetický
|align="right" | 4
|-
|<tt>+</tt>
|sčítání
|aritmetický
|align="right" | 3
|-
|<tt>-</tt>
|odečítání
|aritmetický
|align="right" | 3
|- 
|<tt>==</tt>
|rovnost
|logický
|align="right" | 2
|- 
|<tt>!=</tt>
|nerovnost
|logický
|align="right" | 2
|- 
|<tt>></tt>
|větší než
|logický
|align="right" | 2
|-     
|<tt><</tt>
|menší než
|logický
|align="right" | 2
|- 
|<tt>>=</tt>
|větší rovno než
|logický
|align="right" | 2
|- 
|<tt><=</tt>
|menší rovno než
|logický
|align="right" | 2
|-     
|<tt>&&</tt>
|AND (a)
|logický
|align="right" | 1
|-
|<tt><nowiki>||</nowiki></tt>
|OR  (nebo)
|logický
|align="right" | 1
|- 
|<tt>#</tt>
|rozděluje tabulku barev
|aritmetický
|align="right" | -
|-
|}
 
{|class=wikitable
|+Tab č.2: Funkce
|-
!funkce !! význam !! mat. zápis !! typ
|-
|<tt>abs(x)</tt>
|absolutní hodnota z ''x''
|<tt><nowiki>|x|</nowiki></tt>
|&nbsp;
|-
|<tt>atan(x)</tt>
|arkustangens z ''x'' (výsledek ve stupních)
|<tt>arctan(x)</tt>
|F
|-
|<tt>atan(x,y)</tt>
|arkustangens z ''x/y'' (výsledek ve stupních)
|<tt>arctan(x/y)</tt>
|F
|-
|<tt>cos(x)</tt>
| kosinus z ''x'' (''x'' ve stupních)
|<tt>cos(x)</tt>
|F
|-
|<tt>double(x)</tt>
|převede ''x'' na číslo s plovoucí desetinnou čárkou (výstupní mapa DCELL)
|
|F
|-
|<tt>eval([x,y,...,]z)</tt>
|vyhodnotí výrazy ''x'' a ''y'' a výsledek uloží do ''z''
|&nbsp;
|&nbsp;
|-
|<tt>exp(x)</tt>
|exponenciální funkce z ''x''
|<tt>e<sup>x</sup></tt>
|F
|-
|<tt>exp(x,y)</tt>
| mocnina ''y'' z ''x''
|<tt>x<sup>y</sup></tt>
|F
|-
|<tt>float(x)</tt>
| převede ''x'' na číslo s plovoucí desetinnou čárkou (výstupní mapa FCELL)
|&nbsp;
|F
|-
|<tt>if</tt>
|podmíněný příkaz
|&nbsp;
|&nbsp;
|-
|<tt>if(x)</tt>
|''1'', když ''x'' není ''0'', jinak ''0''
|&nbsp;
|&nbsp;
|-
|<tt>if(x,a)</tt>
| ''a'', když ''x'' není ''0'', jinak ''0''
|&nbsp;
|&nbsp;
|-
|<tt>if(x,a,b)</tt>
|''a'', když ''x'' není ''0'', jinak ''b''
|&nbsp;
|&nbsp;
|-
|<tt>if(x,a,b,c)</tt>
|''a'', když ''x > 0'', ''b'' když ''x'' rovno ''0'', ''c'' když ''x < 0''
|&nbsp;
|&nbsp;
|-
|<tt>int(x)</tt>
|převede ''x'' na celé číslo, zbytek odřízne
|&nbsp;
|&nbsp;
|-
|<tt>isnull()</tt>
|''1'', když ''x'' rovno ''NULL'' (no data)
|&nbsp;
|F
|-
|<tt>log(x)</tt>
|přirozený logaritmus z ''x''
|<tt>ln(x)</tt>
|F
|-
|<tt>log(x,b)</tt>
|logaritmus z ''x'' při základu ''b''
|<tt>log<sub>b</sub>(x)</tt>
|F
|-
|<tt>max(x,y [,z...])</tt>
|vrátí největší ze zadaných hodnot
|&nbsp;
|*
|-
|<tt>median(x,y[,z...])</tt>
|vrátí prostřední hodnotu (medián) ze zadaných hodnot
|&nbsp;
|*
|-
|<tt>min(x,y[,z...])</tt>
|vrátí nejmenší ze zadaných hodnot
|&nbsp;
|*
|-
|<tt>mode(x,y[,z...])</tt>
|vrátí nejčetnější hodnotu (modus) ze zadaných hodnot
|&nbsp;
|*
|-
|<tt>not(x)</tt>
|1 pokud ''x'' je nula, jinak 0
|&nbsp;
|&nbsp;
|-
|<tt>rand(low,high)</tt>
|vrátí náhodnou číselnou hodnotu z intervalu [low; high]
|&nbsp;
|&nbsp;
|-
|<tt>round(x)</tt>
|zaokrouhlí ''x'' na nejbližší celočíselnou hodnotu
|&nbsp;
|I
|-
|<tt>sin(x)</tt>
|sinus z ''x'' (''x'' ve stupních)
|<tt>sin(x)</tt>
|F
|-
|<tt>sqrt(x)</tt>
|odmocnina z ''x''
|<tt>sqrt(x)</tt>
|F
|-
|<tt>tan(x)</tt>
|tangens z ''x'' (''x'' ve stupních)
|<tt>tan(x)</tt>
|F
|}
</center>
 
''Vysvětlivky:''
<div style="margin-left: 10px">
{|
|-
|*
|...
|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.
|-
|}
</div>
 
<center>
{|class=wikitable
|+Tab č.3: Speciální proměnné
|-
! proměnná !! význam
|-
|<tt>x()</tt>
|proměnná pro aktuální souřadnici ''x''
|-
|<tt>y()</tt>
|proměnná pro aktuální souřadnici ''y''
|-
|<tt>col()</tt>
|proměnná pro aktuální číslo sloupce
|-
|<tt>row()</tt>
|proměnná pro aktuální číslo řádku
|-
|<tt>nsres()</tt>
|proměnná pro aktuální rozlišení sever-jih
|-
|<tt>ewres()</tt>
|proměnná pro aktuální rozlišení východ-západ
|-
|}
</center>
 
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í rastrové mapy či map. 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.
 
{{fig|ZOD-cv4-test-tm1_tm2|Vstupní kanály tm1 a tm2}}
 
{{fig|ZOD-cv4-test-operace12|(a) tm1 + tm2; (b) tm1 - tm2; (c) tm1 * tm2; (d) tm1 / tm2; (e) 1.0 * tm1 / tm2}}
 
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 rastrových map a datovém typu konstant uvedených ve výrazu. Výsledkem dělení dvou celočíselných rastrových map je tedy opět mapa 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 <code>float()</code>.
 
Jako příklad můžeme vypočítat podíl sedmého a čtvrtého kanálu družicové scény LandSat-TM5:
 
<source lang="bash">
# podíl obrazových kanálů sedm a čtyři (výsledná mapa DCELL)
#
r.mapcalc 'p74 = float(tm7) / tm4'
d.rast map=p74
</source>
 
{{fig|ZOD-cv4-p74|Podíl obrazu (tm7/tm4) - tabulka barev rainbow (výchozí)}}
 
=== NDVI ===
 
Praktickou úlohou podílu obrazu je výpočet ''normalizovaného diferenčního vegetačního indexu'' ({{wikipedia|NDVI|lang=en}}).
 
<source lang="bash">
# výpočet NDVI
#
r.mapcalc 'ndvi = float(tm4 - tm3) / (tm4 + tm3)'
</source>
 
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.
 
{{fig|ZOD-cv4-ndvi-rainbow|NDVI - tabulka barev rainbow}}
 
Pro lepší vizuální interpretaci obrazu zvolíme vlastní tabulku barev, nejprve zjistíme rozsah dat ({{GrassPrikaz|r.info}} -r), pro názornější představu můžeme vytisknout statistiku rastrové vrstvy pomocí modulu {{GrassPrikaz|r.report}}:
 
<source lang="bash">
# rozsah hodnot
#
r.info -r map=ndvi
</source>
min=-1.000000
max=0.861314
 
Rastrové hodnoty se tedy pohybují v intervalu [-1.0; +0.87].
 
<source lang="bash">
# 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
</source>
+-----------------------------------------------------------------------------+
|                        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
 
<source lang="bash">
# nastavit tabulku barev ze souboru tb-ndvi.txt
#
r.colors map=ndvi color=rules rules=tb-ndvi.txt
</source>
 
{{fig|ZOD-cv4-ndvi-tb|NDVI - vlastní tabulka barev}}
 
Zobrazení legendy je ve {{grassPrikaz|wxGUI}} dostupné z nástrojové lišty mapového okna.
 
{{fig|grass-ndvi-legenda|Zobrazení NDVI a legendy v prostředí wxGUI}}
 
<!--
{{fig|wxgui-legenda|Zobrazení legendy ve {{grassPrikaz|wxGUI}}|size=640}}
 
Současné zobrazení legendy a rastrové vrstvy nabízí (za cenu omezené možnosti ovlivnit podobu legendy) modul {{GrassPrikaz|d.rast.leg}}.
 
''Poznámka:'' Modul {{grassPrikaz|d.rast.leg}} není dostupný pro {{grassPrikaz|wxGUI}}. Následující ukázka nebude v prostředí wxGUI fungovat.
 
<source lang="bash">
# současné zobrazení rastrových dat a odpovídající legendy
#
d.rast.leg map=ndvi
</source>
 
{{fig|ZOD-cv4-ndvi-legenda|Zobrazení NDVI a legendy}}
 
Flexibilnější podoba legendy vyžaduje více příkazů:
 
<source lang="bash">
# 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
</source>
 
{{fig|ZOD-cv4-ndvi-legenda1|Zobrazení NDVI a podruhé}}
-->
 
== 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 {{GrassPrikaz|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.
 
<source lang="bash">
# pomocná vrstva pro r.reclass
#
r.mapcalc 'temp1 = 100 * ndvi'
</source>
 
Pro kontrolu vytiskneme rozsah hodnot rastrové vrstvy 'temp1':
 
<source lang="bash">
r.info -r map=temp1
</source>
 
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.
 
<source lang="bash">
# reklasifikace rastrových dat (reklasifikační pravidla uložena v textovém souboru)
#
r.reclass input=temp1 output=r_ndvi rules=reklas-ndvi.txt
</source>
 
{{fig|grass-r-reclass|Definice reklasifikační tabulky v GUI}}
 
Dále nastavíme vhodnou tabulku barev
 
1 red
2 yellow
3 0 136 26
 
<source lang="bash">
# 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
</source>
 
{{fig|ZOD-cv4-r_ndvi|Reklasifikovaný NDVI}}
 
{{fig|ZOD-cv4-r_ndvi-leg|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:
 
<source lang="bash">
r.report map=r_ndvi units=c,p,h
</source>
+-----------------------------------------------------------------------------+
|                        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 {{GrassPrikaz|g.copy}}) v aktuálním mapsetu
 
<source lang="bash">
# kopie tm6@PERMANENT -> tm6_null@<mapset>
#
g.copy rast=tm6,tm6_null
</source>
 
nebo z {{grassPrikaz|wxGUI}}
 
{{fig|wxgui-copy-map|Vytvoření kopie rastrové či vektorové mapy z kontextového menu Layer Manageru|size=300}}


{{Upravit}}
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 {{GrassPrikaz|r.null}}:


< [[153YZOD Zpracování obrazových dat|Stránky předmětu]] {{bullet}} [[153YZOD Zpracování obrazových dat - cvičení 3|Předchozí cvičení]] {{bullet}} [[153YZOD Zpracování obrazových dat - cvičení 5|Další cvičení]]
<source lang="bash">
g.region rast=tm6_null
#
# nahrazení hodnoty "0" hodnotou NULL
#
r.null map=tm6_null setnull=0
</source>


__TOC__
Výpočet převzat z [http://www.grassbook.org Open source GIS: A&nbsp;GRASS GIS Approach].
== Osnova ==
 
<source lang="bash">
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'
</source>
 
{{fig|wxgui-r-mapcalc-temp|Výpočet mapy teploty povrchu v rastrovém kalkulátoru wxGUI}}
 
Před vizualizací výsledné rastrové vrstvy nastavíme vhodnou tabulku barev:
 
0  blue
15 green
30 yellow
45 red
 
<source lang="bash">
r.colors map=temp_celsius col=rules < tb-teplota.txt
d.rast temp_celsius  
</source>
 
{{fig|ZOD-cv4-temp_celsius|Teplota povrchu}}


Na závěr můžeme vzniknuvší rastrové vrstvy odstranit pomocí modulu {{GrassPrikaz|g.remove}}.


----
<source lang="bash">
< [[153YZOD Zpracování obrazových dat|Stránky předmětu]] {{bullet}} [[153YZOD Zpracování obrazových dat - cvičení 3|Předchozí cvičení]] {{bullet}} [[153YZOD Zpracování obrazových dat - cvičení 5|Další cvičení]]
# odstranění zvolených rastrových vrstev
#
g.remove rast=tm6_null,tm6rad,temp_kelvin
</source>


{{ZOD}}
{{ZOD}}
{{GRASS}}

Aktuální verze z 3. 9. 2014, 09:16


Tento text se vztahuje k předmětu, který již není součástí aktuálních studijních plánů.
Stránky navazujícího předmětu 155ZDDP jsou dostupné zde.
 


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í) - viz 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

Grafické rozhraní mapového kalkulátoru

Obecná forma výrazu: 

nova_mapa = vyraz (mapa1, mapa2, ...)

kde vyraz může zahrnovat názvy rastrových map, 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 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.
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í rastrové mapy či map. 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.

Vstupní kanály tm1 a tm2
(a) tm1 + tm2; (b) tm1 - tm2; (c) tm1 * tm2; (d) tm1 / tm2; (e) 1.0 * tm1 / tm2

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 rastrových map a datovém typu konstant uvedených ve výrazu. Výsledkem dělení dvou celočíselných rastrových map je tedy opět mapa 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 = float(tm7) / tm4'
 d.rast map=p74
Podíl obrazu (tm7/tm4) - tabulka barev rainbow (výchozí)

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 = float(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.

NDVI - tabulka barev rainbow

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
NDVI - vlastní tabulka barev

Zobrazení legendy je ve wxGUI dostupné z nástrojové lišty mapového okna.

Zobrazení NDVI a legendy v prostředí wxGUI


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
Definice reklasifikační tabulky v GUI

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
Reklasifikovaný NDVI
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) v aktuálním mapsetu 

 # kopie tm6@PERMANENT -> tm6_null@<mapset>
 #
 g.copy rast=tm6,tm6_null

nebo z wxGUI

Vytvoření kopie rastrové či vektorové mapy z kontextového menu Layer Manageru

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_null
 #
 # 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'
Výpočet mapy teploty povrchu v rastrovém kalkulátoru wxGUI

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
Teplota povrchu

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
Citováno z „https://geo.fsv.cvut.cz/index.php?title=153ZODH_/_4._cvičení&oldid=29400“