153YZOD Zpracování obrazových dat 2006 - 4. cvičení
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ů
- r.mapcalc
- d.rast
- r.info
- r.report
- r.colors
- d.redraw
- d.legend
- d.rast.leg
- r.reclass
- g.copy
- r.null
- g.remove
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.
| 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 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".
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
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.
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
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
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
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
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
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