Průměrná dlouhodobá ztráta půdy

Teoretické základy

Průměrnou roční ztrátu půdy způsobenou odtokem z pozemku určitého sklonu a způsobu využívaní je možno predikovat pomocí matematického modelu USLE, tzv. univerzální rovnice ztráty půdy:

\[G = R \times K \times L \times S \times C \times P\]

kde:

• G - průmerná dlouhodobá ztráta půdy (\(t.ha^{-1} . rok^{-1}\))

• R - faktor erozní účinnosti deště (\(MJ.ha^{-1} .cm.h^{-1}\))

• K - faktor eroze půdy (\(t.h.MJ^{-1} .cm^{-1} .rok^{-1}\))

• L - faktor délky svahu (\(-\))

• S - faktor sklonu svahu (\(-\))

• C - faktor ochranného vlivu vegetačního krytu (\(-\))

• P - faktor účinnosti protierozních opatření (\(-\))

Vstupní data

• hpj.shp - vektorová vrstva hlavních půdních jednotek z kódů BPEJ

• kpp.shp - vektorová vrstva komplexního průzkumu půd

• landuse.shp - vektorová vrstva využití území

• povodi.shp - vektorová vrstva povodí IV. řádu s návrhovými srážkami \(H_s\) (doba opakovaní 2, 5, 10, 20, 50 a 100 let)

• hpj_k.csv - číselník s kódem K pro hlavní půdní jednotky, Obr. 32 vlevo

• kpp_k.csv - číselník s kódem K pro vrstvu komplexního průzkumu půd, Obr. 32 vpravo

• lu_c.csv - číselník s kódem C pro vrstvu využití území, Obr. 32 vpravo

• dmt.tif - digitální model terénu v rozlišení 10 x 10 m, Obr. 31 vlevo

• maska.pack - oblast území bez liniových a plošných prvků prerušujících odtok, Obr. 31 vpravo

Navrhovaný postup

1. sjednocení vrstvy hlavních půdních jednotek a komplexního průzkumu půd

2. připojení hodnot faktoru K k elementárním plochám

3. průnik vrstvy elementárních ploch a využití území

4. připojení hodnot faktoru C

5. výpočet parametru KC

6. vytvoření rastrové mapy sklonu a akumulace odtoku v každé buňce

7. výpočet parametru LS

8. výpočet G a vytvoření rastru s hodnotami představující průměrnou dlouhodobou ztrátu půdy

9. výpočet průměrných hodnot G pro povodí

10. vytvoření rastrových vrstev LS a G s maskou

11. výpočet průměrných hodnot G pro povodí s maskou

Obr. 30 Grafické schéma postupu.

Znázornění vstupních vektorových dat spolu s atributovými tabulkami je totožné se vstupními vektorovými daty pro metodu SCS CN. Digitální model reliéfu a oblast řešeného území bez liniových a plošných prvků přerušující odtok (maska) jsou zobrazena na Obr. 31. Tabulky s faktory K a C jsou uvedeny na Obr. 32.

Obr. 31 Vrstva digitálního modelu terénu a oblast řešeného území bez prvků přerušujících odtok.

Obr. 32 Číselníky s hodnotami K a C.

Postup zpracování v GRASS GIS

Z digitálního modelu terénu (DMT) vytvoříme rastrovou mapu znázorňující sklonové poměry ve stupních (slope). Ten bude potřebný později pro výpočet topografického faktoru LS. V prvním kroku nastavíme výpočetní region na základě vstupního DMT a následně použijeme modul r.slope.aspect.

Tip

Podrobné informace ohledně výpočetního regionu a topografických analýz ve školení GRASS GIS pro začátečníky.

g.region raster=dmt
r.slope.aspect elevation=dmt slope=svah

Obr. 33 Hypsografické stupně (DMT) v metrech a sklonové poměry v stupních.

Dále vytvoříme vyhlazený DMT (filled), rastrovou mapu směru odtoku do sousední buňky s největším sklonem (direction) a rastrovou mapu znázorňující akumulaci odtoku v každé buňce (accumulation).

Poznámka

Pro vytvoření vyhlazeného DMT možno alternativně použít také Addons modul r.hydrodem, pro výpočet směru odtoku modul r.fill.dir a pro akumulaci odtoku r.watershed.

Todo

Tady by chtělo hlubší analýzu, v čem se moduly liší, to je otázka na kolegy z k143.

Před výpočtem si nastavíme masku podle zájmového území pomocí modulu r.mask.

