"""
 Toolbox Name:Low_relief_finder 
 Tool List:
 Description: 
The purpose of this script is to highlight areas of an Elevation raster that are filled
during the "Fill Sinks". This will essentially produce a DEM of difference. This script will also
output a polygon shapefile of areas that have been filled w/o a representation of depth to guide a
reviewers eyes to areas that could be of interest with engaging in for more careful examination.

INPUTS:
#### INPUTS ####
Unfilled Raster DEM
Filled Raster DEM (Optional if using NHDPlus)

#### OUTPUTS ####
1) Filled Raster
2) Filled Depth Raster
3) Filled Areas Polygon Shapefile


 Author: Matthew C. Morriss
 Created: September 4th, 2020 
 Language: python3
 History: https://pro.arcgis.com/en/pro-app/arcpy/geoprocessing_and_python/a-template-for-python-toolboxes.htm
          https://pro.arcgis.com/en/pro-app/arcpy/geoprocessing_and_python/a-quick-tour-of-python-toolboxes.htm
          
As of October 23,2020, adding an optional input for a filled DEM at the begining and made all the outputs derived not required
"""

# imports
import os
import sys
import arcpy
from arcpy.sa import *

class Fill_polygon(object):
    def __init__(self):
        """"
        Ininitializes a tool class.
        """
        self.label       = "Fill Depth and Polygon"
        self.description = "This tool was written to take an unfilled DEM, fill the sinks" + \
                           "And then calculate the difference between the two (filled-unfilled) " + \
                           "to calculate the overall fill depth for the DEM." +\
                           "Areas highlighted will be 1 standard deviation below the mean relief of the grid."
        self.canRunInBackground = True
        self.callfrom_pyt = True
        self.category = 'Archive'
        
    def getParameterInfo(self):
        """
        Define the parameters for use in arcmap/pro.
        This function is only needed if you are going to run the 
        script as a tool.
        """
        ##### 1) parameter: INPUT RASTER
        in_features = arcpy.Parameter(
            displayName="Input Raster",
            name="in_features",
            datatype="GPRasterLayer",  # data types: https://pro.arcgis.com/en/pro-app/arcpy/geoprocessing_and_python/defining-parameter-data-types-in-a-python-toolbox.htm
            parameterType="Required",  
            direction="Input")
        
        ##### 0) Paremter: Boolean checkbox if providing a filled raster
        in_filled_bool = arcpy.Parameter(
            displayName = "Providing Filled Raster?",
            name = 'in_filled_bool',
            datatype = 'GPBoolean',
            parameterType = 'Optional',
            direction = 'Input')
        in_filled_bool.value = 'False'
        
        ##### A) Parameter: Input Filled Raster (OPTIONAL)) only if using NHD plus or if user already created 
        in_fill_raster = arcpy.Parameter(
            displayName = "Filled Raster",
            name = "in_fill_raster",
            datatype = "GPRasterLayer",
            parameterType = "Optional")
        
        
        #### 2) Second parameter, output raster depth raster name/location
        fName_raster = arcpy.Parameter(
            displayName="Output Filled Raster",
            name="Fill_raster",
            datatype="DERasterDataset",   # see dataytpes above
            parameterType="Required",
            direction="Output")
        #optional_param.value = "default value" 
        
        
        #### 3) Third parameter, output raster depth raster name/location
        Depth_raster = arcpy.Parameter(
            displayName="Output Fill-Depth Raster 1",
            name="Depth_raster",
            datatype="DERasterDataset",   # see dataytpes above
            parameterType="Required",
            direction="Output")
        #optional_param.value = "default value" 
        
        
        #### 3) Fourth parameter, output raster name/location
        fName_polygon = arcpy.Parameter(
            displayName="Output Polygon",
            name="polygon_name",
            datatype="DEFeatureClass",   # see dataytpes above
            parameterType="Required",
            direction="Output")
        #optional_param.value = "default value" 
        
        
        parameters = [in_features,in_filled_bool, in_fill_raster,fName_raster, Depth_raster, fName_polygon]
        
        
        return parameters

    def myfunc1(self, emtpy1, parameters):
        """
        user defined function.
        Can be renamed as desired and inputs can be changed.
        """
        # if arcpy.Exists(feature_layer):
        #     msg = arcpy.Describe(feature_layer).spatialReference.name
        #     arcpy.AddMessage(msg)
        #     print(msg)
        if self.callfrom_pyt:
            in_raster = parameters[0].valueAsText
            ischecked = parameters[1].valueAsText
        else:
            in_raster = parameters[0]
            ischecked = parameters[1]
        arcpy.AddMessage(str(ischecked))
        
        #check if no filled DEM was provided
        if (ischecked) == 'false':
            #If true then proceed through logic normally
            arcpy.AddMessage(str(ischecked))
            if self.callfrom_pyt:
                fName_raster =  parameters[3].valueAsText
                Depth_raster =  parameters[4].valueAsText
                fName_polygon = parameters[5].valueAsText
            else:
                fName_raster =  parameters[3]
                Depth_raster =  parameters[4]
                fName_polygon = parameters[5]
                
            
            #moving onto filling and differencing
            #STEP 1) Fill Raster
            print('filling raster')
            arcpy.AddMessage("Filling Raster")
            OutFill = Fill(in_raster)                   
            
            #STEP 2) Create DEM of Difference with fill and original Rasters
            print('calculating DEM of difference')
            arcpy.AddMessage("Calculating DEM of difference")
            DEM_Diff = OutFill - Raster(in_raster)
            
            
            #### CREATING OUTPUTS ####
            #STEP 3) Remove unfilled areas by setting to null
            threshHeld_DEM = SetNull(DEM_Diff, DEM_Diff, "VALUE = 0")
    
            
            print('Saving Output Polygon')
            arcpy.AddMessage("Saving Output")
            filledShape = arcpy.RasterToPolygon_conversion(Int(threshHeld_DEM), fName_polygon,'NO_SIMPLIFY')
            threshHeld_DEM.save(Depth_raster)
            OutFill.save(fName_raster)
            
            #add fieldand caculate area for this shape
            filledShape = arcpy.AddField_management(filledShape,"area","Double")
            expression1 = ("!SHAPE.area@SQUAREKILOMETERS!")
            filledShape = arcpy.CalculateField_management(filledShape, "area", expression1, "PYTHON")
            
            threshHeld_DEM.save(Depth_raster)
            fName_polygon = "1_" + fName_polygon
            
            # outLocation = r'D:\GIS_Project\Fill_polygon\Test_area\'
            # outFeatureClass = "test_polygon"
            # arcpy.FeatureClassToFeatureClass_conversion(filledShape, outLocation, 
            #                                 outFeatureClass)
            arcpy.AddMessage(fName_polygon)
            # filledShape.save(fName_polygon)
            
        elif str(ischecked == 'true'):
            # arcpy.AddMessage(parameters[1].valueAsText)
            arcpy.AddMessage('Skipping Fill Raster')
            if self.callfrom_pyt:
                OutFill = parameters[2].valueAsText
                fName_raster =  parameters[3].valueAsText
                Depth_raster =  parameters[4].valueAsText
                fName_polygon = parameters[5].valueAsText
            else:
                OutFill = parameters[2]
                fName_raster =  parameters[3]
                Depth_raster =  parameters[4]
                fName_polygon = parameters[5]
                
            
        
            #STEP 1) Fill Raster already provided, go straight to differencing
            print('calculating DEM of difference')
            arcpy.AddMessage("Calculating DEM of difference")
            DEM_Diff = OutFill - Raster(in_raster)
            
            
            #### CREATING OUTPUTS ####
            #STEP 3) Remove unfilled areas by setting to null
            threshHeld_DEM = SetNull(DEM_Diff, DEM_Diff, "VALUE = 0")
    
            
            print('Saving Output Polygon')
            arcpy.AddMessage("Saving Output")
            filledShape = arcpy.RasterToPolygon_conversion(Int(threshHeld_DEM), fName_polygon,'NO_SIMPLIFY')
            filledShape = arcpy.AddField_management(filledShape,"area","Double")
            expression1 = ("!SHAPE.area@SQUAREKILOMETERS!")
            filledShape = arcpy.CalculateField_management(filledShape, "area", expression1, "PYTHON")
            
            threshHeld_DEM.save(Depth_raster)
            fName_polygon = "1_" + fName_polygon
            
            # outLocation = r'D:\GIS_Project\Fill_polygon\Test_area\'
            # outFeatureClass = "test_polygon"
            # arcpy.FeatureClassToFeatureClass_conversion(filledShape, outLocation, 
            #                                 outFeatureClass)
            arcpy.AddMessage(fName_polygon)
            # filledShape.save(fName_polygon)
            

  

    def execute(self, parameters, messages):
        """
        This function must be called "execute" to run in a toolbox.
        This function should reference the above user defined functions and not contain a ton of code.
        """
        # # use validated parameters if ran as a tool
        # if os.path.basename(sys.executable).lower() != 'python.exe':
        #     parameters = [p.valueAsText for p in parameters]
        # self.myfunc1(parameters[0])
        # self.myfunc2(parameters[0], parameters[1])
        
        self.myfunc1(self, parameters)
      

if __name__ == "__main__":
    """ 
    If this script is run on its own the followiung code is executed.
    This block runs the execute function from the Toolname1 class above.
    If possible, environmental settings should be changed in this code block not in 
    the above functions.  
    """
    
    in_features = r'D:\GIS_Project\Fill_polygon\USGS_13_n31w098.tif'
    # in_features = r'C:\GIS_Project\Low_Relief_Script\Test_area\USGS_13_n31w098.tif'
    in_fill_raster = ''
    Depth_raster = r'D:\GIS_Project\Fill_polygon\Test_area\fill_depth.tif'
    fName_raster = r'D:\GIS_Project\Fill_polygon\Test_area\filled_raster.tif'
    fName_polygon = r'D:\GIS_Project\Fill_polygon\Test_area\test.shp' 
    
    parameters = [in_features, in_fill_raster,fName_raster,Depth_raster, fName_polygon]
    tool1 = Fill_polygon()
    tool1.callfrom_pyt = False
    tool1.execute(parameters, None)

    
    
    
    

    

        

