(Version 1.0)
Index sl.aml: Core program, calls others as needed. fil.aml: Fills depressions. dn_slope.aml: Calculates downhill slope angle high_pts.aml: Identifies high points from the DEM. s_length.aml: Calculates cumulative downhill slope length. ls.aml: Calculates USLE LS factor values.
To run the program, just run sl.aml with the following args: /* Usage: sl<elevation grid> <slope length grid> <LS value grid> /* {FEET | METER} {cutoff value} {slope angle grid}
where
sl.aml
/* sl.aml *****************************************************************
/* The input : a grid of elevations.
/* The elevations must be in the same units as the horizontal distance.
/* The unit of measurement for the elevation grid.
/* The change in slope(as a %) that will cause the slope length
/* calculation to stop and start over.
/* The output: a grid of cumulative slope lengths,
/*: a grid of LS values for the soil loss equation,
/*: an optional grid of down hill slope angle.
/* Usage: sl <elevation grid> <slope length grid> <LS value grid>
/* {FEET | METER} {cutoff value} {slope angle grid}
&args sl_elev sl_out LS_out grd_units cutoff_value sl_angle
/* Convert user input to capital letters **********************************
&sv .grid_units = [translate %grd_units%]
/* Set default cutoff value if necessary **********************************
&if [null %cutoff_value%] or [index %cutoff_value% #] eq 1 &then
&sv .slope_cutoff_value = .5
&else
&sv .slope_cutoff_value = [calc [value cutoff_value] / 100]
/* Set the grid environment *********************************************
setcell %sl_elev%
setwindow %sl_elev%
&describe %sl_elev%
/* Create a depressionless DEM ******************************************
&run fil.aml %sl_elev% sl_DEM
/* Create an outflow direction grid **************************************
sl_outflow = flowdirection(sl_DEM)
/* Create a possible inflow grid ******************************************
sl_inflow = focalflow(sl_DEM)
/* Calculate the degree of the down slope for each cell ******************
&run dn_slope.aml sl_DEM sl_slope
/* Calculate the slope length for each cell *******************************
&sv cell_size = [show scalar $$cellsize]
&sv diagonal_length = 1.414216 * %cell_size%
/* Convert to radians for cos--to calculate slope length
if (sl_outflow in {2, 8, 32, 128})
sl_length = %diagonal_length% div cos(sl_slope div deg)
else
sl_length = %cell_size% div cos(sl_slope div deg)
endif
/* Set the window with a one cell buffer to avoid NODATA around the edges **
setwindow [calc [show scalar $$wx0] - [show scalar $$cellsize]] ~
[calc [show scalar $$wy0] - [show scalar $$cellsize]] ~
[calc [show scalar $$wx1] + [show scalar $$cellsize]] ~
[calc [show scalar $$wy1] + [show scalar $$cellsize]]
/* Create a new flow direction grid with a one cell buffer ***************
sl_flow = sl_outflow
kill sl_outflow
sl_outflow = con(isnull(sl_flow), 0, sl_flow)
kill sl_flow
/* Create a grid of the high points and NODATA *************************
/* The high points will have 1/2 their cell slope length for VALUE ***
&run high_pts.aml
/* Create a grid of high points and 0's **********************************
/* This will be added back in after the slope lengths for all other ****
/* cells has been determined for each iteration ***********************
sl_high_pts = con(isnull(sl_cum_l), 0, sl_cum_l)
/* Calculate the cumulative slope length for every cell ******************
&run s_length.aml
/* Calculate the LS value for the soil loss equation **********************
&run ls.aml
/* Reset window and mask **********************************************
setwindow %sl_elev%
setmask %sl_elev%
setmask off
/* Kill temporary grids **************************************************
kill sl_(!DEM outflow inflow length high_pts!)
/* Set the output grid names to the user input **************************
rename sl_cum_l %sl_out%
rename ls_amount %LS_out%
&if [null %sl_angle%] &then
kill sl_slope
&else
rename sl_slope %sl_angle%
&return
fil.aml /* fil.aml **************************************************************** /* The input : a grid of elevations /* The output: a depressionless elevation grid. /* Usage: fil&args DEM_grid fil_DEM /* Copy original elevation grid ******************************************* %fil_DEM% = %DEM_grid% /* Create a depressionless DEM grid ************************************* finished = scalar(0) &do &until [show scalar finished] eq 1 finished = scalar(1) rename %fil_DEM% old_DEM if (focalflow(old_DEM) eq 255) { %fil_DEM% = focalmin(old_DEM, annulus, 1, 1) test_grid = 0 } else { %fil_DEM% = old_DEM test_grid = 1 } endif kill old_DEM /* Test for no more sinks filled ************************************ docell finished {= test_grid end kill test_grid &end &return
dn_slope.aml /* dn_slope.aml ********************************************************* /* The input : a grid of elevations with no depressions /* The output: a grid of down slopes in degrees. /* Usage: dn_slope&args DEM_grid down_slope /* Compute the outflow direction for each cell *************************** dn_outflow = flowdirection(%DEM_grid%) /* Set the window with a one cell buffer ******************************** &describe %DEM_grid% setwindow [calc [show scalar $$wx0] - [show scalar $$cellsize]] ~ [calc [show scalar $$wy0] - [show scalar $$cellsize]] ~ [calc [show scalar $$wx1] + [show scalar $$cellsize]] ~ [calc [show scalar $$wy1] + [show scalar $$cellsize]] /* Create a DEM with a one cell buffer *********************************** /* This prevents NODATA being assigned to the edge cells that flow /* off the DEM. Cells that flow off the DEM will get 0 slope ************ dn_buff_DEM = con(isnull(%DEM_grid%), focalmin(%DEM_grid%), %DEM_grid%) /* Calculate the down slope in degrees ********************************** &sv cell = [show scalar $$cellsize] /* The () pervent problems that occur with using whole numbers **** &sv cell_size = (1.00 * %cell%) &sv diagonal_length = (1.414216 * %cell_size%) /* find down slope cell and calculate slope *************************** if (dn_outflow eq 64) %down_slope% = deg * atan((dn_buff_DEM - dn_buff_DEM(0, -1)) div ~ %cell_size%) else if (dn_outflow eq 128) %down_slope% = deg * atan((dn_buff_DEM - dn_buff_DEM(1, -1)) div ~ %diagonal_length%) else if (dn_outflow eq 1) %down_slope% = deg * atan((dn_buff_DEM - dn_buff_DEM(1, 0)) div ~ %cell_size%) else if (dn_outflow eq 2) %down_slope% = deg * atan((dn_buff_DEM - dn_buff_DEM(1, 1)) div ~ %diagonal_length%) else if (dn_outflow eq 4) %down_slope% = deg * atan((dn_buff_DEM - dn_buff_DEM(0, 1)) div ~ %cell_size%) else if (dn_outflow eq 8) %down_slope% = deg * atan((dn_buff_DEM - dn_buff_DEM(-1, 1)) div ~ %diagonal_length%) else if (dn_outflow eq 16) %down_slope% = deg * atan((dn_buff_DEM - dn_buff_DEM(-1, 0)) div ~ %cell_size%) else if (dn_outflow eq 32) %down_slope% = deg * atan((dn_buff_DEM - dn_buff_DEM(-1, -1)) div ~ %diagonal_length%) else %down_slope% = 0.00 endif /* Reset old settings ***************************************************** setwindow %DEM_grid% /* Clip the output grid *************************************************** dn_slope = %down_slope% kill %down_slope% rename dn_slope %down_slope% /* Kill the temporary grids ********************************************** kill dn_buff_DEM kill dn_outflow &return
high_pts.aml
/* high_pts.aml **********************************************************
/* This is not a stand alone AML ****************************************
/* Grids used from sl.aml:
/* sl_outflow
/* sl_inflow
/* sl_slope
/* Grid produced for sl.aml:
/* sl_cum_l
/* Find the high points and set value to half their own slope length *****
/* A high point is a cell that has no points flowing into it or if the only
/* cells flowing in to it are of equal elevation. ***************************
if ((sl_inflow && 64) and (sl_outflow(0, -1) eq 4))
sl_cum_l = setnull(1 eq 1)
else if ((sl_inflow && 128) and (sl_outflow(1, -1) eq 8))
sl_cum_l = setnull(1 eq 1)
else if ((sl_inflow && 1) and (sl_outflow(1, 0) eq 16))
sl_cum_l = setnull(1 eq 1)
else if ((sl_inflow && 2) and (sl_outflow(1, 1) eq 32))
sl_cum_l = setnull(1 eq 1)
else if ((sl_inflow && 4) and (sl_outflow(0, 1) eq 64))
sl_cum_l = setnull(1 eq 1)
else if ((sl_inflow && 8) and (sl_outflow(-1, 1) eq 128))
sl_cum_l = setnull(1 eq 1)
else if ((sl_inflow && 16) and (sl_outflow(-1, 0) eq 1))
sl_cum_l = setnull(1 eq 1)
else if ((sl_inflow && 32) and (sl_outflow(-1, -1) eq 2))
sl_cum_l = setnull(1 eq 1)
/* Flat high points get 0 instead of 1/2 slope length ******************
else if (sl_slope eq 0)
sl_cum_l = 0.0
else
sl_cum_l = 0.5 * sl_length
endif
&return
s_length.aml
/* s_length.aml **********************************************************
/* This is not a stand alone AML ****************************************
/* Grids used from sl.aml:
/* sl_inflow
/* sl_outflow
/* sl_slope
/* sl_length
/* sl_high_pts
/* sl_DEM
/* Grid produced for sl_aml:
/* sl_cum_l
/* Pervents the testing of the buffer cells *******************************
setmask sl_DEM
/* Calculate the cumulative slope length for each cell *******************
nodata_cell = scalar(1)
&sv finished = .FALSE.
&do &until %finished%
rename sl_cum_l sl_out_old
&sv counter = 0
&do counter = 1 &to 8
/* Set the varibles for the if that follows
&select %counter%
&when 1
&do
&sv from_cell_grid = sl_north_cell
&sv from_cell_direction = 4
&sv possible_cell_direction = 64
&sv column = 0
&sv row = -1
&end
&when 2
&do
&sv from_cell_grid = sl_NE_cell
&sv from_cell_direction = 8
&sv possible_cell_direction = 128
&sv column = 1
&sv row = -1
&end
&when 3
&do
&sv from_cell_grid = sl_east_cell
&sv from_cell_direction = 16
&sv possible_cell_direction = 1
&sv column = 1
&sv row = 0
&end
&when 4
&do
&sv from_cell_grid = sl_SE_cell
&sv from_cell_direction = 32
&sv possible_cell_direction = 2
&sv column = 1
&sv row = 1
&end
&when 5
&do
&sv from_cell_grid = sl_south_cell
&sv from_cell_direction = 64
&sv possible_cell_direction = 4
&sv column = 0
&sv row = 1
&end
&when 6
&do
&sv from_cell_grid = sl_SW_cell
&sv from_cell_direction = 128
&sv possible_cell_direction = 8
&sv column = -1
&sv row = 1
&end
&when 7
&do
&sv from_cell_grid = sl_west_cell
&sv from_cell_direction = 1
&sv possible_cell_direction = 16
&sv column = -1
&sv row = 0
&end
&when 8
&do
&sv from_cell_grid = sl_NW_cell
&sv from_cell_direction = 2
&sv possible_cell_direction = 32
&sv column = -1
&sv row = -1
&end
&end
/* Test for possible flow source cell
if (not(sl_inflow && %possible_cell_direction%))
%from_cell_grid% = 0
/* Test for flow source cell
else if (sl_outflow(%column%, %row%) <> %from_cell_direction%)
%from_cell_grid% = 0
/* Test flow source cell for nodata
else if (isnull(sl_out_old(%column%, %row%)))
%from_cell_grid% = setnull(1 eq 1)
/* Test current cell slope against cutoff value
else if (sl_slope >= (sl_slope(%column%, %row%) * %.slope_cutoff_value%))
%from_cell_grid% = sl_out_old(%column%, %row%) + ~
sl_length(%column%, %row%)
else
%from_cell_grid% = 0
endif
&end
/* Select the longest slope length
sl_cum_l = max(sl_north_cell, sl_NE_cell, sl_east_cell, sl_SE_cell, ~
sl_south_cell, sl_SW_cell, sl_west_cell, sl_NW_cell, ~
sl_high_pts)
/* Kill the temporary grids
kill (!sl_north_cell sl_NE_cell sl_east_cell sl_SE_cell ~
sl_south_cell sl_SW_cell sl_west_cell sl_NW_cell!)
kill sl_out_old
/* Test for the last iteration filling in all cells with data
&sv no_data = [show scalar nodata_cell]
&if %no_data% eq 0 &then
&sv finished = .TRUE.
/* Test for any nodata cells
if (isnull(sl_cum_l) and not isnull(sl_outflow))
sl_nodata = 1
else
sl_nodata = 0
endif
nodata_cell = scalar(0)
docell
nodata_cell }= sl_nodata
end
kill sl_nodata
&end
/* Reset original window and clip the cumulative slope length grid *****
setwindow sl_DEM
rename sl_cum_l sl_out_old
sl_cum_l = sl_out_old
kill sl_out_old
&return
/* ls.aml *****************************************************************
/* This is not a stand alone AML ****************************************
/* Grids used from sl.aml:
/* sl_cum_l
/* sl_slope
/* Grid produced for sl.aml:
/* ls_amount
/* Convert meters to feet if necessary ***********************************
&if %.grid_units% eq METERS &then
ls_length = sl_cum_l div 0.3048
&else
ls_length = sl_cum_l
/* Calculate LS for the soil loss equation ********************************
/* For cells of depostion
if (ls_length eq 0)
ls_amount = 0
/* For slopes 5% and over
else if (sl_slope >= 2.862405)
ls_amount = pow((ls_length div 72.6), 0.5) * ~
(65.41 * pow(sin(sl_slope div deg), 2) + ~
4.56 * sin(sl_slope div deg) + 0.065)
/* For slopes 3% to less than 5%
else if ((sl_slope >= 1.718358) and (sl_slope < 2.862405))
ls_amount = pow((ls_length div 72.6), 0.4) * ~
(65.41 * pow(sin(sl_slope div deg), 2) + ~
4.56 * sin(sl_slope div deg) + 0.065)
/* For slopes 1% to less than 3%
else if ((sl_slope >= 0.572939) and (sl_slope < 1.718358))
ls_amount = pow((ls_length div 72.6), 0.3) * ~
(65.41 * pow(sin(sl_slope div deg), 2) + ~
4.56 * sin(sl_slope div deg) + 0.065)
/* For slopes under 1%
else
ls_amount = pow((ls_length div 72.6), 0.2) * ~
(65.41 * pow(sin(sl_slope div deg), 2) + ~
4.56 * sin(sl_slope div deg) + 0.065)
endif
/* kill temporary grids **************************************************
kill ls_length
&return
Return to Bob's slope page.