r.mask raster=dmt
r.terraflow elevation=dmt filled=dmt_fill direction=dir swatershed=sink accumulation=accu tci=tci

Obr. 34 Směr odtoku ve stupních a akumulace odtoku v \(m^2\) vytvořené modulem r.terraflow.

LS faktor

LS faktor (topografický faktor) možno vypočítat podle vztahu:

\[LS = 1.6 \times (accu \times \frac{res}{22.13})^{0.6} \times (\frac{sin(slope \times \frac{pi}{180})}{0.09})^{1.3}\]

kde:

• accu je rastrová mapa akumulace odtoku

• res je prostorové rozlišení DMT

• slope je rastrová mapa míry sklonu

Poznámka

Rovnice vychází z metody Mitášová.

Pro tyto účely využijeme nástroj r.mapcalc jako hlavní nástroj mapové algebry v systému GRASS.

Tip

Více na téma mapové algebry ve školení GRASS GIS pro začátečníky.

V zápisu pro tento nástroj bude rovnice vypadat následovně:

r.mapcalc expr="ls = 1.6 * pow(accu * (10.0 / 22.13), 0.6) * pow(sin(svah * (3.1415926/180)) / 0.09, 1.3)"

Poznámka

Nastavíme vhodnou tabulku barev:

r.colors map=ls color=colors.txt

např.

0.00 128:64:64
0.01 255:128:64
0.05 0:255:0
0.10 0:128:128
0.20 0:128:255

Obr. 35 Topografický faktor LS zahrnující vliv délky a sklonu svahu.

K a C faktor

Vytvoříme vektorovou vrstvu elementárních ploch hpj_kpp_land (viz. návod na její vytvoření).

v.overlay ainput=hpj binput=kpp operator=or output=hpj_kpp
v.overlay ainput=hpj_kpp binput=landuse operator=and output=hpj_kpp_land

Do její atributové tabulky přidáme dva nové sloupce K a C. To provedeme pomocí správce atributových dat anebo modulu v.db.addcolumn.

v.db.addcolumn map=hpj_kpp_land columns="K double"
v.db.addcolumn map=hpj_kpp_land columns="C double"

Hodnotu K faktoru pro jednotlivé elementární plochy přiřadíme pomocí číselníku hpj_k.csv. Pro plochy bez hodnoty K faktoru doplníme údaje na základě půdních typů a subtypů podle komplexního průzkumu půd (tabulka kpp_k.csv). Hodnota C faktoru zemědělsky využívaných oblastí zjistíme z průměrných hodnot pro jednotlivé plodiny z tabulky lu_c.csv. Na spojení tabulek použijeme modul v.db.join.

Převodové tabulky je potřebné nejprve naimportovat do prostředí GRASS GIS. Použijeme modul db.in.ogr:

db.in.ogr in=kpp_k.csv out=kpp_k
db.in.ogr in=hpj_k.csv out=hpj_k
db.in.ogr in=lu_c.csv out=lu_c

Potom přistoupíme k připojení tabulky hpj_k k atributům vektorové vrstvy hpj_kpp_land, klíčem bude atribut HPJ.

v.db.join map=hpj_kpp_land column=a_HPJ other_table=hpj_k other_column=HPJ

Chýbějící informace hodnoty faktoru K doplníme z tabulky kpp_k SQL dotazem prostřednictvím modulu db.execute.

db.execute sql="UPDATE hpj_kpp_land SET K = (
SELECT b.K FROM hpj_kpp_land AS a JOIN kpp_k as b ON a.a_b_KPP = b.KPP)
WHERE K IS NULL"

V dalším kroku doplníme hodnoty C faktoru z importované tabulky lu_c.

v.db.join map=hpj_kpp_land column=b_LandUse other_table=lu_c other_column=LU

Údaje v atributové tabulky si zkontrolujeme, zda jsou vyplněny správně. Použijeme SQL dotaz db.select, přitom vybere jen první tři záznamy (resp. elementární plochy).

db.select sql="select cat,K,C from hpj_kpp_land where cat <= 3"

Výsledek může vypadat například takto:

cat|K|C
1|0.13|0.19
2|0.13|0.19
3|0.13|0.19
...

Poznámka

Atribut cat je hodnota, která propojuje geometrii prvků se záznamem v atributové tabulce.

Dále do atributové tabulky přidáme nový atribut KC, do kterého uložíme součin faktorů K * C. To můžeme vykonat pomocí správce atributových dat anebo modulem v.db.addcolumn v kombinaci s v.db.update.

