GRASS GIS / Tutoriál / GIS 2 / 2. cvičení

Z FreeGIS portál
Přejít na: navigace, hledání
předchozí cvičenístránky tutoriáludalší cvičení

Reklasifikace rastrových dat, mapová algebra

Stránka obsahuje řešení úloh z 2.cvičení předmětu GIS 2 v prostředí GRASS GIS. Pro řešení úloh (pokud není uvedeno jinak) jsou použity datové vrstvy z ArcČR 500.

Poznámky

Reklasifikace rastrových dat

Reklasifikaci rastrových dat umožňuje modul r.reclass dostupný z menu Raster → Change category values and labels → Reclassify.

Ukázka reklasifikace DMT s nadmořskou výškou nad 700m

Reklasifikační pravidla se zadávají ve formátu

<oldvalue> = <newvalue>

anebo pro rozsah

<oldvalue1> thru <oldvalue2> = <newvalue>

Mapová algebra

Pro mapovou algebru slouží modul r.mapcalc, jeho grafická nadstavba je dostupná z menu Raster → Raster map calculator.

Ukázka odvození rastrové mapy s nadmořskou výškou nad 700m na základě DMT

V rámci tohoto cvičení budeme používat výhradně funkci if() v zápisu if(a, b, c), kde 'a' je podmínka, 'b' je výsledná hodnota za předpokladu, že je podmínka splněna, jinak je výslednou hodnotou 'c'. Např. if(dmt <= 700, 0, 1), vytvoří na základně rastrové mapy dmt novou rastrovou mapu s hodnotami 1 a 0, kde hodnoty 1 označují místa s nadmořskou výškou větší než 700m. Místa, která tuto podmínku nesplňují jsou označeny hodnotou 0.

Určení výměry

Určení výměry pro jednotlivé kategorie rastrové mapy poskytuje modul r.report, který je mimojiné dostupný také z kontextového menu ve správci vrstev (viz obrázek níže).

Spuštení modulu r.report z kontextového menu správce vrstev
Určení výměry území s nadmořskou výškou větší než 700 m

Úlohy

1.

Vytvořte DMT na základě vrstevnic a výškových kót pro Ústecký kraj s prostorovým rozlišením 100m. Jaká je výměra území v ha s nadmořskou výškou větší než 700m?
Datové vrstvy: KrajePolygony (AC); Vrstevnice, VyskoveKoty → dmt
Výsledek: 50 214 ha

Nejprve vybereme polygon Ústeckého kraje

Vytvoření polygonové mapy Ústeckého kraje jako výsledek atributového dotazu

nebo z příkazové řádky

v.extract input=KrajePolygony output=uk where="nk = 'US'"

DMT na základě vrstevnic a výškových kót vytvoříme pomocí modulu v.surf.rst, který na vstupu vyžaduje bodová data. Nejprve tedy vrstevnice a výškové kóty pro Ústecký kraj (v.overlay + v.select) převedeme na 3D bodová data (v.to.3d) a poté je spojíme modulem v.patch do jedné vektorové mapy:

v.overlay ainput=Vrstevnice binput=uk@cv2 operator=and output=vrstevnice_uk olayer=0,1,0
v.select ainput=VyskoveKoty binput=uk@cv2 output=koty_uk operator=within
v.to.points input=vrstevnice_uk output=vrstevnice_body_uk
v.to.3d input=vrstevnice_body_uk output=vrstevnice_body_uk_3d column=vyska
v.to.3d input=koty_uk output=koty_uk_3d column=vyska
v.patch input=koty_uk_3d,vrstevnice_body_uk_3d output=vrstevnice_koty_3d

Před interpolací DMT z 3D bodových dat nastavíme výpočetní region (g.region) na polygon Ústeckého kraje s prostorovým rozlišením 100m. Masku nastavíme pomocí modulu r.mask, po vytvoření DMT masku odstraníme:

g.region vect=uk res=100
r.mask vector=uk
v.surf.rst input=vrstevnice_koty_3d elevation=dmtcv2
r.mask -r
Vytvořený DMT s legendou

Reklasifikační tabulka pro r.reclass:

125 trhu 700 = 0
700 thru 1221 = 1

anebo výraz pro r.mapcalc:

if(dmt <= 700, 0, 1)

2.

Jaká je výměra území v ha se sklonem svahu větším než 15 stupňů?
Datové vrstvy: slope
Výsledek: 15 931 ha

Výpočet míry svahu zajišťuje modul r.slope.aspect s parametrem slope.

r.slope.aspect elev=dmt slope=slope
Mapa míry svahu ve stupních

Reklasifikační tabulka pro r.reclass:

0 thru 15 = 0
15 thru 38 = 1

anebo výraz pro r.mapcalc:

if(slope <= 15, 0, 1)

3.

Jaká je výměra území v ha s orientací svahu na sever a zároveň se sklonem větším než 15 stupňů?
Datové vrstvy: aspect, slope
Výsledek: 1 529 ha

Výpočet orientace svahu zajišťuje modul r.slope.aspect s parametrem aspect.

r.slope.aspect elev=dmt aspect=aspect
Mapa orienace svahu s legendou

Za orientaci svahu na sever budeme považovat rozsah (67.5; 112.5).

Reklasifikační tabulka pro r.reclass:

0 thru 67.5 = 0
67.5 thru 112.5 = 1
112.5 thru 360 = 0

anebo výraz pro r.mapcalc:

if(aspect > 66.5 && aspect < 112.5, 1, 0)

Řešením bude kombinace mapy míry svahu nad 15° a orientace svahu na sever aplikací r.mapcalc:

slope_15 * aspect_sever

4.

