Resources: ESRI virtual campus Task: Compute soil detachment and net erosion/deposition using simple GIS-based models
Create the assignment working directory (folder) .\Erosion
Check out a Spatial Analyst Extension license:
Under Tools->Extensions make sure there is a check next to 'Spatial Analyst' Select View->Toolbars->Spatial Analyst to activate the extensionSet your workspace and scratch workspace to C:\path\.Erosion
Compute soil detachment rate using USLE3DCompute a map of soil detachment rate using USLE3D equation which assumes that there is no deposition.
This represents detachment limited sediment transport by overland flow.
E = RKLSCP
- E is the annual soil detachment rate (soil loss)
- R is the rainfall
- K is the soil
- C is the land cover
- P is the preventive measures factor
- LS = 1.4 * (U/22.1)**0.4 * (sin (slope)/0.09)**1.2
LS is the topography factor, U is contributing area per unit width (function of flow accumulation)
For units and explanation of constants see Erosion modeling notes
First compute slope, and flow accumulation for the elevation raster map e_lid792 (NOT elid_792):
In ArcMap Add the raster feature layer: ncrast.mdb/e_lid792 (NOT elid_792) Open the ArcToolbox
Compute slope_1m map using the Slope function, we will use it later:
Select 'Spatial Analyst Tools->Surface->Slope' input: e_lid792 output: slope_1m keep DEGREES for units
Compute the flow accumulation map Flowacc:
Compute Fill input=e_lid792 output=Fill_e_lid Compute FlowDirection input=Fill_e_lid ouput=FlowDir_fil Compute FlowAccumulation input=FlowDir_fil output=Flowacc
LS factorBuild an expression in the Spatial Analyst Raster Calculator for lsfac map using the Flowacc and slope_1m layers you just generated.
Note that the lsfac map represents the topographic factor LS used in the USLE3D model.
Flowacc map multiplied by resolution is the contributing area per unit width U.
In the ArcToolbox, Select 'Spatial Analyst->Map Algebra->Raster Calculator
Type in the expression for lsfac layer (make sure you have spaces around operators).
Slope should be in radians, flowaccumulation is in number of cells so we multiply by resolution:
1.4 * (Power("Flowacc" * 1. / 22.1 , 0.4) * Power(Sin("slope_1m" * math.pi/180.0) / 0.09, 1.2)) Set the output raster to lsfac Click 'OK'
Soil lossNow build the expression for soil detachment rate E (soil loss) s_loss. From the formula
E = RKLSCPWe computed LS factor R factor is considered to be constant and set to 270.0 K factor values can be added as a layer ncrast.mdb/soils_kf C factior values can be added as a layer ncrast.mdb/cfactb P factor (prevention measures) is here taken to = 1.
Multiplying all the factors in Raster Calculator we will compute the raster map of erosion rate
In the ArcToolbox, Select 'Spatial Analyst->Map Algebra->Raster Calculator expression= 270.0 * "cfactb" * "lsfac" * "soils_kf" output=s_lossAdjust color table for the s_loss map - histogram equalized color ramp may be needed Find average soil loss rate in the area from the s_loss map (in Properties:Source).
USPED modelCompute net erosion and deposition as a divergence of sediment flow vector field.
Note that custom color tables that cover entire range of values including negative are important in this section.
We have already computed slope_1m and flowaccumulation, but we will also need direction of flow.
Compute aspect_1m map (direction of flow, direction of gradient vector) using the Aspect function:
'Spatial Analyst Tools->Surface->Aspect' input=e_lid792, output=aspect_1m
Build an expression for sflowtopo (topographic factor LS for sediment flow T)
using the Raster Calculator.
For the exponents use m=n=1, resolution is 1.
In the ArcToolbox, Select 'Spatial Analyst->Map Algebra->Raster Calculator expression= "Flowacc" * 1. * Sin("slope_1m" * math.pi/180.0) output raster = sflowtopo
OR use m=1.3 and n=1.2 for study areas with extensive rills.
In this case, channels/streams will have large erosion rates due to high values of flowaccumulation.
In the ArcToolbox, Select 'Spatial Analyst->Map Algebra->Raster Calculator expression= (Power("Flowacc" * 1. , 1.3) * Power((Sin("slope_1m" * math.pi/180.0)), 1.2)) output raster = sflowtopo2
Compute divergence of sediment flowFirst compute components of sediment flow in x and y direction: Calculate qsx
In the ArcToolbox, Select 'Spatial Analyst->Map Algebra->Raster Calculator expression= "sflowtopo" * "soils_kf" * "cfactb" * 270. * Cos((- "aspect_1m" + 450.) * math.pi / 180.0) Output raster = qsxCalculate qsy
In the ArcToolbox, Select 'Spatial Analyst->Map Algebra->Raster Calculator expression= "sflowtopo" * "soils_kf" * "cfactb" * 270. * Sin((- "aspect_1m" + 450.) * math.pi / 180.0) Output raster = qsy
Compute components of change in sediment flow in x and y direction as partial derivatives of sediment flow field, derived from slope and aspect - see eqs 1,2,3 from here.
In ArcToolbox 'Spatial Analyst Tools->Surface->Slope' input=qsx, output=Slope_qsx
In the same way compute Aspect_qsx, Slope_qsy and Aspect_qsy
In ArcToolbox 'Spatial Analyst Tools->Surface->Aspect' input=qsx, output=Aspect_qsx 'Spatial Analyst Tools->Surface->Slope' input=qsy, output=Slope_qsy In ArcToolbox 'Spatial Analyst Tools->Surface->Aspect' input=qsy, output=Aspect_qsyCompute qsx_dx and qsy_dy
'Spatial Analyst->Map Algebra->Raster Calculator: expression= Cos((- "Aspect_qsx" + 450) * math.pi / 180.0) * Tan("Slope_qsx" * math.pi/180.0) output=qsx_dx expression= Sin((- "Aspect_qsy" + 450) * math.pi / 180.0) * Tan("Slope_qsy" * math.pi/180.0) output=qsy_dy
Finally compute net erosion deposition:
'Spatial Analyst->Map Algebra->Raster Calculator: expression="qsx_dx" + "qsy_dy" output=erdep
Assign an appropriate color classification scheme to the erdep raster and generate a map for the report.
Note: the differences in the erdep are subtle and the default color table will almost assuredly render
the erosion deposition map in a way that little detail is revealed.
You are free to experiment with the the tools in Arc in order to define classification scheme which will reveal the erosion/deposition pattern.
Here's an option that should yield pretty good contrast between erosion and deposition:
Open Layer Properties for the erdep layer Under the Symbology tab‚ Select Classified from options (Unique Values, Classified, Stretched, Discrete) in the left column Set Classes to 11 Choose a broad, divergent color ramp from the Color Ramp drop down field that ranges from brown grading through white to green Click on Classify Edit the eleven break values (located in column on far right in Classification dialog to: -250000.00 (data minimum) -50.00 -5.00 -1.00 -0.10 0.10 1.00 5.00 50.00 330000.00 (data maximum) Click OK Click ApplyYou should now see the erosion and deposition areas rendered with pretty good contrast.