/*--------------------------------------------------------------------------------------
/*  USDA Forest Service - Rocky Mountain Research Station - FSL. Moscow, ID
/*--------------------------------------------------------------------------------------
/*  Program: CTI
/*  Purpose: Caculates Compound Topographic Index
/*
/*--------------------------------------------------------------------------------------
/*  Usage: CTI <DEM> <OUTGRID> {CONTRIBUTING AREA}
/*
/*  Arguments: DEM - Digital Elevtion Model
/*             OUTGRID - Output name of final CTI grid
/*             CONTRIBUTING AREA - Optional contributing area Grid
/*
/*--------------------------------------------------------------------------------------
/*  Notes: CTI is a steady state wetness index. The CTI is a function of both the
/*         slope and the upstream contributing area per unit width orthigonal to
/*         the flow direction. CTI was desigined for hillslope catenas. Accumulation
/*         numbers in flat areas will be very large and CTI will not be a relevant
/*         variable. CTI is highly correlated with several soil attributes such as
/*         horizon depth(r=0.55), silt percentage(r=0.61), organic matter content(r=0.57),
/*         and phosphorus(r=0.53) (Moore et al. 1993).
/*
/*                       The implementation of CTI can be shown as:
/*
/*                              CTI = ln (As / (tan (beta))
/*
/*                        where As = Area Value calculated as
/*                    (flow accumulation + 1 ) * (pixel area in m2)
/*                     and beta is the slope expressed in radians.
/*
/*
/*         The ArcInfo approach to calculating Flow Direction uses the D8 algorithm
/*         producing very unrealistic results. Several other methods are available
/*         for calculating flow directions. One of the more robust approaches is the
/*         D infinity algorithm (Tarboton 1997). There is a freeware download
/*         and documentation for TARDEM at http://www.engineering.usu.edu/dtarb/
/*         or TAUDEM http://moose.cee.usu.edu/taudem/taudem.html
/*         The deritive of the CONTRIBUTING AREA caculation from either of these two programs
/*         can be used in the CTI AML or you can caculate CONTRIBUTING AREA in GRID.
/*
/*         NOTE: If you use CTI without the contributing area input all caculations
/*               will be conducted in ArcInfo - GRID using the D8 algorithm.
/*
/*--------------------------------------------------------------------------------------
/*  Required Input: <Digital Elevation Model> <OUT GRID>
/*  Optional Input  (contributing area)
/*  Output:         Compound Topographic Index (CTI)
/*
/*--------------------------------------------------------------------------------------
/*  History: Jeffrey Evans - 12/10/01 - Original coding
/*           1221 South Main, Moscow, ID 83843
/*           (208) 882-1221
/*
/*           05/06/2003
/*           Changed coding; to accept optional contributing area grid
/*                           to recognize the analysis cell size
/*
/*           Set up to run as a ATOOL
/*
/*======================================================================================
/*  References:
/*
/*    Gessler, P.E., I.D. Moore, N.J. McKenzie, and P.J. Ryan. 1995.
/*    Soil-landscape modeling and spatial prediction of soil attributes.
/*    International Journal of GIS. Vol 9, No 4, 421-432.
/*
/*    Moore, ID., Gessler, P.E., Nielsen, G.A., and Petersen, G.A. 1993
/*    Terrain attributes: estimation methods and scale effects.
/*    In Modeling Change in Environmental Systems, edited by A.J. Jakeman
/*    M.B. Beck and M. McAleer (London: Wiley), pp. 189 - 214.
/*
/*    Tarboton, D. G., (1997), A New Method for the Determination of Flow
/*    Directions and Contributing Areas in Grid Digital Elevation Models,
/*    Water Resources Research, 33(2): 309-319.
/*
/*======================================================================================
/*
&args dem outgrid ca

&if [show PROGRAM] <> GRID &then
  &do
     grid
       &type Can Only Be run From GRID
          &type Starting GRID
             &type Please re-run CTI
     &end

&if [NULL %dem%] = .TRUE. &then
  &return &inform Usage:CTI <DEM> <OUTGRID> {CONTRIBUTING AREA}

&if [NULL %outgrid%] = .TRUE. &then
  &return &inform Usage:CTI <DEM> <OUTGRID> {CONTRIBUTING AREA}

&if [exists %dem% -grid] = .FALSE. &then
  &return &inform Grid [upcase %dem%] does not exist!

&if [exists %outgrid% -grid] = .TRUE. &then
  &return &inform Grid [upcase %outgrid%] already exist!

&if [NULL %ca%] = .TRUE. &then
  &do

setcell %dem%
setwindow %dem%

&s tmp1 [scratchname -prefix xx1] /* dem pit fill
&s tmp2 [scratchname -prefix xx2] /* flow direction
&s tmp3 [scratchname -prefix xx3] /* flow accumulation
&s tmp4 [scratchname -prefix xx4] /* slope degrees
&s tmp5 [scratchname -prefix xx5] /* slope radians conversion
&s tmp6 [scratchname -prefix xx6] /* slope tan
&s tmp7 [scratchname -prefix xx7] /* slope reclass
&s tmp8 [scratchname -prefix xx8] /* area value

&messages &off

describe %dem%

&s cellsize = %GRD$DX%

&messages &on

&type /& Pit filling DEM and Caculating Flow Directions /&
fill %dem% %tmp1% sink 50 %tmp2%

&type /& Caculating Flow Accumulation /&
%tmp3% = flowaccumulation(%tmp2%) 

&type /& Caculating Slope in Degrees /&
%tmp4% = slope(%tmp1%, degree)

&type /& Converting slope degrees to radians using: slope * (pi / 2) / 90 /&
%tmp5% = (%tmp4% * 1.570796) / 90

&type /& Caculating tangent of slope(radians) /&
%tmp6% = tan(%tmp5%)

&type /& Reclassing 0 slope values to .001 /&
%tmp7% = con(%tmp6% == 0, .001, %tmp6%)

&type /& Caculating upslope surface contributing area /&
%tmp8% = (%tmp3% + 1) * %cellsize%

&type /& Caculating Compound Topographic Index /&
%outgrid% = ln((%tmp8% / %tmp7%))

&type /& Cleaning Up Temporary GRIDS /&

&messages &off

kill (!%tmp1% %tmp2% %tmp3% %tmp4% %tmp5% %tmp6% %tmp7% %tmp8%!) all

&messages &on

&type /& Compound Topographic Index written to [upcase %outgrid%] /&

&end

&else &do

&if [exists %ca% -grid] = .FALSE. &then
  &return &inform Grid [upcase %ca%] does not exist!

&messages &off

describe %ca%

&messages &on

&s cellsize = %GRD$DX%

setcell %cellsize%
setwindow %ca%

&s tmp1 [scratchname -prefix xx1] /* Slope Radians
&s tmp2 [scratchname -prefix xx2] /* Slope Reclass
&s tmp3 [scratchname -prefix xx3] /* Slope Tangent
&s tmp4 [scratchname -prefix xx4] /* Area Value

&type /& Converting slope degrees to radians using: slope * (pi / 2) / 90 /&
%tmp1% = slope(%dem%) * 1.570796 / 90

&type /& Reclassing 0 slope values to .001 /&
%tmp2% = con(%tmp1% == 0, .001, %tmp1%)

&type /& Caculating tangent of slope(radians) /&
%tmp3% = tan(%tmp2%)

&type /& Caculating upslope surface contributing area /&
%tmp4% = (%ca% + 1) * %cellsize%

&type /& Caculating Compound Topographic Index /&
%outgrid% = ln((%tmp4% / %tmp3%))

&type /& Cleaning Up Temporary GRIDS /&

&messages &off

kill (!%tmp1% %tmp2% %tmp3% %tmp4%!) all

&messages &on

&type /& Compound Topographic Index written to [upcase %outgrid%] /&

&end

