/*----------------------------------------------------------------------------
/*
/* NAME: lsfactor.aml
/*
/*-------------------------------------------------------------------------------------
/* AUTHOR: Jacek S. Blaszczynski, Physical Scientist. BLM NARSC, Original coding, 04/99
/*         Revision and additions of LSFACTOR2, 02/2000
/*P----------------------------- Purpose ----------------------------------------------
/*
/* To calculate grids of values for the terrain factor of the Revised 
/* Universal Soil Loss Equation from digital elevation model data using the
/* empirically based LS factor equation described in the BLM RUSLE diskette 
/* (Simanton, 1987) and physical-process based LS equation developed 
/* by Moore and Wilson (1992) and Mitasova (1993).  
/*
/* The results of the statistical analysis of about 10000 plot 
/* years of data was used in developing what is now a well known 
/* model for sheet and rill soil loss prediction called the 
/* Universal Soil Loss Equation.  Soil loss equations are 
/* empirical (based on field collected data) models developed to 
/* enable conservation planners to project limited erosion data 
/* to the many localities and conditions that have never been 
/* directly represented in field research.  The Revised USLE was 
/* developed by the Agricultural Research Service to further extend the 
/* USLE (prepared for croplands) to wild areas of rangelands, particularly in 
/* arid and semi-arid climates.  The RUSLE reflects the basic USLE structure: 
/*
/*                A = R * K * L * S * C * P     
/*
/*                            where
/*
/*   A  is the computed soil loss per unit area, usually in tons per
/*       acre per year;
/*
/*   R  is the rainfall and runoff factor and is the number of rainfall 
/*       erosion index units;
/*
/*   K  is the soil erodibility factor, is the soil loss rate per 
/*       erosion index unit for a specified soil as measured on a unit plot,
/*       which is defined as a 72.6-ft length of uniform 9-percent slope 
/*       continuously in clean tilled fallow;
/*         
/*   L  is the slope-length factor, is the ratio of soil loss from the 
/*       field slope length to that from a 72.6-ft length under identical 
/*       conditions; 
/*                    
/*   S  is the slope-steepness factor, is the ratio of soil loss from 
/*       the field slope gradient to that from a 9-percent slope under 
/*       otherwise identical conditions;
/*
/*   C  is the cover and management factor, is the ratio of soil loss 
/*       from an area with specified cover and management to that from an 
/*       identical area in tilled continuous fallow;
/*      
/*   P  is the support practice factor, is the ratio of soil loss with 
/*       a support practice like contouring, stripcropping, or terracing to 
/*       that with straight-row farming up and down the slope (Wischmeier and 
/*       Smith, 1978).
/*
/* The automated AML capability provided here works in ARC GRID to determine
/* the topographic,  terrain,  or the slope gradient (S)-slope length (L) /* 
/* factors of the RUSLE equation.  The mathematical methods for determination 
/* of the terrain factor used in preparation of the AML were obtained from the 
/* RUSLE diskette prepared for the BLM in 1987.  The AML first prepares two 
/* grids from Digital Elevation Model (DEM): a slope gradient grid, and a 
/* slope length grid, after which the two are mathematically combined to 
/* prepare the LS-factor grid.  Because the equation does not work for slopes 
/* over 800 ft long, all the slope lengths in the slope length grid do not 
/* exceed 800 ft.  After the LS-factor grid has been prepared it can then be 
/* classified into low, medium, and high levels of potential erosion.
/*
/* The LS-factor only provides a method for analysis of the influence of 
/* terrain, while the overall soil loss potential is calculated not only from
/* terrain, but also other data since the RUSLE model is multidisciplinary. 
/* Therefore the terrain factor map provides only partial information on 
/* potential erosion.  The information will be more accurate if soils data are 
/* available which contain values for soil erodibility or K-factor for each of 
/* the soil types (link to soils page), and even more accurate if vegetation 
/* data (link to vegetation cover map) are available which can permit 
/* calculation of the surface cover, or the C-factor for each surface cover 
/* type.  Multiplying grids containing K-factor values per soil type and C-
/* factor values per vegetation type by the LS factor map will result in a 
/* most accurate model of potential erosion, but in absence of these 
/* additional data sets the terrain factor can serve as a tool for at least 
/* partial evaluation of the erosiveness of the landscape.
/*
/*
/*
/*A---------------------------- Arguments ------------------------------------
/*
/* input_dem   raster map of elevations (Digital Elevation Model)
/*
/*H----------------------------- History -------------------------------------
/*
/* AUTHOR:  Jacek Blaszczynski, Physical Scientist, BLM NARSC 
/* DATE:    4/99
/* EVENT:   Original Coding
/* 
/*E===========================================================================
&severity &error &routine bailout 
&args input_dem cell_size 

/* Check arguments
&call chkargs

/* Do the work
&call main

&call exit
&return /* End of LSfactor.aml

/*--------------
&routine chkargs
/*--------------
&if [SHOW PROGRAM] <> GRID &then &do
  &type This aml must be run from GRID.
  &call bailout
&end
&do arg &list input_dem cell_size
  &if [NULL [VALUE %arg%]] &then &do
    &call usage
    &call bailout
  &end
&end
&if not [EXISTS %input_dem% -GRID] &then &do
  &type Input grid: %input_dem%, does not exist.
  &call bailout
&end
&return /* End of routine chkargs
/*-----------
&routine main
/*-----------
&wat lsfactor.wat

/*-----------------------------------------------------------------------------
/* 
/* CLEAN MAPS PRODUCED ON THE PREVIOUS LFACTOR.AML RUN 
/*
/*----------------------------------------------------------------------------

&if [EXISTS mmap -GRID] &then &do  
  kill MMAP all
&end
&if [EXISTS smap -GRID] &then &do 
  kill SMAP all
&end
&if [EXISTS smap2 -GRID] &then &do 
   kill SMAP2 all
&end
&if [EXISTS bmap -GRID] &then &do 
  kill BMAP all
&end
&if [EXISTS slopes0to8 -GRID] &then &do   
    kill SLOPES0TO8 all
&end
&if [EXISTS slopes8to50 -GRID] &then &do 
  kill SLOPES8TO50 all
&end
&if [EXISTS lfactor1 -GRID] &then &do  
  kill LFACTOR1 all
&end
&if [EXISTS sfact0to8 -GRID] &then &do
  kill sfact0to8 all
&end
&if [EXISTS sfact8to50 -GRID] &then &do
    kill sfact8to50 all
&end
&if [EXISTS sfactor1 -GRID] &then &do 
  kill SFACTOR1 all
&end
&if [EXISTS lsfactor1 -GRID] &then &do 
  kill LSFACTOR1 all
&end
&if [EXISTS lsfactor2 -GRID] &then &do 
    kill lsfactor2 all
&end
&if [EXISTS lsfactor3 -GRID] &then &do 
    kill lsfactor3 all
&end
&if [EXISTS lfactor2 -GRID] &then &do 
    kill lfactor2 all
&end
&if [EXISTS sfactor2 -GRID] &then &do 
    kill sfactor2 all
&end
&if [EXISTS lfactor3 -GRID] &then &do 
    kill lfactor3 all
&end
&if [EXISTS sfactor3 -GRID] &then &do 
    kill sfactor3 all
&end
&if [EXISTS sloplenm -GRID] &then &do 
    kill sloplenm all
&end
&if [EXISTS temp3 -GRID] &then &do 
    kill TEMP3 all
&end
&if [EXISTS sloplenft -GRID] &then &do
    kill sloplenft all
&end
&if [EXISTS sloplenft2 -GRID] &then &do
    kill sloplenft2 all
&end
&if [EXISTS flowcounta2 -GRID] &then &do 
    kill flowcounta2 all
&end
&if [EXISTS flowlengtha2 -GRID] &then &do 
    kill flowlengtha2 all
&end
&if [EXISTS flowcountz2 -GRID] &then &do 
    kill flowcountz2 all
&end


/*---------------------------------------------------------------------------
/* Set the environment
/*---------------------------------------------------------------------------

&terminal 9999
display 9999 2 position cr
&echo &off
&echo &brief
setwindow %input_dem% %input_dem%
mape %input_dem%

/*----------------------------------------------------------------------------
/* Calculate LSFACTOR1 based on equations from Simanton, 1987
/*----------------------------------------------------------------------------

&if not [EXISTS flowdira -GRID] &then &do
  &type Preparing a flowdirection map for an unfilled DEM...
  flowdira = flowdirection (%input_dem%)
  buildsta flowdira
&end
gridpaint flowdira value linear nowrap gray

/*------------------------------------------------------------------------
/* CALCULATE THE SLOPE LENGTH FACTOR
/*------------------------------------------------------------------------

&if not [EXISTS flowcountz2 -GRID] &then &do
  &if not [EXISTS flowcountz -GRID] &then &do
     &if not [EXISTS flowdirz -GRID] &then &do
        &if not [EXISTS filled -GRID] &then &do
           fill %input_dem% filled
           buildsta filled
        &end
        flowdirz = flowdirection (filled)
        buildsta flowdirz
     &end
     flowcountz = flowaccumulation (flowdirz)
     buildsta flowcountz
  &end      
/*---------------------------------------------------------------
/* Generate a map consisting only of hillslope areas by setting
/* all filled dem flow accumulation (flowcountz) values above
/* 50 to a null or NODATA value.  By editing this value we can
/* change the extent of what we consider hillslopes.  Another
/* alternative is to use an existing stream coverage to create
/* a FLOWCOUNTZ2 map prior to running this AML (contact author)
/*--------------------------------------------------------------
  flowcountz2 = con (flowcountz < 50, 1, setnull (flowcountz) )
  buildsta flowcountz2
&end

/*-------------------------------------------------------------------
/* Build a flow length map calculating the distance along steepest
/* flow path.  The flow length is calculated for the unfilled DEM
/* so that slope lengths represented are between surface depressions
/* in the terrain. Slope lenghts are calculated only for the areas
/* that have been identified in map FLOWCOUNTZ2 as hillslopes, 
/* i.e. not part of a drainage.
/*------------------------------------------------------------------- 

&if not [EXISTS flowlengtha2 -GRID] &then &do
   &if not [EXISTS flowlengtha -GRID] &then &do
      flowlengtha = flowlength (flowdira, #, upstream)
      buildsta flowlengtha
   &end
   flowlengtha2 = flowlengtha *  ( int ( flowcountz2 ) )
   buildsta flowlengtha2
&end 


gridpaint flowlengtha2 value linear nowrap gray

/*---------------------------------------------------------------
/* This LS-factor equation uses feet and has the limitation of 
/* working only on slope gradients less than 50 percent and
/* slope lengths less than 800 feet.  We maintained use of feet
/* rather than meters as units for slope length to keep it 
/* as much the same as the equation in the BLM RUSLE program 
/* as possible. We renumbered all values greater than 800 feet 
/* to 800.
/*--------------------------------------------------------------

&if not [EXISTS sloplenft -GRID] &then &do
    sloplenft = flowlengtha2 * 3.28
    buildsta sloplenft   
&end

gridpaint sloplenft value linear nowrap gray

&if not [EXISTS sloplenft2 -GRID] &then &do
    sloplenft2 = con(sloplenft > 800, 800, sloplenft)
    buildsta sloplenft2
&end

gridpaint sloplenft2 value linear nowrap gray

&if not [EXISTS slopeperc -GRID] &then &do
  slopeperc = slope (%input_dem%, percentrise)
  buildsta slopeperc
&end

gridpaint slopeperc value linear nowrap gray

SMAP = ( SIN( ATAN ( SLOPEPERC / 100 ) ) )
SMAP2 = POW (SMAP, 0.79)
BMAP =  ( ( ( SMAP / 0.0895999 ) * 0.5 ) / ( (2.96 * SMAP2 ) + 0.56 ) )
kill smap all
kill smap2 all
MMAP = BMAP / (1 + BMAP)
kill bmap all
LFACTOR1 = POW ( ( SLOPLENFT / 72.6 ), MMAP)
/*buildsta lfactor1 
gridpaint lfactor1 value linear nowrap gray
kill mmap all

/*--------------------------------------------------------------------------
/*  calculate the slope gradient factor (s) using two equations: one for
/*  slopes less than 8 percent and the other for slopes 8-50 percent. 
/*  these maps are combined to create an s-factor map in which all areas
/*  over 50 percent is slope gradient are excluded (given a nodata value)
/*--------------------------------------------------------------------------

SLOPES0TO8 = CON (SLOPEPERC > 8, SETNULL (SLOPEPERC), SLOPEPERC) 

if (SLOPEPERC lt 8 or SLOPEPERC gt 50) 
  slopes8to50 = setnull ( SLOPEPERC )
  else 
  slopes8to50 = SLOPEPERC
endif

SFACT0TO8 = ( ( SIN ( ATAN ( SLOPES0TO8 / 100 ) ) ) * 10 ) + 0.027 
SFACT8TO50 =  ( ( SIN ( ATAN ( SLOPES8TO50 / 100 ) ) ) * 17.2 ) - 0.56 
SFACTOR1 = MERGE (SFACT0TO8, SFACT8TO50)
/* buildsta sfactor1
kill sfact0to8 all
kill sfact8to50 all
kill slopes0to8 all
kill slopes8to50 all
gridpaint sfactor1 value linear nowrap gray

/*-------------------------------------------------------------------------
/* Calculate the LS (terrain) factor according to the equation used
/* in the BLM RUSLE program (contact author for the program).  The output
/* units are equivalent to the units calculated using the second method
/* that follows, even though here we are using feet and slope percent and
/* in the second method we are using slope degrees and meters.  Since
/* this equation is limited to calculations within 800 feet slope length
/* and not greater than 50 percent slopes, the LS-factor values here
/* are lower than in the second equation which calculates the factor for
/* all slope gradients and slope lengths.
/*-------------------------------------------------------------------------

LSFACTOR1 = LFACTOR1 * SFACTOR1
/* buildsta lsfactor1

describe lsfactor1
gridpaint lsfactor1 value linear nowrap gray


/*---------------------------------------------------------------------------------
/* Calculate LSFACTOR2 using methods developed by Moore, I. D. and J. P. Wilson, 1992,
/* Length-slope factors for RUSLE: Simplified method of estimation, Soil Science Society
/* of America Journal, vol. 50, pp. 1294-1298,  and Mitasova, H., 1993, Surfaces and
/* modeling, Grassclippings, Spring 1993.  The calculation of the factor is based on
/* physical process analysis and the results are equivalent to LS-factor values
/* calculate through other means.  Since this method calculates the LS-factor for
/* all slope gradients and slope lengths it might produce extremely high values for
/* some locations (compare the output maps produced by this program).
/* This method calculates LS-factor for flow accumulations not greater than 50 cells
/* which is equivalent to a drainage area of about 122 acres.  That means that 
/* accumulation of flow from an area of 122 acres has to happen before we consider
/* these locations as drainages (channels) rather than hillslopes. Therefore the method
/* uses a flow count map to define hillslope erosional areas
/*---------------------------------------------------------------------------------

lfactor2 = POW ( ( flowlengtha2 / 22.13 ) , 1.3 ) 
slopedeg = slope (%input_dem%, 1, degree)
sfactor2 = POW ( ( ( sin (slopedeg DIV DEG) ) / 0.0896 ), 0.6 )
buildsta sfactor2
lsfactor2 = lfactor2 * sfactor2
buildsta lsfactor2
describe lsfactor2
gridpaint lsfactor2 value linear nowrap gray
&type TO BETTER UNDERSTAND RESULTS (NO DATA AREAS) CREATE A STREM MAP WITH 50 FLOW THRESHOLD
&type USING STREAMNETWORK.AML AND DISPLAY IT (OR ANOTHER APPROPRIATE STREAM MAP) OVER THE LSFACTOR MAPS
&wat &off
/*
/*
&return /* End of routine main
 
/*------------
&routine usage
/*------------
&type Usage: &r lsfactor <input_dem> <cell_size> 
&type (Overwrites output maps prepared on the previous run)

&return /* End of routine usage

/*-----------
&routine exit
/*-----------
&do grd &list lsfactor
  &if [VARIABLE %grd%] &then &do
    &if [EXISTS [VALUE %grd%] -GRID] &then &do
      &type LSFACTOR grid is ready.
    &end
  &end
&end
&return /* End of routine exit

/*--------------
&routine bailout
/*--------------
&severity &error &ignore
&call exit
&return &error Bailing out of lsfactor.aml...