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:
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:
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 (Obr. 38), 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. 39)
Obr. 39 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. 39 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. 40.
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. 40 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. 41).
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. 41 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. 42 a Obr. 43).
Obr. 42 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. 43 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. 44.
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. 44 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. 45.
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. 45 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.