Cumulative Downhill Slope Length AMLs

(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.