v.db.addcolumn map=hpj_kpp_land columns="KC double"
v.db.update map=hpj_kpp_land column=KC value="K * C"

Výsledek opět zkontrolujeme.

db.select sql="select cat,K,C,KC from hpj_kpp_land where cat <= 3"

cat|K|C|KC
1|0.13|0.19|0.0247
2|0.13|0.19|0.0247
3|0.13|0.19|0.0247
...

V dalším kroku vektorovou mapu převedeme na rastrovou reprezentaci modulem v.to.rast. Pro zachovaní informací použijeme prostorové rozlišení 1 m (g.region, viz. výpočetní region ze školení GRASS GIS pro začátečníky).

g.region raster=dmt res=1
v.to.rast input=hpj_kpp_land output=hpj_kpp_land_kc use=attr attribute_column=KC

Pomocí modulu r.resamp.stats provedeme převzorkovaní na prostorové rozlišení DMT 10m a to na základě průměru hodnot vypočítaného z hodnot okolních buněk. Tímto postupom zabráníme ztrátě informací, ke kterému by došlo při přímém převodu na rastr s rozlišením 10m. Při rasterizaci se totiž hodnota buňky rastru odvozuje na základě polygonu, který prochází středem buňky anebo na základě polygonu, který zabírá největší část plochy buňky.

g.region raster=dmt
r.resamp.stats input=hpj_kpp_land_kc output=hpj_kpp_land_kc10

Na obrázku Obr. 36 je znázorněná část zájmového území, kde možno vidět rastrovou vrstvu hpj_kpp_land_kc před (vlevo dole) a po použití modulu r.resamp.stats.

Obr. 36 Část zájmového území s faktorem KC před a po převzorkovaní.

Kvůli vizualizaci nastavíme vhodnou tabulku barev a kvůli přehlednosti mapu přejmenujeme na kc modulem g.rename. Výsledek je na Obr. 37.

r.colors map=hpj_kpp_land_kc10 color=wave
g.rename raster=hpj_kpp_land_kc10,kc

Obr. 37 Faktor KC zahrnující vliv eroze půdy a vliv ochranného vlivu vegetačního pokrytu.

R a P faktor

Hodnoty těchto parametrů nebudeme odvozovat jako ty předcházející. V tomto případě jednoduše použijeme průmernou hodnotu R a P faktoru pro Českou republiku, t.j R = 40 a P = 1.

Výpočet průmerné dlouhodobé ztráty půdy

Ztrátu půdy G vypočítame modulem r.mapcalc (rmapcalc), přičemž vycházíme ze vztahu, který byl uvedený v teoretické časti školení.

r.mapcalc expr="g = 40 ∗ ls ∗ kc ∗ 1"
r.colors -n -e map=g color=corine

Pro výslednou vrstvu zvolíme vhodnou barevnou škálu, přidáme legendu, měřítko a mapu zobrazíme (Obr. 38)

Obr. 38 Vrstva s hodnotami představujícími průměrnou dlouhodobou ztrátu půdy G v jednotkách \(t.ha^{-1} . rok^{-1}\).

Poznámka

Na Obr. 38 je maximální hodnota v legendě 1. Je to pouze z důvodu, aby byl výsledek přehledný a korespondoval s barvami v mapě. V skutečnosti parametr G nabývá hodnot až 230, při takovémto rozsahu by byla stupnice v legendě jednobarevná (v našem případě červená). Změnit rozsah intervalu v legendě bylo možné nastavením parametru range, konkrétněji příkazem d.legend raster=g range=0,1.

Průměrná hodnota ztráty pro povodí

Na určení průměrné hodnoty a sumy ztráty pro každé částečné povodí využijeme modul v.rast.stats. Klíčovou vrstvou je vektorová mapa povodí povodi_iv, kde nastavíme pro nově vytvořený sloupec prefix g_. Z těchto hodnot potom modulem v.db.univar vypočítáme statistiky průměrných hodnot ztráty půdy.

v.rast.stats map=povodi_iv raster=g column_prefix=g method=average
v.db.univar map=povodi_iv column=g_average

Poznámka

Vektorová vrstva povodí musí být umístěna v aktuálním mapsetu. Pokud například pracujeme v jiném mapsetu, stačí když ji přidáme z mapsetu PERMANENT a následně v menu pravým kliknutím na mapě zvolíme Make a copy in the current mapset.

