NCSU GIS/MEA582:
Geospatial Modeling and Analysis

This is an unmaintained course material, please see current material at:

Erosion modeling

Resources:

Text files with recode rules and color rules:

Start GRASS GIS

In startup pannel set GIS Data Directory to path to datasets, for example on MS Windows, C:\Users\myname\grassdata. For Project location select nc_spm_08_grass7 (North Carolina, State Plane, meters) and for Accessible mapset create a new mapset (called e.g. HW_erosion) and click Start GRASS. 
grass70

Change working directory:
Settings > GRASS working environment > Change working directory > select/create any directory
or type cd (stands for change directory) into the GUI command console and hit Enter: 

cd
Download all text files with recode rules, color rules and site locations (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.

Compute soil detachment using USLE3D

Compute topographic potential (LS factor)

Compare impact of the power function exponents on the erosion pattern.
Which equation represents conditions when water flow has greater impact and which models stronger influence of slope? 
g.region raster=elev_lid792_1m -p
r.slope.aspect elev_lid792_1m slope=slope_1m aspect=aspect_1m
r.flow elev_lid792_1m flowaccumulation=flowacc_1m
r.mapcalc "lsfac3d_1m = 1.2 * pow(flowacc_1m * 1./22.1,0.2) * pow(sin(slope_1m)/0.09,1.3)"
r.mapcalc "lsfac3d_s1_1m = 1.5 * pow(flowacc_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.rast lsfac3d_1m
d.vect elev_lid792_cont1m
d.legend lsfac3d_1m at=2,50,2,6
d.out.file lsfac_s13

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

Compute soil detachment for spatially variable land cover

Set region to rural area and recode landcover_1m to cfactor using the r.recode module. Assign special color table and display result.
g.region rural_1m -p
r.recode landcover_1m output=cfactor_1m rules=cfac_rules.txt
r.colors map=cfactor_1m rules=cfac_color.txt
d.rast cfactor_1m

Compute the USLE3D equation using map algebra, cfactorbare_1m is the same as cfactor_1m, cfactorgrow_1m has landuse recoded for summer conditions.
Specify units of the raster maps using the r.support module.
Compare erosion rates and distribution for winter (bare) and summer 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 and display the new ones: 
d.rast soillossbare_1m
d.legend soillossbare_1m at=2,50,2,6
d.out.file soillossbare
d.rast soillossgrow_1m
d.legend soillossgrow_1m at=2,50,2,6
d.out.file soillosgrow
r.univar soillossbare_1m
r.univar soillossgrow_1m

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

Compute net erosion/deposition maps (using USPED)

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

First compute sediment flow and its components in x, y directions: 

r.mapcalc "sedflow_1m = 270. * soils_Kfactor * cfactorgrow_1m * flowacc_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 and then change in the direction of flow using divergence.

r.slope.aspect qsx dx=qsx_dx
r.slope.aspect 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=-100,100

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

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.

Compute univariate statistics: 

r.univar erdep

Look for line with sum: 

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 do the division.