Metoda SCS CN

Teoretický popis metody, základní symboly použité v dalším textu, popis vstupních dat a navrhovaný postup je součástí školení GRASS GIS pro pokročilé.

Postup zpracovaní v QGIS

Krok 1

Sjednocení hlavních půdních jednotek a komplexního průzkumu půd

V prvním kroku založíme projekt a pomocí add_vector Přidaní vektorové vrstvy a add_csv Přidat vrstvu s odděleným textem vložíme do panelu vrstev vstupní data - vektorová data ve formátu Esri Shapefile hpj.shp, kpp.shp, landuse.shp, povodi.shp a soubory s odděleným textem hpj_hydrsk.csv, kpp_hydrsk.csv, lu_hydrsk_cn.csv, viz školení QGIS pro začátečníky.

Následně sjednotíme vrstvu hlavních půdních jednotek (hpj) a komplexního průzkumu půd (kpp). Využijeme nástroj geoprocessingu union Sjednotit, který najdeme v záložce Vektor ‣ Nástroje geoprocessingu. Vznikne nová vektorová vrstva hpj_kpp.

Krok 2

Připojení informací o hydrologické skupině

Tabulku hpj_hydrsk můžeme připojit k atributům nové vektorové vrstvy hpj_kpp pomocí společného sloupce HPJ (Obr. 113). Pravým tlačítkem myši v panelu vrstev u hpj_kpp zvolíme Vlastnosti a v dialogovém okně přejdeme do záložky join Připojení. Kliknutím na ikonku plus spustíme dialogové okno s nastavením pro připojení (Obr. 114), viz kapitola Připojení tabulkových dat ve školení QGIS pro začátečníky.

../_images/at_pred_join.png

Obr. 113 Společný atribut HPJ a hydrologické skupiny hlavních půdních jednotek.

../_images/at_join.png

Obr. 114 Připojení tabulky k vektorové vrstvě na zákládě společného atributu.

Tímto způsobem připojíme tabulky s informacemi o hydrologických skupinách (Obr. 115).

../_images/tab_pripojene.png

Obr. 115 Zobrazení připojených vrstev ve vlastnostech vrstvy hpj_kpp.

Poznámka

V některých případech nemusí připojení tabulek proběhnout korektně, např. pokud se liší datový typ sloupečků použitých pro připojení. Typicky textové vs. číselné pole, např. ‚09‘ vs. ‚9‘. V tomto případě je nutné data před přípojením sjednotit.

Poté otevřeme atributovou tabulku hpj_kpp, zapneme editační mód ikonkou edit Přepnout editaci a pomocí kalkulačky polí kalk Otevřít kalkulačku polí vytvoříme nový atribut. Použijeme připojené atributy o hydrologické skupině (hpj_HydrSk z hlavních půdních jednotek a kpp_HydrSk z komplexního průzkumu půd). Primárně použijeme hydrologickou skupinu pro hlavní půdní jednotky. Kde informace není dosupná - hodnota NULL, tam použijeme kpp_HydrSk (Obr. 116) a výsledek zobrazíme (Obr. 117).

CASE WHEN "hpj_HydrSk" IS NULL THEN "kpp_HydrSk" ELSE "hpj_HydrSk" END
../_images/at_hydrsk_kalk.png

Obr. 116 Vytvoření atributu s informacemi o hydrologické skupině pro elementární plochy.

../_images/hydrsk.png

Obr. 117 Hydrologické skupiny elementárních ploch v zájmovém území.

Při pohledu na legendu na Obr. 117 je možno si všimnout, že kódy hydrologických skupin jako (A)B, A(B), AB a podobně by bylo vhodné sjednotit. K tomu použijeme editační mód a atributové dotazy. V hlavní liště anebo v liště atributové tabulky zvolíme volbu select-attr Vybrat prvky pomocí vzorce pomocí které vybereme elementární plochy s hydrologickou skupinou (A)B a A(B), potom zapneme editační režim, spustíme kalk Kalkulačka polí a aktualizujeme existujíce atributy hydrsk vybraných prvků (Obr. 118). Obdobně postupujeme pro další kódy. Výsledek je prezentován na Obr. 119.

