153YZOD Zpracování obrazových dat 2006 - 4. cvičení

Z GeoWikiCZ
Verze z 21. 11. 2005, 15:19, kterou vytvořil Landa (diskuse | příspěvky) (→‎seznam použitých příkazů: - naprava chybnych odkazu)

Mapová algebra, sčítání, odčítání rastrových dat, podíl obrazu, NDVI

osnova

Úvod do rastrové algebry, výpočty orientované na rastrové buňky - podíl obrazu, normalizovaný diferenční vegetační index (NDVI).

seznam použitých příkazů

rastrová algebra

Rastrová (mapová) algebra je poměrně komplexní téma. Vzhledem k rozsahu problematiky si v rámci tohoto kurzu nastíníme pouze její základy. Naše znalosti v této oblasti zpracování obrazu budeme postupně během následujících cvičení dále prohlubovat. Rastrovou algebru má v GRASSu ve své moci komplexní modul r.mapcalc. Jde o "mocný" modul s rozsáhlou manuálovou stránkou (viz g.manual r.mapcalc). K tématu rastrové algebry je dostupný i speciální tutoriál. Rastrové algebře je věnována kapitola 6.16 v GRASS příručce a kapitola 15 v GIS GRASS 6.0 - Praktická rukověť začínajících uživatelů. Pro lepší pochopení probírané látky, doporučuji obě dvě kapitoly alespoň zběžně prostudovat.

Chování modulu r.mapcalc je celkem standardní. Lze ho ovládat ve dvou různým módech - tzv. interaktivním, kdy spustíte modul bez parametrů a zadáváte jednotlivé výrazy, pro ukončení slouží klíčové slovo "end". Tedy podobně jako při zadávání pravidel pro vlastní tabulku barev v modulu r.colors.

#interaktivní mód modulu r.mapcalc
#
GRASS:~ > 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 (uzavřít do uvozovek!), nebo je uložit do textového souboru (například při opětovném použití) a přesměrovat do modulu přes standardní vstup (podobně jako tabulku barev do r.colors). Např. (podily.txt):

#zadání rovnice jako parametr
#
GRASS:~ > r.mapcalc 'p12 = tm1 / tm2'
GRASS:~ > r.mapcalc 'p23 = tm2 / tm3'
#rovnice uložené v textovém souboru
#
GRASS:~ > r.mapcalc < podily.txt

My bude upřednostňovat způsob zadávání výrazů jako parametru.

Poznámka: Mějte na paměti, že i tento modul provádí výpočet pouze v aktuální regionu!

Obecně výraz pro r.mapcalc vypadá následovně:

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

kde vyraz může zahrnovat názvy rastrových mapových vrstev, celočíselné či reálné konstanty a různé 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 této části budeme hovořit pouze o výpočtech pixelově orientovaných (hodnota výsledné buňky vzniká na základě hodnot jednotlivých pixelů vstupní mapy/map). V některém z příštích kurzů narazíme na výpočty, které se provádí v rámci množiny pixelů pokryté tzv. kernelem - plovoucím oknem.

Vybereme-li dostatečný detail, můžeme velmi snadno výsledky zvolených operací "překontrolovat".

Soubor:ZOD-cv4-test-tm1.png
Obr č.1a Vstupní kanál tm1
Soubor:ZOD-cv4-test-tm2.png
Obr č.1b: Vstupní kanál tm2
Soubor:ZOD-cv4/test-operace12.png
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 k chování modulu:

  • Typ výsledné mapy (celočíselná; hodnoty s plovoucí desetinnou čárkou) závisí na typu vstupních map a datovém typu konstant uvedených ve výrazu. Výsledkem dělení dvou celočíselných map je tedy opět mapa celočíselná!
  • Vstupují-li do výpočtu rastrové buňky s hodnotou NULL (žádná data) je výsledkem opět hodnota NULL.

Asi nemá význam se dále potácet na poli pusté teorie, přejděme tedy raději k praktickým úlohám. Modul r.mapcalc budeme v dalších kurzech hojně používat, naši pozornost si rozhodně zaslouží...

podíl obrazů

Podíl obrazu je základní a velmi jednoduchou algebraickou metodou používanou pro detekci příznaků, zvýraznění obrazu, výpočet vegetačních indexů a řadu dalších úloh. Prakticky se podíl obrazu provádí pomocí komplexního modulu pro mapovou algebru - r.mapcalc. Velmi důležité je doplnit rovnici o násobení konstantou "1.0". Dělíme totiž celočíselná rastrová data. Výsledek je tedy, jak již víme, opět celočíselný! Trik s vynásobením hodnotou "1.0" má za následek vytvoření výstupní rastrové vrstvy obsahující hodnoty s plovoucí desetinnou čárkou (DCELL).

