NCSU GIS/MEA582:
Geospatial Modeling and Analysis

Erosion modeling

Resources:

Text files with recode, color, and label rules:

Start GRASS GIS

Create a new mapset (called e.g. HW_erosion) in nc_spm_08_grass7 project and change working directory:
Settings > GRASS working environment > Change working directory > select/create any directory
or type cd (stands for change directory) into the GUI Console and hit Enter: 
cd
Download all text files with recode rules, color rules and labels (see above) to the selected directory. Now you can use the commands from the assignment requiring the text file without the need to specify the full path to the file.

Equations and units used in this assignment are explained in this Brief theoretical background for GIS-based erosion modeling with RUSLE3D and USPED

Compute soil detachment using USLE3D

First, we compute topographic potential (LS factor) for soil erosion and compare impact of the power function exponents on the erosion pattern.
Which equation represents conditions when contributing area has greater impact and which models stronger influence of slope? Which result predicts higher erosion rates? 
g.region raster=elev_lid792_1m -p
r.slope.aspect elevation=elev_lid792_1m slope=slope_1m aspect=aspect_1m
r.flow elevation=elev_lid792_1m flowaccumulation=flowaccum_1m
r.mapcalc "lsfac3d_1m = 1.2 * pow(flowaccum_1m * 1./22.1,0.2) * pow(sin(slope_1m)/0.09,1.3)"
r.mapcalc "lsfac3d_s1_1m = 1.5 * pow(flowaccum_1m * 1./22.1,0.5) * pow(sin(slope_1m)/0.09,1.0)"
r.colors lsfac3d_s1_1m color=gyr -e
r.colors lsfac3d_1m raster=lsfac3d_s1_1m
Display layers and save outputs: 
d.erase
d.rast lsfac3d_1m
d.vect elev_lid792_cont1m
d.legend -d lsfac3d_1m at=2,30,2,6
d.barscale
d.out.file lsfac_s13

d.rast lsfac3d_s1_1m
d.vect elev_lid792_cont1m
d.legend -d lsfac3d_s1_1m at=2,30,2,6
d.out.file lsfac_s10

Compute soil detachment for spatially variable land cover

Recode landcover_1m to land cover erosion factor C values given in the cfac_rules.txt table using the r.recode module. Assign a custom color table and display result.
r.recode landcover_1m output=cfactor_1m rules=cfac_rules.txt
r.colors map=cfactor_1m rules=cfac_color.txt
d.rast cfactor_1m
d.legend cfactor_1m
d.out.file cfactor_bare

Compute the USLE3D equation using map algebra, cfactorbare_1m is the same as cfactor_1m with bare agricultural fields, cfactorgrow_1m has landuse recoded for growing season conditions. We use rainfall intensity factor R=270 and soil erodibility factor K extracted from the soils database.
Assign custom color ramp using the provided rules and specify units of the raster maps using the r.support module.
Compare erosion rates and distribution for bare and growing season conditions. 

r.mapcalc "soillossbare_1m = 270. * soils_Kfactor * lsfac3d_1m * cfactorbare_1m"
r.mapcalc "soillossgrow_1m = 270. * soils_Kfactor * lsfac3d_1m * cfactorgrow_1m"
r.colors soillossbare_1m rules=soilloss_color.txt
r.colors soillossgrow_1m raster=soillossbare_1m
r.support map=soillossbare_1m units="ton/(acre.year)"
r.support map=soillossgrow_1m units="ton/(acre.year)"
Remove previous layers, display the soil loss maps for two different seasons and compute the summary statistics for erosion rates: 
d.erase
d.rast soillossbare_1m
d.legend soillossbare_1m at=2,40,2,6
d.out.file soillossbare
d.rast soillossgrow_1m
d.legend soillossgrow_1m at=2,40,2,6
d.out.file soillosgrow
r.univar soillossbare_1m
r.univar soillossgrow_1m

Compute net erosion/deposition maps using USPED model

Compute net erosion/deposition maps as divergence of sediment flow (USPED).

First compute sediment flow and its x, y components in the steepest slope (aspect) direction: 

r.mapcalc "sedflow_1m = 270. * soils_Kfactor * cfactorgrow_1m * flowaccum_1m * sin(slope_1m)"
r.colors sedflow_1m raster=soillossbare_1m
d.rast sedflow_1m
r.mapcalc "qsx = sedflow_1m * cos(aspect_1m)"
r.mapcalc "qsy = sedflow_1m * sin(aspect_1m)"

Compute change of sediment flow in the x and y directions as partial derivatives of sediment flow component surfaces. We are using qsx and dsy surfaces (rasters) instead of elevation raster as input and partial derivatives dx and dy as output of the r.slope.aspect tool. The change of sediment flow in the direction of flow (aspect) is then computed as divergence resulting in net erosion and deposition map erdep. Display the map of erosion and deposition pattern using custom color ramp.

r.slope.aspect elevation=qsx dx=qsx_dx
r.slope.aspect elevation=qsy dy=qsy_dy
r.mapcalc "erdep = qsx_dx + qsy_dy"
r.info -r erdep
r.colors erdep rules=erdep_color.txt
d.rast erdep
d.legend erdep at=2,50,2,6 range=-5,5
d.out.file erosion_deposition

Optional: Display elev_lid792_1m in 3D and drape over erdep as color (switch off all layers except for elev_lid792_1m).

Compute summary statistics

Use r.recode to classify erosion/deposition and r.category to add labels (stable, high erosion, etc) to individual categories: 
r.recode erdep output=erdep_class rules=erdep_class.txt
r.category erdep_class rules=erdep_label.txt sep=:
r.report erdep_class unit=p,h,a
r.colors erdep_class color=differences
d.rast erdep_class
d.legend -c erdep_class at=60,80,2,6
d.out.file erdep_classes

Example output: 

[...]
| #| description         |  %  | hectares |  acres  |
|-4| severe erosion . . .| 0.19|  0.101300|  0.25031|
|-3| high erosion . . . .| 1.34|  0.701600|  1.73365|
|-2| moderate erosion . .| 3.89|  2.042600|  5.04726|
|-1| low erosion . . . . |19.74| 10.366000| 25.61438|
| 0| stable . . . . . . .|61.32| 32.192000| 79.54643|
| 1| low deposition . . .| 8.40|  4.407600| 10.89118|
| 2| moderate deposition | 2.49|  1.307500|  3.23083|
| 3| high deposition . . | 1.29|  0.676900|  1.67262|
| 4| severe deposition . | 0.24|  0.126100|  0.31159|
| *|no data. . . . . . . | 1.10|  0.578400|  1.42922|
|---------------------------------------------------|
|TOTAL                   |100.00| 52.500000|129.7275|
Comment on advantages, disadvantages and risks of classifying erosion/deposition data into categories.

Optional

Estimate amount of excess sediment transported out of the studied site.

Compute univariate statistics for the erdep map and extract the line with sum: 
r.univar erdep

sum: -2374.473814

The units are tons/(acre.year), but our cells are 1 m2. Therefore we have to divide by 4046 (1 acre = 4046 m2) to get total ton per year transported from the watershed. You can use the Python shell in the Python tab do the division.

Compute new DEM with erosion carved-in

r.mapcalc "elev_erodedb_1m = elev_lid792_1m-(soillossbare_1m/100.)"
Display elev_erodedb_1m in 3D and drape over soillossbare_1m as color.
To view it in 3D switch off everything except elev_erodedb_1m.
Drape soillossbare_1m as color and don't forget to set fine resolution to 1.
Use lighting from SW, z-exag at least 2.0