../_images/kalk_AB.png

Obr. 118 Sjednocení hodnot atributů pomocí kalkulátoru polí.

Poznámka

Pro sjednocení hodnot je možno použít také výraz CASE:

CASE WHEN "hydrsk"  =  'B(C)' THEN replace("hydrsk",'B(C)','BC') ELSE "hydrsk" END

a

CASE WHEN "hydrsk"  =  'C(D)' THEN replace("hydrsk",'C(D)','CD') ELSE "hydrsk" END
../_images/hydrsk_ok.png

Obr. 119 Sjednocené hydrologické skupiny elementárních ploch v zájmovém území.

Do této fáze je možné používat QGIS relativně bez problémů. Dále však budeme přidávat informace o využití území pro každou elementární plochu pomocí operace průniku. Při větších objemech dat mohou být nástroje geoprocessingu časově náročné a nestabilní. Pro další řešení tedy použijeme výpočetně stabilnější nástroje systému GRASS GIS.

Tip

Více o systému GRASS v rámci školení GRASS GIS pro začátečníky.

Vytvoření lokace a mapsetu

Data, ke kterým GRASS přistupuje, udržuje v pevné 3-úrovňové struktuře (databáze, lokace a mapset), viz Struktura dat - koncept lokací a mapsetů ze školení GRASS GIS pro začátečníky. Z hlavní lišty menu vybereme Zásuvné moduly ‣ GRASS ‣ Nový mapset.

../_images/menu_mapset.png

Obr. 120 Zásuvný modul GRASS - vytvoření nového mapsetu.

V dialogovém okně se objeví předvolená cesta k hlavnímu adresáři databáze GRASS (obvykle adresář s názvem grassdata). V případě, že tento adresář obsahuje již nějaké lokace, vybereme tu, ve které chceme pracovat anebo si vytvoříme novou. Nastavíme souřadnicový systém a výpočetní region (viz školení GRASS GIS pro začátečníky). Kromě mapsetu PERMANENT, který se vytvoří automaticky, je vhodné zadat i název nového mapsetu, ve kterém budou probíhat výpočty. Mapset se automaticky otevře jako aktuální mapset. V záložce Region dialogového okna nástrojů GRASS je možné měnit rozsah výpočetní oblasti výběrem v mapovém okně QGIS pomocí Select the extent by dragging on canvas (Obr. 121). Zároveň je zde možno nastavit prostorové rozlišení.

../_images/n_mapset.png

Obr. 121 Vytvoření lokace a mapsetu, nastavení výpočetní oblasti a prostorového rozlišení.

Krok 3

Průnik vrstvy hydrologických skupin s vrstvou využití území

Zájmové území potřebujeme rozdělit na více elementárních ploch. Vrstvy hpj_kpp a landuse, pro které vytvoříme průnik, musíme nejprve naimportovat do mapsetu. Import dat zajišťuje více nástrojů, tzv. modulů (import). Použijeme například modul v.in.ogr.qgis, který umožňuje načítat vrstvy (jakoby) z prostředí QGIS. Názvy vrstev zachováme stejné.

../_images/v_in_ogr_qgis.png

Obr. 122 Možnosti importu vektorových vrstev do GRASS mapsetu v prostředí QGIS.

Tip

Pokud chceme oveřit, zda se zadané vrstvy po importu v mapsetu nacházejí, použijeme shell. Kliknutím na grass_shell GRASS shell spustíme příkazový řádek. Obsah konkrétního mapsetu vypiše modul g.list. Pro výpis vektorových vrstev v aktuálním mapsetu zadáme g.list vector. Parametry nástroje lze zadat i interaktivně v grafickém dialogu nástroje, který se vyvolá pomocí přepínače --ui, např. g.list --ui.

Dokumentaci a povinné parametry každého modulu lze zobrazit zadáním příkazu man před název modulu, například man g.list.