Jako příklad můžeme vypočítat podíl sedmého a čtvrtého kanálu družicové scény LandSat-TM5:

#jistota...
#
GRASS:~ > g.region rast=tm4
#podíl obrazových kanálů sedm a čtyři (výsledná mapa DCELL)
#
GRASS:~ > r.mapcalc 'p74 = 1.0 * tm7 / tm4'
GRASS:~ > d.rast p74
Obr č.3: Podíl obrazu (tm7/tm4) - tabulka barev rainbow

Praktickou úlohou na podíl obrazu je výpočet normalizovaného diferenčního vegetačního indexu (NDVI). Příkaz pro tvorbu mapy s NDVI vypadá následovně:

#výpočet NDVI
#
GRASS:~ > 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.

Obr č.4: NDVI - tabulka barev rainbow

Pro lepší vizuální interpretaci obrazu zvolíme vlastní tabulku barev (tb-ndvi.txt), nejprve zjistíme rozsah dat (r.info), pro názornější představu můžeme vytisknout statistiku dat (r.report):

#rozsah hodnot
#
GRASS:~ > r.info ndvi


+----------------------------------------------------------------------------+
| Layer:    ndvi                           Date: Sat Oct 29 14:42:15 2005    |
| Mapset:   student                        Login of Creator: martin          |
| Location: sevcech                                                          |
| DataBase: /home/martin/grassdata                                           |
| Title:     ( ndvi )                                                        |
|----------------------------------------------------------------------------|
|                                                                            |
|   Type of Map:  cell                Number of Categories: 255              |
|   Data Type:    DCELL                                                      |
|   Rows:         1661                                                       |
|   Columns:      2223                                                       |
|   Total Cells:  3692403                                                    |
|        Projection: krovak (zone 0)                                         |
|            N:    -957500    S:   -1007318   Res: 29.99277544               |
|            E: -763855.0602047    W:    -830529   Res: 29.99277544          |
|   Range of data:    min =  -1.000000 max = 0.861314                        |
|                                                                            |
|   Data Source:                                                             |
|                                                                            |
|                                                                            |
|                                                                            |
|   Data Description:                                                        |
|    generated by r.mapcalc                                                  |
|                                                                            |
|   Comments:                                                                |
|    1.000000 * (tm4 - tm3) / (tm4 + tm3)                                    |
|                                                                            |
+----------------------------------------------------------------------------+

Rastrové hodnoty se tedy pohybují v intervalu od [-1.0; +0.87]. Všimněme si navíc typu dat - DCELL (rastrová data s plovoucí desetinnou čárkou). V komentáři (Comments) je navíc uveden výraz, na základě něhož tato data vznikla.

#vytisknout statistiku (počet buněk, procenta, plocha v hektarech); přerozdělit na 25 intervalů
#
GRASS:~ > 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
#
GRASS:~ > r.colors ndvi col=rules < tb-ndvi.txt
#... a překreslit obsah grafického okna
#
GRASS:~ > d.redraw
Obr č.5: NDVI - vlastní tabulka barev

Na tomto místě si ukážeme, jak zobrazit legendu k rastrové vrstvě. Pro zobrazení samostatné legendy slouží modul d.legend, současné zobrazení legendy a rastrové mapové 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
#
GRASS:~ > d.rast.leg map=ndvi
Obr č.6: Zobrazení NDVI a legendy

Poznámka: Jelikož d.rast.leg rozdělí GRASS monitor na dvě různé části (rámce), je nutné jej posléze vyčistit příkazem d.erase -f.

Požadujeme-li více flexibilní výstup, je to pěkná hromádka příkazů, pouze pro ilustraci:

#rozdělit grafické okno na dva rámce
GRASS:~ > d.frame f1 at=0,100,0,15;d.frame f2 at=0,100,15,100</tr>
#v jednom zobrazit ndvi a popisek
#
GRASS:~ > echo "NDVI" | d.text size=5 col=black;d.rast ndvi
#a v druhém legendu (pouze pro rozsah [0;1])
#
GRASS:~ > d.frame -s f1;d.legend ndvi range=0,1 labelnum=7
Obr č.7: Zobrazení NDVI a podruhé

reklasifikace rastrových dat

Chceme-li na základě NDVI provést primitivní klasifikaci na tři základní třídy (bez vegetace, vodní plochy / plochy s minimální vegetací / plochy pokryté vegetací), postupujeme takto:

Protože 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.

#trik pro r.reclass
#
GRASS:~ > r.mapcalc 'temp1 = 100 * ndvi'