Jaká je výměra území v ha s orientací svahu na sever anebo se sklonem větším než 15 stupňů?
Datové vrstvy: aspect, slope
Výsledek: 72 354 ha

Výraz pro r.mapcalc:

if(slope_15 == 1 || aspect_sever == 1, 1, 0)

5.

Jaká je výměra území v ha, které je viditelné z vrcholu Milešovky [S-JTSK: 986668, 770118] a zároveň má orientaci svahu na sever?
Datové vrstvy: viewshed, aspect
Výsledek: 26 629ha

Nejdříve vytvoříme bodovou mapu s vrcholem Milešovky pomocí modulu v.in.ascii, viz obrázek níže.

Vytvoření bodové vektorové mapy s vrcholem Milešovky

Dále provedene analýzu viditelnosti r.viewshed:

r.viewshed -b input=dmt output=viewshed coordinates=-770118,-986668

a modulem r.mapcalc vytvoříme výslednou rastrovou mapu:

viewshed * aspect_sever

6.

Jaká je výměra území v ha, kde jsou splněny alespoň 2 z následujících podmínek - nadmořská výška nad 700 m, sklon větší než 15 stupňů, orientace na sever?
Datové vrstvy: dmt, slope, aspect
Výsledek: 8 958 ha

Výraz pro r.mapcalc:

if((dmt_700 + slope_15 + aspect_sever) >= 2, 1, 0)

7.

Jaká je výměra území v ha, které je do 500 m od nejbližší silnice a zároveň má sklon větší než 15 stupňů?
Datové vrstvy: Silnice, slope
Výsledek: 7 198 ha

Na základě polygonu kraje (viz úloha č. 1) ořežeme liniovou vrstvu silnic (v.overlay):

v.overlay ainput=Silnice binput=uk operator=and output=silnice_uk olayer=0,1,0

Vytvoříme obalovou zónu ve vzdálenosti 500m pro silnice na území Ústeckeho kraje (v.buffer) a tuto vektorovou vrstvu nastavíme jako masku (r.mask):

v.buffer input=silnice_uk distance=500 output=silnice_uk_500
r.mask vector=silnice_uk_500
Vizualizace mapy míry svahu nad 15° s nastavením masky obalové zóny 500m od silnic

V posledním kroku masku deaktivujeme, tak abychom neovlivňovali další výpočty:

r.mask -r

8.

Jak dlouhý úsek silnice E55 v km je vidět z vrcholu Milešovky [S-JTSK: 986668, 770118]?
Datové vrstvy: dmt, Silnice
Výsledek: 65 km

Nejprve vybereme úseky silnic, které tvoří E55 (v.extract):

v.extract input=Silnice output=silnice_e55 where="mezinarodni_oznaceni =  'E55'"

Pro další řešení využijeme výsledek analýzy viditelnosti z úlohy č. 5, který modulem r.to.vect převedeme do vektorové reprezentace a na základě vzniklé polygonové mapy silnice ořežeme tak, abychom dostaly pouze úseky viditelné z vrcholu Milešovky:

r.to.vect input=viewshed output=viewshed_poly type=area
v.overlay ainput=silnice_e55 binput=viewshed_poly6 operator=and output=silnice_e55_view olayer=0,1,0

Výslednou délku můžeme určit například modulem v.to.db:

v.to.db -p -c map=silnice_e55_view option=length units=kilometers
cat|length
408|0.19434493507374
409|0.274179715197274
...
17073|3.28684997718662
total length|65.0384302001245

9.

Jaká ve výměra území v ha, kde nadmořská výška je menší než výraz "10 krát sklon svahu ve stupních"?
Datové vrstvy: dmt, slope
Výsledek: 864 ha

Výraz pro r.mapcalc:

if(dmt < 10 * slope, 1, 0)

10.

Jaká je nadmořská výška vrcholu Milešovky [S-JTSK: 986668, 770118] odvozená z DMT (správně je 836,5 m)?
Datové vrstvy: dmt, milesovka
Výsledek: 819,7 m

Odvozenou výšku vrcholu Milešovky (bodovou mapu použijeme z úlohy č. 5) určíme modulem v.what.rast:

v.what.rast -p map=milesovka raster=dmt

11.

Jaký je rozdíl celkové délky v metrech silnic 1.třídy měřeného po povrhu a jeho průmětu do roviny?
Datové vrstvy: dmt, silnice
Výsledek: 225 m

Pro řešení použijeme mapu silnic Ústeckého kraje z úlohy č. 7, z nich pomocí modulu v.buffer vybere pouze ty spadající do 1. třídy:

v.extract input=silnice_uk where="trida = 3" output=silnice_uk_t1

Nejrpve aktualizujeme délku silnic jako průmětu do roviny (v.to.db):

v.to.db map=silnice_uk option=length columns=shape_length

anebo pomocí správce atributových dat jak ukazuje obrázek níže.

Výpočet délky linovových geoprvků ve wxGUI

Délku silnic po povrchu vypočítame pomocí modulu v.drape a její hodnotu uložíme do atributové tabulky jako sloupeček 'Shape_length_dmt' (v.db.addcolumn anebo správce atributových dat, viz obrazek níže):

v.drape input=silnice_uk output=silnice_uk_dmt elevation=dmt
v.db.addcolumn map=silnice_uk_dmt column="shape_length_dmt double"                    
v.to.db map=silnice_uk_dmt option=length columns=shape_length_dmt
Přidání nového sloupečku do atributové tabulky ve wxGUI

Výsledný rozdíl můžeme určit jednoduchým SQL dotazem (db.select):

db.select sql="select sum(shape_length_dmt) - sum(shape_length) from silnice_uk_dmt"