Operaci překrytí, resp. určení průniku vektorových vrstev, zajišťuje modul v.overlay.and v.overlay.and, který spustíme z Vektor ‣ Prostorová analýza ‣ Překrytí (Obr. 123). Výslednou vrstvu průniku nazveme hpj_kpp_landuse.

Poznámka

Počet záznamů v atributové tabulce se průnikem prvků výrazně zvýší. Což je zapříčiněno hlavně tím, že QGIS zobrazuje záznamy pro multiprvky jako duplicitní.

../_images/v_overlay_and.png

Obr. 123 Modul pro určení průniku dvou vektorových vrstev.

Tip

V příkazovém řádku můžeme vypsat například:

  • seznam tabulek v aktuálním mapsetu, resp. jejich názvy: db.tables
  • seznam atributů konkrétní tabulky: db.columns table=NAZEVTABULKY
  • počet záznamů v tabulce: db.select sql='select count(*) from NAZEVTABULKY'

Příklad použití GRASS shell je znázorněn na Obr. 124. Modul v.db.select v.db.select vypíšeme hodnoty atributů, modulem v.db.select v.db.select.where je možné zadat i podmínku.

../_images/gshell_db_columns.png

Obr. 124 Zobrazení tabulek a záznamů v příkazovém řádku.

Modul v.out.ogr umožňuje exportovat atributovou tabulku do různých formátů a dále s ni pracovat. Na Obr. 125 je znázorněn export do bežného formátu CSV.

../_images/db_export.png

Obr. 125 Export atributů do formátu CSV.

Krok 4

Připojení hodnot odtokové křivky \(CN\)

V dalším kroku je potřeba vytvořit atribut, který bude obsahovat údaje o využití území a o hydrologické skupině půdy dané elementární plochy ve tvaru VyužitíÚzemí_HydrologickáSkupina.

Vytvoříme nový atribut pomocí modulu v.db.addcolumn v.db.add.column, který nazveme landuse_hydrsk (Obr. 126). Poté doplníme hodnoty atributu s využitím modulu v.db.update v.db.update_op jako výsledek operace v rámci jedné atributové tabulky. Hodnotu zadáme ve tvaru b_LandUse||'_'||a_hydrsk.

../_images/v_db_addcolumn.png

Obr. 126 Přidání atributu do atributové tabulky s datovým typem text.

Poznámka

Výsledek můžeme zkontrolovat v příkazovém řádku zadaním

v.db.select map=hpj_kpp_landuse columns=cat,b_LandUse,a_hydrsk,landuse_hydrsk where=cat=1

cat|b_LandUse|a_hydrsk|landuse_hydrsk
1|OP|B|OP_B

Dále do mapsetu modulem db.in.ogr importujeme tabulku s hodnotami CN. Nazveme ji lu_hydrsk_cn.

Následně použijeme modul v.db.join v.db.join, pomocí kterého připojíme importovanou tabulku k vektorové vrstvě hpj_kpp_landuse a to kvůli přiřazení hodnot CN ke každé elementární ploše řešeného území, viz. Obr. 127. Obsah výsledné tabulky je možno oveřit v příkazovém řádku pomocí v.db.select map=hpj_kpp_landuse where=cat=1.

Důležité

Jednotlivé atributy v tabulkách, které spojujeme, nesmí obsahovat stejné názvy atributů. Tento problém lze vyřešit zavoláním modulu v.db.join z GUI systému GRASS a volbou subset_columns, která ale v dialogu GRASS pluginu chybí. Tento problém lze obejít nativním dialogem nástroje vyvolaného příkazem v.db.join --ui.

../_images/v_db_join.png

Obr. 127 Připojení tabulky k existující atributové tabulce vektorové vrstvy.

Poznámka

Tento způsob spojení atributových dat je alternativou k operaci záložky join Připojení ve vlastnostech vektorové vrstvy, viz Krok 2.

Krok 5

Sjednocení průniku vrstvy hydrologických skupin a využití území s vrstvou povodí

Hodnoty návrhových srážek s různou dobou opakovaní do vrstvy přidáme pomocí modulu v.overlay.or v.overlay.or. Sjednocení předchází import vrstvy povodí s informacemi o srážkách do mapsetu, přičemž postup je obdobný jako při importu vektorových vrstev v úvodní části.

