## Knickzone Extraction Tool (KET) Part1: Calculation of relative steepness
## Importing modules
import arcpy
import os
from arcpy.sa import *

arcpy.CheckOutExtension("Spatial")
arcpy.env.overwriteOutput = True

##Input Variables
##the input for workspace
work = arcpy.env.workspace = arcpy.GetParameterAsText(0)

##The input DEM (must be an uncompressed TIFF)
inZfile = arcpy.GetParameterAsText(1)
desc = arcpy.Describe(inZfile)
DEM = str(desc.catalogPath)


##The range and step value for averaging distance 'd'
FromValue = int(arcpy.GetParameterAsText(2)) ##the starting mimimum
ToValue = int(arcpy.GetParameterAsText(3)) ##the limiting maximum
ByValue = int(arcpy.GetParameterAsText(4)) ##the stepping value

## Local variables:
Distance = FromValue
SumRaster = 0
SumAllProduct = 0
num = "8"
FinalRaster = "Finalraster.tif"
DEMfill = os.path.join(work+"€€DEMfill"+".tif") ##Hydrologically corrected DEM

## Process: Pit Remove (MPIEXEC compilation of TauDEM is used in this Step)
cmd = 'mpiexec -n ' + num + ' pitremove -z ' + '"' + DEM + '"' + ' -fel ' + '"' + DEMfill + '"'
os.system(cmd)

desc2 = arcpy.Describe(DEMfill)
fel = str(desc2.catalogPath)

p = os.path.join(work+"€€flowdr"+".tif")  ##declaring flow direction variable
sd8 = os.path.join(work+"€€flslop"+".tif") ##declaring slope variable

## Process: D8 Flow Directions

cmd2 = 'mpiexec -n ' + num + ' D8FlowDir -fel ' + '"' + fel + '"' + ' -p ' + '"' + p + '"' + ' -sd8 ' + '"' + sd8 + '"'
os.system(cmd2)

desc3 = arcpy.Describe(p)
p = str(desc3.catalogPath)

##Loop: looping through the incremental averaging distance 'd'

count = (ToValue- FromValue)/ByValue + 1

for i in range(count):
    n = i + 1
    distance = str(Distance)
    slpd = os.path.join(work+"€€slpd"+distance+".tif")

    ## Process: Slope Average Down
    cmd3 = 'mpiexec -n ' + num + ' SlopeAveDown -p ' + '"' + p + '"' + ' -fel ' + '"' + fel + '"' + ' -slpd ' + '"' + slpd + '"' + ' -dn ' + distance
    os.system(cmd3)

    Distance = Distance + ByValue ##the incremental distance
    arcpy.AddMessage("€nCurrently processing:€n" + distance)
    SumRaster = SumRaster + Raster(slpd) ##for the summation of all y(n)
    RasterProduct = n * Raster(slpd) ## for the product of y(n) * n
    SumAllProduct = SumAllProduct + RasterProduct ##for the summation of all y(n) * n
    i = i + 1

## Finalising the calculations
RasterY = ((12 * SumAllProduct) - (6 * (count + 1) * SumRaster)) / (ByValue * count * (count + 1) * (count - 1))
##Saving the final raster
RasterY.save(FinalRaster)