Pro kontrolu vytiskneme rozsah hodnot rastrové vrstvy temp1:

GRASS:~ > r.info temp1 | grep Range
|   Range of data:    min =  -100.000000 max = 86.131387

Připravíme si tzv. reklasifikační tabulku (reklas-ndvi.txt), pro bližší informace odkazuji na manuálové stránky modulu r.reclass:

-100 thru 5   = 1 bez vegetace, vodni plochy
 5   thru 35  = 2 plochy s minimalni vegetaci
 35  thru 87  = 3 plochy pokryte vegetaci

Na tomto místě se sluší krátký popis: na každém řádku je uvedeno jedno reklasifikační pravidlo, první řádek například "říká": proveď přiřazení hodnot vstupních dat v intervalu od -100 do 5 kategorii (hodnotě) 1, poté následuje textový popisek dané kategorie. Ve svém důsledku bude reklasifikovaná, výstupní rastrová mapa obsahovat celkem tři kategorie (tj. hodnoty). Reklasifikaci provedeme příkazem, předpokládejme, že máme reklasifikační tabulku předpřipravenu v textovém souboru:

#reklasifikace rastrových dat
#
GRASS:~ > r.reclass input=temp1 output=r_ndvi < reklas-ndvi.txt</tr>

Před samotným zobrazením reklasifikovaných dat nastavíme vhodnou tabulku barev (tb-r_ndvi.txt):

1 red
2 yellow
3 0 136 26
#nastavení vhodné tabulky barev
#
GRASS:~ > r.colors r_ndvi col=rules < tb-r_ndvi.txt
#zobrazení reklasifikovaných dat
GRASS:~ > d.rast r_ndvi
GRASS:~ > d.rast.leg r_ndvi
Obr č.8: Reklasifikovaný NDVI
Obr č.9: Reklasifikovaný NDVI s legendou

Důležité je poznamenat, že reklasifikační modul nevytváří ve skutečnosti plnohodnotnou rastrovou vrstvu (šlo by o zbytečné plýtvání diskového prostoru), pouze zapíše zvolená reklasifikační pravidla. Tento fakt je dobré mít na paměti, například to vede k tomu, že nelze smazat podkladovou vrstvu, která sloužila jako vstup do reklasifikace před samotnou reklasifikovanou mapou.

A ještě krátkou statistiku reklasifikované mapy:

GRASS:~ > 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í mapy s teplotou povrchu

Před samotným výpočtem vytvoříme kopii šestého kanálu LandSat-TM5 (seznámíme se tedy s novým modulem - g.copy určeným pro kopírování jednotlivých mapových vrstev) - tj. kopii tm6 z mapsetu PERMANENT, kterou uložíme v aktuálním mapsetu (student) pod jménem tm6_null:

#kopie tm6@PERMANENT -> tm6_null@student
#
GRASS:~ > g.copy rast=tm6,tm6_null

Nyní nahradíme hodnotu "0" hodnotou NULL (žádná data) - nulové hodnoty obsahuje levý dolní roh snímku, který nenese žá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:

#pro jistotu nastavit aktivní region
#
GRASS:~ > g.region rast=tm6
#nahrazení hodnoty "0" hodnotou NULL
#
GRASS:~ > r.null map=tm6_null setnull=0

Dle zadaných vzorců (převzato z Open source GIS: A GRASS GIS Approach) provedeme sérii výpočtů:

GRASS:~ > r.mapcalc 'tm6rad = 1.238 + (15.6 - 1.238) * tm6_null / 255'</tr>
GRASS:~ > r.mapcalc 'temp_kelvin = 1260.56 / (log (607.76 / tm6rad + 1))'</tr>
GRASS:~ > r.mapcalc 'temp_celsius = temp_kelvin - 273.15'</tr>

Před vizualizací výsledné mapy nastavíme vhodnou tabulku barev (tb-teplota.txt):

GRASS:~ > r.colors map=temp_celsius col=rules < tb-teplota.txt
GRASS:~ > d.rast temp_celsius	  
Obr č.10: Teplota povrchu

Zde se hodí představit další všeobecný modul a to g.remove. Jak jeho název napovídá, slouží pro mazaní jednotlivých mapových vrstev. Nás v podstatě budou zajímat pouze parametry rast pro rastrová a vect pro vektorová data. Nic nám tedy nebrání odstranit již nepotřebné rastrové vrstvy, které vznikly během výpočtu temp_celsius:

#odstranění zvolených rastrových mapových vrstev
#
GRASS:~ > g.remove rast=tm6_null,tm6rad,temp_kelvin