Ukázka záznamu (vybrané atributy) atributové tabulky nově vytvořené vektorové vrstvy hpj_kpp_lu_pov pro 2-letý úhrn srážek v mm s dobou trvaní srážky 120 min:

v.db.select map=hpj_kpp_lu_pov columns=cat,a_CN,b_H_002_120 where="cat=1"

cat|a_CN|b_H_002_120
1|80|21.6804582207

Přehled o tom, jak se změnil počet plošných prvků ve vrstvě hpj_kpp_landuse po sjednocení s vrstvou povodí, dostaneme jako výstup modulu v.info, viz. Vektor ‣ Zprávy a statistiky. Standardní zobrazení informací je uvedeno Obr. 128.

../_images/v_info.png

Obr. 128 Výpis základních informací o vektorové vrstvě pomocí modulu v.info.

Tip

Z příkazového řádku je možno spustit nativní grafické uživatelké rozhraní systému GRASS příkazem g.gui. Taktéž je možné zapnout mapové okno (příkaz d.mon), vykreslit v něm konkrétní rastrovou (d.rast) anebo vektorovou (d.vect) vrstvu, přidat měřítko (d.barscale) či legendu (d.legend). Příkazem d.rast.leg vykreslíme rastrovou vrstvu i s legendou.

Dále budeme pracovat především s hodnotami CN. Pro další operace je potřeba, aby typ tohoto atributu byl číselný, na což použijeme funkci cast(). Vytvoříme tedy nový atribut CN s datovým typem integer.

Poznámka

Vektorovou vrstvu hpj_kpp_landuse je možno převést na rastrovou vrstvu s hodnotami CN a zobrazit v mapovém okně systému GRASS. Začneme vytvořením nového atributu typu integer (modul v.db.addcolumn), pokračujeme jeho editací v.db.update_op a následně spustíme modul v.to.rast.attr v.to.rast.attr, viz. Obr. 129. Pomocí příkazů d.mon wx0, d.rast.leg cn, d.barscale a d.vect povodi type=boundary zobrazíme mapu s CN včetně měřítka legendy v překrytí s vektorovou vrstvou povodí.

../_images/v_to_rast_cn.png

Obr. 129 Konverze vektorové mapy na rastrovou na základě atributu.

Krok 6

Výpočet výměry elementárních ploch, parametru \(A\) a parametru \(I_a\)

Pro každou elementární plochu vypočítame její výměru, parametr \(A\) a \(I_a\).

\[A = 25.4 \times (\frac{1000}{CN} - 10)\]
\[I_a = 0.2 \times A\]

Do atributové tabulky hpj_kpp_lu_pov přidáme nové atributy typu double, konkrétně vymera, A, I_a. Poté vypočítame jejich příslušné hodnoty. Postupujeme obdobně jako při tvorbě atributu s hodnotami o využití území a hydrologické skupině (landuse_hydrsk), přičemž pro jejich výpočet použijeme matematické operáce jako sčítaní, odčítaní, násobení a podobně (Obr. 130 a Obr. 131). Pro určení plochy každé elementární plochy využijeme modul z kategorie Vektor ‣ Zprávy a statistiky, modul v.to.db v.to.db.

../_images/add_columns.png

Obr. 130 Vytvoření více atributů najednou s využitím v.db.addcolumn.

../_images/area_A.png

Obr. 131 Výpočet výměry modulem v.to.db a parametru A modulem v.db.update_op.

Poznámka pro pokročilé

V příkazovém řádku by tyto kroky vypadaly následovně:

v.db.addcolumn map=hpj_kpp_lu_pov columns="vymera double,A double,Ia double"
v.to.db map=hpj_kpp_lu_pov option=area columns=vymera
v.db.update map=hpj_kpp_lu_pov column=A value="24.5 * (1000 / CN - 10)"
v.db.update map=hpj_kpp_lu_pov column=I_a value="0.2 * A"

Krok 7

Výpočet parametru \(H_o\) a parametru \(O_p\) pro každou elementární plochu