Pro účely vizualizace vektorovou vrstvu převedeme na rastr, pomocí modulu r.colors nastavíme vhodnou tabulku barev a výsledek prezentujeme, viz. Obr. 39.

v.to.rast input=povodi_iv output=pov_avg_G use=attr attribute_column=g_average
r.colors -e map=pov_avg_G color=bgyr

Obr. 39 Povodí s průměrnými hodnotami ztráty půdy

Poznámka

Z důvodu přehlednosti je opět interval v legendě upravený. Maximální hodnota průmerné ztráty půdy na povodí je až 0.74 (v jednotkách \(t.ha^{-1} . rok^{-1}\))

Zahrnutí prvků prerušujících odtok

Pro výpočet uvedený výše vychází ztráta půdy v některých místech enormně vysoká. To je způsobené tím, že ve výpočtech nejsou zahrnuté liniové a plošné prvky přerušující povrchový odtok. Těmito prvky jsou především budovy, příkopy dálnic a silnic, železniční tratě anebo ploty lemující pozemky.

Abysme zjistili přesnější hodnoty, je nutné tyto prvky do výpočtu zahrnout. Pro tento účel použijeme masku liniových a plošných prvků přerušujících odtok maska.pack a vypočítame nové hodnoty LS faktoru a ztráty půdy. Vstupem bude dmt bez prvků přerušujících odtok (Obr. 40).

Todo

Tuto část je potřeba rozšířit. Maska by se dala určit z RÚIAN, a pod.

r.unpack -o input=MASK.pack output=maska
r.mask raster=maska
r.terraflow elevation=dmt filled=dmt_fill_m direction=dir_m swatershed=sink_maccumulation=accu_m tci=tci_m

Obr. 40 Vrstva digitálního modelu terénu vstupujícího do výpočtu bez prvků přerušujících odtok.

Dále dopočítame faktor LS a následně G.

r.mapcalc expr="ls_m = pow(accu_m * (10.0 / 22.13), 0.6) * pow(sin(svah * (3.1415926/180)) / 0.09, 1.3)"
r.mapcalc expr="g_m = 40 ∗ ls_m ∗ kc ∗ 1"

r.colors map=ls_m color=wave
r.colors -n -e map=g_m color=corine

V posledním kroku vymažeme masku, výsledky zobrazíme a porovnáme (Obr. 41 a Obr. 42).

Obr. 41 Porovnání hodnot faktoru LS bez ohledu na prvky přerušující odtok (vlevo) a s prvky přerušujícími odtok (vpravo).

Obr. 42 Porovnaní výsledků USLE bez ohledu na prvky přerušující odtok (vlevo) a s prvky přerušujícími odtok (vpravo).

Průměrná hodnota ztráty pro povodí s prvky přerušujícími odtok

Opět využijeme modul v.rast.stats. Vektorové mapě povodí povodi_iv nastavíme prefix g_m pro nově vytvořený sloupec a potom modulem v.db.univar zobrazíme statistiky průměrných hodnot ztráty půdy. Výsledek v rastrové podobě je na Obr. 43.

v.rast.stats map=povodi_iv raster=g_m column_prefix=g_m method=average
v.db.univar map=povodi_iv column=g_m_average

v.to.rast input=povodi_iv output=pov_avg_G_m use=attr attribute_column=g_m_average
r.colors -e map=pov_avg_G_m color=bgyr

Obr. 43 Povodí s průměrnými hodnotami ztráty půdy s uvážením prvků, které přerušují odtok.

Na závěr vypočítáme rozdíly (modul r.mapcalc) výsledných vrstev bez a s uvážením prvků, které přerušují odtok pro faktor LS, hodnoty představující průměrnou dlouhodobou ztrátu půdy G a povodí s průměrnými hodnotami ztráty půdy G_pov. Nazveme je delta_ls, delta_g a delta_pov_avg a nastavíme barevnou stupnici differences. Viz. Obr. 44.

r.mapcalc expression=delta_ls = ls - ls_m
r.mapcalc expression=delta_g = g - g_m
r.mapcalc expression=delta_pov_avg = pov_avg_G - pov_avg_G_m

r.colors map=delta_ls color=differences
r.colors map=delta_g color=differences
r.colors map=delta_pov_avg color=differences

Obr. 44 Znázornění rozdílů rastrových vrstev LS, G a G_pov, které vznikly bez uvážení a s uvážením prvků, které přerušují odtok.

Poznámky

GRASS nabízí na výpočet USLE dva užitečné moduly r.uslek a r.usler.