"""
 Toolbox Name:Low_relief_finder 
 Tool List:
 Description: This tool was written to identify areas of low-relief within a raster of elevation data. The goal is to create a polygon layer that highlights these areas at a given threshold (set as 1-sigma below the mean relief). I hope to create a user input for the window of relief, but we'll see how that works.'
 Author: Matthew C. Morriss
 Created:September 3rd, 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
          
This tool was updated on 12/02/2020 to include an interative run. Meaning that if the user so chose, they could click a button that would run the low relief tool - again on areas identified in the first pass of the tool. 
"""

# imports
import arcpy
from arcpy.sa import *
from utils import isolate_env

class Low_relief_finder(object):
    def __init__(self):
        """"
        Ininitializes a tool class.
        """
        self.label       = "G5 Low Relief Finder"
        self.description = ("This tool was written to find areas of low relief and "
                           "highlight these areas by circling them in a new polygon.\n"
                           "Areas highlighted will be 1 standard deviation below the mean relief of the grid. "
                           "This will highlight the lowest 15% of relief.\n"
                           "The User can provide an additional threshold (in terms of standard deviations).")
        self.canRunInBackground = False
        self.callfrom_pyt = True
        self.category = 'Geometry checks'
        
    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.
        """
        ##### 1st 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")
        
        
        #### 2) Second parameter, moving window radius (recommended 1 km)
        wRadius = arcpy.Parameter(
            # value = 1000,
            displayName="Moving Window Radius (m)",
            name="win_radius",
            datatype="GPDouble",   # see dataytpes above
            parameterType="Required",
            direction="Input")
        # optional_param.value = "default value" 
        #parameters = [in_features, optional_param]
        wRadius.value = 1000
        
        #### 3) Check box if providing a lower (custom) thresold than default)
        in_filled_bool = arcpy.Parameter(
            displayName = "Providing Custom Threshold",
            name = 'in_filled_bool',
            datatype = 'GPBoolean',
            parameterType = 'Optional',
            direction = 'Input')
        in_filled_bool.value = 'False'
        
        #### 4) optional third parameter, lowering the threshold for low relief by standard deviations below the mean
        cutoffThreshold = arcpy.Parameter(
            displayName = "Cutoff Threshold (standard deviations)",
            name = 'cutoffThreshold',
            datatype = "GPDouble",
            parameterType = "Optional",
            direction = "Input")
        
        #### 5) optional fourth parameter, run iteratively will add a second pass of the low-relief tool, looking to extract those areas within the region identified in the first run of the tool
        runIter_bool = arcpy.Parameter(
            displayName = "Run Second Iteration?",
            name = 'runIter_bool',
            datatype = 'GPBoolean',
            parameterType = 'Optional',
            direction = 'Input')
        runIter_bool.value  = "false" #set condition unchecked to start
        
        #### 6) Third parameter, output raster name/location
        fName_raster = arcpy.Parameter(
            displayName="Output Raster",
            name="raster_name",
            datatype="DERasterDataset",   # see dataytpes above
            parameterType="Required",
            direction="Output")
        #optional_param.value = "default value" 
        
        
        #### 7) 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,
                      wRadius,
                      in_filled_bool,
                      cutoffThreshold,
                      runIter_bool,
                      fName_raster,
                      fName_polygon]
        
        
        return parameters

    @isolate_env
    def myfunc1(self, empty1, parameters):
        """
        This function will follow the steps below after a raster dataset has been loaded
        #### STEPS ####
        1) Reproject raster to have meters as the linear unit (North American Albers Equal Area Conic)
        2) Using focal statistics calculate elevation range in moving window (standard is 1 km radius)
        3) From this resulting relief raster, find areas of relief 1 standard deviation below the mean
        4) Convert these areas of the raster to polygons and export.
        """
        arcpy.env.parallelProcessingFactor = "100%"

        if self.callfrom_pyt:
            arcpy.AddMessage('Running from Arcmap/Arcpro')
            in_raster = parameters[0].valueAsText
            wRadius = parameters[1].valueAsText
            ischecked = parameters[2].valueAsText
            newThreshold = parameters[3].valueAsText
            runIter_bool = parameters[4].valueAsText
            fName_raster = parameters[5].valueAsText
            fName_polygon = parameters[6].valueAsText
            
        else:
            
            in_raster = parameters[0]
            wRadius = parameters[1]
            ischecked = parameters[2]
            newThreshold = parameters[3]
            runIter_bool = parameters[4]
            fName_raster = parameters[5]
            fName_polygon = parameters[6]
        
        #if box is not checked then proceed through normal routine.
        if (ischecked == "false"):
                
            spatial_ref = arcpy.Describe(in_raster).spatialReference
                #### FIRST STEP #######
                #Check projection
            if spatial_ref.linearUnitName == 'Meter':
                
                
                #evaluate first iteration as normal
                projected_raster = Raster(in_raster)
                print('Projected correctly, skipping reproject')
                arcpy.AddMessage("Projected correctly, skipping reproject")
                
                arcpy.AddMessage('Now calculating relief in moving window')
                print("Calculated Relief in Moving Window")
                neighborhood = NbrCircle(wRadius, "MAP")
                outFocalStats = FocalStatistics(in_raster, neighborhood, 'RANGE',"")
    
                
                ### THIRD STEP #### 
                 #Calculate thresholds for relief
                # rather than have the user input the threshold set it at something quantitative
                # like 1 or 2 s.d. below mean relief value
                arcpy.AddMessage('Setting threshold on relief to be 1 standard deviation')
                print('Setting threshold on relief to be 1 standard deviation')
                mean = arcpy.GetRasterProperties_management(outFocalStats,'MEAN')
                stdev = arcpy.GetRasterProperties_management(outFocalStats,'STD')
                n = 1                                                       #<-------- Fraction of Stdev could be increased/decreased
                #currently set to 1 so threshold is mean - 1 standard deviation. 
                threshold  = float(mean[0])-(float(stdev[0])/n)
                
                # extract portions of the raster that at 1 s.d. less rough than the mean
                threshHeld = (outFocalStats) <= (threshold)
                
                #Set all areas below threshold to null values
                threshHeld = SetNull(threshHeld, threshHeld, "VALUE = 0")
                
                if (runIter_bool == "true"): #if the box to runIter is checked then do one more step
                    outFocalStats2 = FocalStatistics(outFocalStats, neighborhood, 'RANGE',"")
        
                    
                    ### THIRD STEP #### 
                     #Calculate thresholds for relief
                    # rather than have the user input the threshold set it at something quantitative
                    # like 1 or 2 s.d. below mean relief value
                    arcpy.AddMessage('Setting threshold on relief to be 1 standard deviation')
                    print('Running a second iteration')
                    print('Setting threshold on relief to be 1 standard deviation')
                    mean = arcpy.GetRasterProperties_management(outFocalStats2,'MEAN')
                    stdev = arcpy.GetRasterProperties_management(outFocalStats2,'STD')
                    n = 1                                                       #<-------- Fraction of Stdev could be increased/decreased
                    #currently set to 1 so threshold is mean - 1 standard deviation. 
                    threshold  = float(mean[0])-(float(stdev[0])/n)
                    
                    # extract portions of the raster that at 1 s.d. less rough than the mean
                    threshHeld = (outFocalStats2) <= (threshold)
                    
                    #Set all areas below threshold to null values
                    threshHeld = SetNull(threshHeld, threshHeld, "VALUE = 0")
                    
                    print('Saving Output')
                    arcpy.AddMessage("Saving Output")
                    arcpy.RasterToPolygon_conversion(threshHeld, fName_polygon,'SIMPLIFY')
                    outFocalStats.save(fName_raster)
                   
                    
                else: #else save output as normal
                    print('Saving Output')
                    arcpy.AddMessage("Saving Output")
                    arcpy.RasterToPolygon_conversion(threshHeld, fName_polygon,'SIMPLIFY')
                    outFocalStats.save(fName_raster)
                    
    
                
                
            else:
                print('NOT projected with Meters as linear unit, please reproject')
                arcpy.AddError("NOT projected with Meters as linear unit, please reproject")
                # projected_raster = arcpy.ProjectRaster_management( raster_path,
                                                # out_raster,
                                                # proj_file_path,
                                                # 'BILINEAR')
        #This loop will execute if the user provides an alternate threshold 
        elif str(ischecked == 'true'):
            arcpy.AddMessage('Lower Threshold Given branching code')
            spatial_ref = arcpy.Describe(in_raster).spatialReference
                #### FIRST STEP #######
                #Check projection
            if spatial_ref.linearUnitName == 'Meter':
                projected_raster = Raster(in_raster)
                print('Projected correctly, skipping reproject')
                arcpy.AddMessage("Projected correctly, skipping reproject")
                
                arcpy.AddMessage('Now calculating relief in moving window')
                print("Calculated Relief in Moving Window")
                neighborhood = NbrCircle(wRadius, "MAP")
                outFocalStats = FocalStatistics(in_raster, neighborhood, 'RANGE',"")
    
                
                ### THIRD STEP #### 
                 #Calculate thresholds for relief
                # rather than have the user input the threshold set it at something quantitative
                # like 1 or 2 s.d. below mean relief value
                arcpy.AddMessage('Setting threshold on relief to be 1 standard deviation')
                print('Setting threshold on relief to be 1 standard deviation')
                mean = arcpy.GetRasterProperties_management(outFocalStats,'MEAN')
                stdev = arcpy.GetRasterProperties_management(outFocalStats,'STD')
                n = newThreshold                                                 #<-------- Fraction of Stdev could be increased/decreased
                #currently set to 1 so threshold is mean - 1 standard deviation. 
                threshold  = float(mean[0])-(float(stdev[0])*(float(n)))
                
                # extract portions of the raster that at 1 s.d. less rough than the mean
                threshHeld = (outFocalStats) <= (threshold)
                
                #Set all areas below threshold to null values
                threshHeld = SetNull(threshHeld, threshHeld, "VALUE = 0")
                
                
                if (runIter_bool == "true"): #if the box to runIter is checked then do one more step
                    outFocalStats2 = FocalStatistics(outFocalStats, neighborhood, 'RANGE',"")
        
                    
                    ### THIRD STEP #### 
                     #Calculate thresholds for relief
                    # rather than have the user input the threshold set it at something quantitative
                    # like 1 or 2 s.d. below mean relief value
                    arcpy.AddMessage('Setting threshold on relief to be 1 standard deviation')
                    print('Running a second iteration')
                    print('Setting threshold on relief to be 1 standard deviation')
                    mean = arcpy.GetRasterProperties_management(outFocalStats2,'MEAN')
                    stdev = arcpy.GetRasterProperties_management(outFocalStats2,'STD')
                    n = 1                                                       #<-------- Fraction of Stdev could be increased/decreased
                    #currently set to 1 so threshold is mean - 1 standard deviation. 
                    threshold  = float(mean[0])-(float(stdev[0])/n)
                    
                    # extract portions of the raster that at 1 s.d. less rough than the mean
                    threshHeld = (outFocalStats2) <= (threshold)
                    
                    #Set all areas below threshold to null values
                    threshHeld = SetNull(threshHeld, threshHeld, "VALUE = 0")
                    
                    print('Saving Output')
                    arcpy.AddMessage("Saving Output")
                    arcpy.RasterToPolygon_conversion(threshHeld, fName_polygon,'SIMPLIFY')
                    outFocalStats.save(fName_raster)
                   
                    
                else: #else save output as normal
                    print('Saving Output')
                    arcpy.AddMessage("Saving Output")
                    arcpy.RasterToPolygon_conversion(threshHeld, fName_polygon,'SIMPLIFY')
                    outFocalStats.save(fName_raster)
    
                
                
            else:
                print('NOT projected with Meters as linear unit, please reproject')
                arcpy.AddError("NOT projected with Meters as linear unit, please reproject")
                # projected_raster = arcpy.ProjectRaster_management( raster_path,
                                                # out_raster,
                                                # proj_file_path,
                                                # 'BILINEAR')

    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.
        """
        self.myfunc1(self, parameters)
        return
      

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\Low_Relief_Script\Test_area\USGS_12_n31_m.tif'
    # in_features = r'C:\GIS_Project\Low_Relief_Script\Test_area\USGS_13_n31w098.tif'
      wRadius = 1000
      in_filled_bool = 'true'
      cutoffThreshold = 1.0
      runIter_bool = "true"
      fName_raster = r'D:\GIS_Project\Low_Relief_Script\Test_area\test5.tif'
      fName_polygon = r'D:\GIS_Project\Low_Relief_Script\Test_area\test5.shp' 
     
      parameters = [in_features, wRadius,in_filled_bool, cutoffThreshold, runIter_bool, fName_raster, fName_polygon]
     
      tool1 = Low_relief_finder()
      tool1.callfrom_pyt = False
      tool1.execute(parameters, None)
     

    

        