Znázornění vektorové vrstvy povodí s návrhovými srážkami v prostředí QGIS je uvedeno na Obr. 132 (maximální hodnota atributů H_002_120 představuje 23 mm). Histogramy je možné vykreslit v záložce diagram Diagramy ve vlastnostech konkrétní vrstvy.

../_images/navrhove_zrazky.png

Obr. 132 Zobrazení povodí IV. řádu s návrhovými srážkami.

Vypočítáme výšku přímého odtoku v mm jako parametr \(H_o\) a objem jako parametr \(O_{p}\).

\[ \begin{align}\begin{aligned}H_O = \frac{(H_S − 0.2 \times A)^2}{H_S + 0.8 \times A}\\O_P = P_P \times \frac{H_O}{1000}\end{aligned}\end{align} \]

V dalších krocích budeme uvažovat průměrný úhrn návrhové srážky \(H_{s}\) = 32 mm. Při úhrnu s dobou opakovaní 2 roky (atribut H_002_120) či dobou 5, 10, 20, 50 anebo 100 roků by byl postup obdobný.

Důležité

Hodnota v čitateli vztahu pro \(H_o\) musí být kladná, resp. nelze umocňovat záporné číslo. V případě, že čitatel je záporný, výška přímého odtoku je rovná nule. Pomůžeme si novým atributem v atributové tabulce, který nazveme HOklad.

Postupujeme obdobně jako na Obr. 130 a Obr. 131 anebo pomocí příkazového řádku.

v.db.addcolumn map=hpj_kpp_lu_pov columns="HOklad double, HO double, OP double"
v.db.update map=hpj_kpp_lu_pov column=HOklad value="(32 - (0.2 * A))"

Záporným hodnotám HOklad přiřadíme konstantu 0 modulem v.db.update v.db.update_query (Obr. 133). Atributy HO a OP vyplníme modulem v.db.update v.db.update_op.

v.db.update map=hpj_kpp_lu_pov column=HO value='(HOklad * HOklad)/(32 + (0.8 * A))'
v.db.update map=hpj_kpp_lu_pov column=OP value="vymera * (HO / 1000)"
../_images/v_db_update_query.png

Obr. 133 Přiřazení konstatní hodnoty atributu v případě splnění podmínky dotazu modulem v.db.update_query.

Tip

Přiřazení konstanty 0 pro záporné HOklad je možno zkontrolovat tak jako je prezentovano na Obr. 134.

../_images/ho_klad.png

Obr. 134 Kontrola editace záporných hodnot v příkazovém řádku.

Krok 8

Vytvoření rastrových vrstev výšky a objemu přímého odtoku

Modulem v.to.rast.attr v.to.rast.attr vytvoříme z vektorové vrstvy hpj_kpp_lu_pov rastrové vrstvy ho a op. Výsledky vizualizované v prostředí QGIS jsou uvedeny na Obr. 135.

Důležité

Před samotnou rasterizací je nutné korektně nastavit výpočetní region.

../_images/ho_op.png

Obr. 135 Zobrazení výšky a objemu přímého odtoku pro elementární plochy v prostředí QGIS.

Krok 9

Výpočet průměrných hodnot výšky a objemu přímého odtoku pro povodí

V dalším kroku vypočítáme průměrné hodnoty přímého odtoku pro každé povodí v řešeném území. Modul v.rast.stats v.rast.stats počítá základní statistické informace rastrové vrstvy na základě vektorové vrstvy a ty ukladá do nových atributů v atributové tabulce. Dialogové okno je uvedeno na Obr. 136.

../_images/v_rast_stats.png

Obr. 136 Dialogové okno modulu v.rast.stats.

Vektorovou vrstvu povodí potom převedeme do podoby rastrové vrstvy, přičem jako klíčový atribut použijeme ho_average, resp. op_average. Výstup zobrazený v prostředí QGIS je na Obr. 137.

../_images/ho_op_avg.png

Obr. 137 Zobrazení průměrné výšky a objemu přímého odtoku pro povodí v prostředí QGIS.