# -*- coding: utf-8 -*-
"""
/***************************************************************************
 CoMMaDelin V1.2
                                 A QGIS plugin
 Delineate confined features on bathymetry data
 Generated by Plugin Builder: http://g-sherman.github.io/Qgis-Plugin-Builder/
                              -------------------
        begin                : 2025-01-07
        git sha              : $Format:%H$
        copyright            : (C) 2025 by Tim Le Bas, Riccardo Arosio, Joana Gafeira
        email                : tlb@noc.ac.uk
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
"""
from qgis.PyQt.QtCore import QSettings, QTranslator, QCoreApplication
from qgis.PyQt.QtGui import QIcon
from qgis.PyQt.QtWidgets import QAction,QFileDialog,QMessageBox 
from qgis.core import QgsProject 

#from .resources import *
# Import the code for the dialog
from .CoMMaDelin_dialog import CoMMaDelinDialog
import os.path

from qgis.core import Qgis,QgsMessageLog
from qgis.gui import QgsMessageBar

from qgis.core import QgsProcessing
from qgis.core import QgsProcessingAlgorithm
from qgis.core import QgsProcessingMultiStepFeedback
from qgis.core import QgsProcessingParameterRasterLayer
from qgis.core import QgsProcessingParameterRasterDestination
from qgis.core import QgsProcessingParameterDefinition
from qgis.core import QgsVectorLayer
from qgis.core import (QgsSymbol,QgsSimpleFillSymbolLayer,QgsRendererCategory,QgsCategorizedSymbolRenderer)
from qgis.core import QgsRasterLayer
from qgis.core import QgsRasterBandStats
from qgis.core import QgsLayerTreeLayer
from qgis.core import QgsProperty
from qgis.PyQt.QtWidgets import QMessageBox 

import processing

def remove_tif_files(tempfilename):
            import glob,os
            #templayer = QgsRasterLayer(tempfilename,"templayer")
            #QgsProject.instance().addMapLayer(templayer, False)
            #QgsProject.instance().removeMapLayer(templayer.id())
            for f in glob.glob(tempfilename.replace(".tif",".*")):
                os.remove(f)
            return
        
def remove_shp_files(tempfilename):  # not working yet
            import glob,os
            from osgeo import ogr
            try:
                templayer = QgsVectorLayer(tempfilename,"templayer")
                QgsProject.instance().addMapLayer(templayer, True)
                QgsProject.instance().takeMapLayer(templayer)
                QgsProject.instance().removeMapLayer(templayer)
                #for f in glob.glob(tempfilename.replace(".shp",".*")):
                #    print (str(f))
                #    os.remove(f)
                
                QgsProject.instance().reloadAllLayers()  # removes 2 of the 4 files
                driver = ogr.GetDriverByName("ESRI Shapefile")
                if os.path.exists(tempfilename):
                     driver.DeleteDataSource(tempfilename)
            except:
                print(tempfilename," not deleted")

            return
        
class LoadingScreenDlg:
    """Loading screen animation."""
    from qgis.PyQt.QtWidgets import QDialog, QLabel 
    from qgis.PyQt.QtGui import QMovie, QPalette, QColor

    def __init__(self, gif_path):
        self.dlg = self.QDialog()
        self.dlg.setWindowTitle("Please Wait")
        # self.dlg.setWindowModality(False)
        self.dlg.setFixedSize(200, 100)
        pal = self.QPalette()
        # role = self.QPalette.Background
        # pal.setColor(role, self.QColor(255, 255, 255))
        self.dlg.setPalette(pal)
        self.label_animation = self.QLabel(self.dlg)
        self.movie = self.QMovie(gif_path)
        self.label_animation.setMovie(self.movie)

    def start_animation(self):
        self.movie.start()
        self.dlg.show()
        return

    def stop_animation(self):
        self.movie.stop()
        self.dlg.done(0)       

class helpers(object):
    
    # This function converts backslash (accepted through the ArcGIS tool) to forwardslash (needed in python script) in a path
    def convert_backslash_forwardslash(self,inText):
        # inText: input path
        
        #inText = fr"{inText}"
        if inText.find('\t'):
            inText = inText.replace('\t', '\\t')
        elif inText.find('\n'):
            inText = inText.replace('\n', '\\n')
        elif inText.find('\r'):
            inText = inText.replace('\r', '\\r')

        inText = inText.replace('\\','/')
        
        return inText
    
class Arc2QGIS:
    def Negate(self,inputRaster):
        import os
        tempName = os.path.split(inputRaster)[1]
        if tempName.count('.') > 1:
            name = tempName[:-4]
        else:
            name = tempName.split('.')[0]
        layerRef = name + '@1'
        expr = '"'+ layerRef + '" * (-1.0)'
        infile = QgsRasterLayer(inputRaster)
        crs = infile.crs()
        extent = infile.extent()
        print(expr)
        result = processing.run("qgis:rastercalculator", {'EXPRESSION':expr,'LAYERS':[inputRaster],
                                                 'CELLSIZE':None,'EXTENT':extent,'CRS':crs,'OUTPUT':'TEMPORARY_OUTPUT'})
        invertedRaster = result['OUTPUT']
        return invertedRaster
    
    def Minus(self,inputRaster1,inputRaster2):
        import os
        tempName1 = os.path.split(inputRaster1)[1]
        if tempName1.count('.') > 1:
            name = tempName1[:-4]
        else:
            name = tempName1.split('.')[0]
        layerRef1 = name + '@1'
        tempName2 = os.path.split(inputRaster2)[1]
        if tempName2.count('.') > 1:
            name = tempName2[:-4]
        else:
            name = tempName2.split('.')[0]
        layerRef2 = name + '@1'
        expr = '"'+ layerRef1 + '" - "' + layerRef2 + '"'
        infile = QgsRasterLayer(inputRaster1)
        crs = infile.crs()
        extent = infile.extent()
        result = processing.run("qgis:rastercalculator", {'EXPRESSION':expr,'LAYERS':[inputRaster1,inputRaster2],
                                                 'CELLSIZE':None,'EXTENT':extent,'CRS':crs,'OUTPUT':'TEMPORARY_OUTPUT'})
        minusRaster = result['OUTPUT']
        return minusRaster
        

class CoMMaDelin:
    """QGIS Plugin Implementation."""


    def __init__(self, iface):

        # Save reference to the QGIS interface
        self.iface = iface
        # initialize plugin directory
        self.plugin_dir = os.path.dirname(__file__)
        # initialize locale
        locale = QSettings().value('locale/userLocale')[0:2]
        locale_path = os.path.join(
            self.plugin_dir,
            'i18n',
            'CoMMaDelin_{}.qm'.format(locale))

        if os.path.exists(locale_path):
            self.translator = QTranslator()
            self.translator.load(locale_path)
            QCoreApplication.installTranslator(self.translator)

        # Declare instance attributes
        self.actions = []
        self.menu = self.tr(u'&Confined morphologies delineators (CoMMa)')

        # Check if plugin was started the first time in current QGIS session
        # Must be set in initGui() to survive plugin reloads
        self.first_start = None

    # noinspection PyMethodMayBeStatic
    def tr(self, message):

        # noinspection PyTypeChecker,PyArgumentList,PyCallByClass
        return QCoreApplication.translate('CoMMaDelin', message)


    def add_action(
        self,
        icon_path,
        text,
        callback,
        enabled_flag=True,
        add_to_menu=True,
        add_to_toolbar=True,
        status_tip=None,
        whats_this=None,
        parent=None):

        icon = QIcon(icon_path)
        action = QAction(icon, text, parent)
        action.triggered.connect(callback)
        action.setEnabled(enabled_flag)

        if status_tip is not None:
            action.setStatusTip(status_tip)

        if whats_this is not None:
            action.setWhatsThis(whats_this)

        if add_to_toolbar:
            # Adds plugin icon to Plugins toolbar
            self.iface.addToolBarIcon(action)

        if add_to_menu:
            self.iface.addPluginToMenu(
                self.menu,
                action)

        self.actions.append(action)

        return action

    def initGui(self):
        """Create the menu entries and toolbar icons inside the QGIS GUI."""

        icon_path = ':/plugins/CoMMaDelin/icon.png'
        self.add_action(
            icon_path,
            text=self.tr(u'CoMMaDelin_v1.2'),
            callback=self.run,
            parent=self.iface.mainWindow())

        # will be set False in run()
        self.first_start = True


    def unload(self):
        """Removes the plugin menu item and icon from QGIS GUI."""
        for action in self.actions:
            self.iface.removePluginMenu(
                self.tr(u'&Confined morphologies delineators (CoMMa)'),
                action)
            self.iface.removeToolBarIcon(action)



    def select_input_file(self): 
        filename, _filter = QFileDialog.getOpenFileName(selfMT.dlg, "Select Input Raster DEM ","", '*.img *.tif') 
        # Add layer to frame, find last in list, add to end of list, create all new lists
        selfMT.dlg.inputDEM.clear() 
        selfMT.dlg.inputDEM.insertItem(0,filename)
        selfMT.dlg.inputDEM.setCurrentIndex(0)

    def select_inputLPTder_file(self): 
        filename, _filter = QFileDialog.getOpenFileName(selfMT.dlg, "Select LPT derivative raster ","", '*.img *.tif') 
        selfMT.dlg.inputDer.clear() 
        selfMT.dlg.inputDer.insertItem(0,filename)
        selfMT.dlg.inputDer.setCurrentIndex(0)


    def select_output_file(self): 
        filename, _filter = QFileDialog.getSaveFileName(selfMT.dlg, "Select Output Feature file ","", '*.shp') 
        selfMT.dlg.outFeat.setText(filename) 
        if os.path.exists(filename):
            selfMT.dlg.exists1.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists1.setText("")
         
    def select_inputGeo_file(self): 
        filename, _filter = QFileDialog.getOpenFileName(selfMT.dlg, "Select the geomorphons raster file ","", '*.img *.tif') 
        selfMT.dlg.inputGeo.clear() 
        selfMT.dlg.inputGeo.insertItem(0,filename)
        selfMT.dlg.inputGeo.setCurrentIndex(0)
        
    def indexChanged(self): 
        selectedLayerIndex = selfMT.dlg.inputDEM.currentIndex()
        currentText = selfMT.dlg.inputDEM.currentText()
        if selectedLayerIndex < 0:
            return
        layers = QgsProject.instance().mapLayers().values()
        a=0
        filename="NULL"
        for layer in (layer1 for layer1 in layers if str(layer1.type())== "1"):
            if a == selectedLayerIndex:
                filename = str(layer.source())
            a=a+1
        filename1= selfMT.dlg.outFeat.text()[0:len(currentText[:-4])]
        #if filename1[0:2] == "_D" or currentText[:-4] == filename1[0:len(currentText[:-4])]:
        if filename1[0:2] == "_D":
            filename = currentText
            
        infile = QgsRasterLayer(filename)
        result = processing.run("gdal:gdalinfo", {'INPUT':infile,'MIN_MAX':False,'STATS':False,'NOGCP':False,
                                              'NO_METADATA':False,'EXTRA':'','OUTPUT':'TEMPORARY_OUTPUT'})
        f= open(result['OUTPUT'])
        lines = f.read()
        f.close()
        projCrs = str(lines).find("PROJCRS")
        #if projCrs == -1 and currentText != "":
        #    QMessageBox.information(None, "Info:", "Check all 5 Parameters\nGeographic ")


        #dict = processing.run("native:rasterlayerproperties", {'INPUT':currentText,'BAND':1})
        #res = float(dict['PIXEL_WIDTH'])
        #check = selfMT.dlg.outFeat.text().find(currentText[:-4])
        dict = processing.run("native:rasterlayerproperties", {'INPUT':filename,'BAND':1})
        res = float(dict['PIXEL_WIDTH'])
        check = selfMT.dlg.outFeat.text().find(filename[:-4])
        #autofill    
        cutoffVR = selfMT.dlg.cutoffVR.text()
        if (cutoffVR == "" or check == -1):
            cutoffVR = str(res/5.0)
            selfMT.dlg.cutoffVR.setText('{0:.2f}'.format(float(cutoffVR)))   
        minVR = selfMT.dlg.minVR.text()
        if (minVR == "" or check == -1):
            minVR = str(res/2.0)
            selfMT.dlg.minVR.setText('{0:.2f}'.format(float(minVR)))   
        minWidth = selfMT.dlg.minWidth.text()
        if (minWidth == "" or check == -1):
            minWidth = str(res*5.0)
            selfMT.dlg.minWidth.setText('{0:.2f}'.format(float(minWidth)))   
        minRatio = selfMT.dlg.minRatio.text()
        if (minRatio == ""):
            minRatio = "0"
            selfMT.dlg.minRatio.setText("0")   
        delBuffer = selfMT.dlg.delBuffer.text()
        if (delBuffer == "" or check == -1):
            delBuffer = str(res)
            selfMT.dlg.delBuffer.setText('{0:.2f}'.format(float(delBuffer)))   
        newname = filename[:-4]+"_Delin_"+'{0:.2f}'.format(float(cutoffVR))+"_"+ \
                  '{0:.2f}'.format(float(minVR))+"_"+ \
                  '{0:.2f}'.format(float(minWidth))+"_"+ \
                  minRatio+"_"+'{0:.2f}'.format(float(minWidth))+".shp"
        selfMT.dlg.outFeat.setText(newname)

        if os.path.exists(newname):
            selfMT.dlg.exists1.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists1.setText("")
        return 
                
    def help(self): 
        import webbrowser
        import marinetools
        MThelp = os.path.dirname(marinetools.__file__) + "\\CoMMaDelin\\CoMMaDelin.pdf"
        webbrowser.open(MThelp)  
 
    def run(self):
        """Run method that performs all the real work"""

        # Create the dialog with elements (after translation) and keep reference
        # Only create GUI ONCE in callback, so that it will only load when the plugin is started
        import tempfile,glob,os
        import random
        import shutil
        import math
        global selfMT

        # Create the dialog with elements (after translation) and keep reference
        # Only create GUI ONCE in callback, so that it will only load when the plugin is started
        if self.first_start == True:
            self.first_start = False
            self.dlg = CoMMaDelinDialog()
            selfMT = self
            self.dlg.inputDEMpushButton.clicked.connect(CoMMaDelin.select_input_file) 
            self.dlg.inputDerpushButton.clicked.connect(CoMMaDelin.select_inputLPTder_file) 
            self.dlg.outFeatpushButton.clicked.connect(CoMMaDelin.select_output_file) 
            self.dlg.inputGeopushButton.clicked.connect(CoMMaDelin.select_inputGeo_file) 
            self.dlg.cutoffVR.textChanged.connect(CoMMaDelin.indexChanged) 
            self.dlg.minVR.textChanged.connect(CoMMaDelin.indexChanged) 
            self.dlg.minWidth.textChanged.connect(CoMMaDelin.indexChanged) 
            self.dlg.minRatio.textChanged.connect(CoMMaDelin.indexChanged) 
            self.dlg.delBuffer.textChanged.connect(CoMMaDelin.indexChanged) 
            self.dlg.helpButton.clicked.connect(CoMMaDelin.help)
            self.dlg.inputDEM.currentIndexChanged.connect(CoMMaDelin.indexChanged)
        
        layers = QgsProject.instance().mapLayers().values()
        self.dlg.inputDEM.clear() 
        self.dlg.inputDEM.addItems([layer.name() for layer in layers if str(layer.type())== "1"])
        self.dlg.inputDer.clear() 
        #self.dlg.inputDer.addItems([layer.name() for layer in layers])
        self.dlg.inputDer.addItems([layer.name() for layer in layers if str(layer.type())== "1"])
        self.dlg.inputDer.insertItem(0,"")
        max = selfMT.dlg.inputDer.setCurrentIndex(0)
        self.dlg.inputGeo.clear() 
        #self.dlg.inputGeo.addItems([layer.name() for layer in layers])
        self.dlg.inputGeo.addItems([layer.name() for layer in layers if str(layer.type())== "1"])
        self.dlg.inputGeo.insertItem(0,"")
        max = selfMT.dlg.inputGeo.setCurrentIndex(0)
        
        names = ["Positive","Negative"]
        self.dlg.fillDirec.clear() 
        self.dlg.fillDirec.addItems(names)
        
        names = ["No simplification or smoothing", "Light simplification and smoothing", "Heavy simplification and smoothing"]
        self.dlg.in_smoosimp.clear() 
        self.dlg.in_smoosimp.addItems(names)

        CoMMaDelin.indexChanged(self) 
        # show the dialog
        self.dlg.show()
        # Run the dialog event loop
        result = self.dlg.exec()
        
        # See if OK was pressed
        if result:
            selectedLayerIndex = self.dlg.inputDEM.currentIndex()
            currentText = selfMT.dlg.inputDEM.currentText()
            layers = QgsProject.instance().mapLayers().values()
            a=0
            filename="NULL"
            for layer in (layer1 for layer1 in layers if str(layer1.type())== "1"):
                if a == selectedLayerIndex:
                    filename = str(layer.source())
                a=a+1
            if currentText not in filename:
                filename = currentText
            inputDEM = filename 

            selectedLayerIndex = self.dlg.inputDer.currentIndex()
            currentText = selfMT.dlg.inputDer.currentText()
            layers = QgsProject.instance().mapLayers().values()
            a=1
            for layer in (layer1 for layer1 in layers if str(layer1.type())== "1"):
                if a == selectedLayerIndex:
                    filename = str(layer.source())
                a=a+1
            if currentText not in filename:
                filename = currentText
            inputDer = filename 
            if currentText == "":
                inputDer = ""
            
            fillDirec = selfMT.dlg.fillDirec.currentText()
            cutoffVR = selfMT.dlg.cutoffVR.text()
            if cutoffVR == "":
                cutoffVR = "2.5"
            minVR = selfMT.dlg.minVR.text()
            if minVR == "":
                minVR = "0.5"
            minWidth = selfMT.dlg.minWidth.text()
            if minWidth == "":
                minWidth = "5"
            minRatio = selfMT.dlg.minRatio.text()
            if minRatio == "":
                minRatio = "0"
            delBuffer = selfMT.dlg.delBuffer.text()
            if delBuffer == "":
                delBuffer = "1"
            in_smoosimp = selfMT.dlg.in_smoosimp.currentText()
            
            newdir = str(os.path.dirname(inputDEM) + "/tempMT")
            if not os.path.exists(newdir):
                os.mkdir(newdir)
            workspace = newdir
            outFeat = self.dlg.outFeat.text()
            if outFeat == "":
                self.iface.messageBar().pushMessage("\nNo outputfile was given")
                return
            if os.path.exists(outFeat):
                os.remove(outFeat)

            selectedLayerIndex = self.dlg.inputGeo.currentIndex()
            currentText = selfMT.dlg.inputGeo.currentText()
            layers = QgsProject.instance().mapLayers().values()
            a=1
            for layer in (layer1 for layer1 in layers if str(layer1.type())== "1"):
                if a == selectedLayerIndex:
                    filename = str(layer.source())
                a=a+1
            if currentText not in filename:
                filename = currentText
            inputGeo = filename 
            if currentText == "":
                inputGeo = ""

            minGeoRatio = selfMT.dlg.minGeoRatio.text()
            WBTyn = True
            if minGeoRatio == "":
                WBTyn = False
            delHoles = False
            if selfMT.dlg.delHoles.isChecked():
                delHoles = True
            delTemp = False
            if selfMT.dlg.delTemp.isChecked():
                delTemp = True
            print("inputDEM  ="+str(inputDEM))
            print("inputDer  ="+str(inputDer))
            print("fillDirec ="+str(fillDirec))
            print("cutoffVR  ="+str(cutoffVR))
            print("minVR     ="+str(minVR))
            print("minWidth  ="+str(minWidth))
            print("minRatio  ="+str(minRatio))
            print("delBuffer ="+str(delBuffer))
            print("in_smoosimp="+str(in_smoosimp))
            print("workspace  ="+str(workspace))
            print("outFeat   ="+str(outFeat))
            print("WBTyn     ="+str(WBTyn))
            print("inputGeo  ="+str(inputGeo))
            print("minGeoRatio  ="+str(minGeoRatio))
            print("delHoles  ="+str(delHoles))
            print("delTemp  ="+str(delTemp))
            
            plugin_dir = os.path.dirname(__file__)
            gif_path = os.path.join(plugin_dir, "loading.gif")
            self.loading_screen = LoadingScreenDlg(gif_path)  # init loading dlg
            self.loading_screen.start_animation()  # start loading dlg

            # Do something useful here - delete the line containing pass and
            # substitute with your code.
            # enable the helper functions
            helper = helpers()
            ArcGIS2QGIS = Arc2QGIS()
            inputDEM = helper.convert_backslash_forwardslash(inputDEM)
            inputDer = helper.convert_backslash_forwardslash(inputDer)
            inputGeo = helper.convert_backslash_forwardslash(inputGeo)
            workspace = helper.convert_backslash_forwardslash(workspace)
                        
            #Cellsize_Xresult = arcpy.GetRasterProperties_management(inputDEM, "CELLSIZEX")
            #Cellsize_X = float(Cellsize_Xresult.getOutput(0))
            dict = processing.run("native:rasterlayerproperties", {'INPUT':inputDEM,'BAND':1})
            Cellsize_X = float(dict['PIXEL_WIDTH'])

            tempWS = workspace
                              
            #list of fields that will be saved
            keepFields = []

            # Process: Invert DEM and fill
            print("Filling bathymetry ...")
            # Invert the input DEM to map positive topographic features
            if fillDirec == "Positive":           
                # Execute Negate
                featuresDEM = ArcGIS2QGIS.Negate(inputDEM)
                #featuresDEM = invertedInputDEM
                # Process: Fill
                filledDEM = tempWS + "/a1f.tif"
                result = processing.run("grass:r.fill.dir", {'input':featuresDEM,'format':0,'-f':False,
                                                             'output':filledDEM,'direction':'TEMPORARY_OUTPUT','areas':'TEMPORARY_OUTPUT'})
                #result = processing.run("sagang:fillsinksplanchondarboux2001", {'DEM':featuresDEM,'RESULT':filledDEM,'MINSLOPE':0.00})
                
            else:
                # Process: Fill
                featuresDEM = inputDEM
                filledDEM = tempWS + "/a1f.tif"
                result = processing.run("grass:r.fill.dir", {'input':featuresDEM,'format':0,'-f':False,
                                                             'output':filledDEM,'direction':'TEMPORARY_OUTPUT','areas':'TEMPORARY_OUTPUT'})
                #result = processing.run("sagang:fillsinksplanchondarboux2001", {'DEM':featuresDEM,'RESULT':filledDEM,'MINSLOPE':0.00})

            # Process: Minus 
            print("Subtracting ...")
            featuresDEMout2 = tempWS + "/featuresDEM.tif"
            #result = processing.run("sagang:copyraster", {'GRID':featuresDEM,'COPY':featuresDEMout2})
            shutil.copy(featuresDEM,featuresDEMout2)

            minusDEM = ArcGIS2QGIS.Minus(filledDEM, featuresDEMout2)
            
            print("Delineating features ...")
            #Delineating features
            #RECLASSIFY the raster using the Minimum Vertical Threshold.
            #Areas with values < cutoffVR will be classed as 1 and areas with values > Cutoff will be classed as 2
            reclassifyRanges = "" + cutoffVR + " 9000 1"
            reclassifyRanges = [str(-9999.0),str(cutoffVR),"1",str(cutoffVR),str(9999.0),"2"]
            print(reclassifyRanges)
            if inputDer != "":
                print("Preparing derivative ...")
                if fillDirec == "Negative":  
                    # Execute Negate
                    featuresDer = ArcGIS2QGIS.Negate(inputDer)     
                else:
                    featuresDer = inputDer
                result = processing.run("native:reclassifybytable", {'INPUT_RASTER':featuresDer,
                                                                     'RASTER_BAND':1,
                                                                     'TABLE':reclassifyRanges,
                                                                     'NO_DATA':0,
                                                                     'RANGE_BOUNDARIES':3,
                                                                     'NODATA_FOR_MISSING':True,
                                                                     'DATA_TYPE':0,'CREATE_OPTIONS':None,
                                                                     'OUTPUT':'TEMPORARY_OUTPUT'})
                OutReclassed = result['OUTPUT']
            else:
                result = processing.run("native:reclassifybytable", {'INPUT_RASTER':minusDEM,
                                                                     'RASTER_BAND':1,
                                                                     'TABLE':reclassifyRanges,
                                                                     'NO_DATA':0,
                                                                     'RANGE_BOUNDARIES':3,
                                                                     'NODATA_FOR_MISSING':True,
                                                                     'DATA_TYPE':0,'CREATE_OPTIONS':None,
                                                                     'OUTPUT':'TEMPORARY_OUTPUT'})
                OutReclassed = result['OUTPUT']
            
            
            #Convert to feature class
            #arcpy.RasterToPolygon_conversion(OutReclassed, tempWS + "b3p.shp", "NO_SIMPLIFY")
            Ras2Poly = tempWS + "/b3p.shp"
            result = processing.run("native:pixelstopolygons", {'INPUT_RASTER':OutReclassed,'RASTER_BAND':1,
                                                       'FIELD_NAME':'VALUE','OUTPUT':Ras2Poly})
            #result = processing.run("native:multiparttosingleparts", {'INPUT':Ras2Poly,'OUTPUT':'TEMPORARY_OUTPUT'})
            #Ras2PolySingle = result['OUTPUT']
            #Dissolved = tempWS + "/b3d.shp"
            #print(Dissolved)
            result = processing.run("gdal:dissolve", {'INPUT':Ras2Poly,'FIELD':'VALUE',
                                                      'GEOMETRY':'geometry','EXPLODE_COLLECTIONS':False,
                                                      'KEEP_ATTRIBUTES':False,'COUNT_FEATURES':False,
                                                      'COMPUTE_AREA':False,'COMPUTE_STATISTICS':False,
                                                      'STATISTICS_ATTRIBUTE':'','OPTIONS':'',
                                                      'OUTPUT':'TEMPORARY_OUTPUT'})
            Dissolved = result['OUTPUT']
            extracted = tempWS + "/b3p.shp"
            result = processing.run("native:extractbyattribute", {'INPUT':Dissolved,'FIELD':'VALUE','OPERATOR':2,
                                                                  'VALUE':'1','OUTPUT':extracted})
 
            result = processing.run("qgis:basicstatisticsforfields", {'INPUT_LAYER':extracted,'FIELD_NAME':'VALUE',
                                                                        'OUTPUT_HTML_FILE':'TEMPORARY_OUTPUT'})
            rows = int(result['FILLED'])

            if rows == 0:
                
                print("No features detected, please change the parameters ...")
            
            else:
                
                if float(delBuffer) > 0:
                    print("Buffering the delineated features ...")
                    result = processing.run("native:buffer", {'INPUT':extracted,'DISTANCE':delBuffer,
                                                              'SEGMENTS':5,'END_CAP_STYLE':0,'JOIN_STYLE':0,
                                                              'MITER_LIMIT':2,'DISSOLVE':False,'SEPARATE_DISJOINT':False,
                                                              'OUTPUT':'TEMPORARY_OUTPUT'})
                    Buffered = result['OUTPUT']
                    #arcpy.analysis.Buffer(tempWS + "b3p.shp", tempWS + "b3pbu.shp", delBuffer + " Meters", 
                    #                      "FULL", "ROUND", "NONE")
                    Dissolved2 = tempWS + "/b3pf.shp"

                    processing.run("native:multiparttosingleparts", {'INPUT':Buffered,'OUTPUT':Dissolved2})
                    #arcpy.management.Dissolve(tempWS + "b3pbu.shp", tempWS + "b3pf.shp", "", "", 
                    #                          "SINGLE_PART")
                    
                else:
                    Dissolved2 = tempWS + "/b3pf.shp"
                    #arcpy.management.Rename(tempWS + "b3p.shp", tempWS + "b3pf.shp")
                    os.rename(extracted,Dissolved2)
                result = processing.run("native:fieldcalculator", {'INPUT':Dissolved2,'FIELD_NAME':'row','FIELD_TYPE':1,
                                                          'FIELD_LENGTH':0,'FIELD_PRECISION':0,'FORMULA':' $id ',
                                                          'OUTPUT':'TEMPORARY_OUTPUT'})
                DissolvedWithRow = result['OUTPUT']
                
                #Characterise the features' geometry _ Part 2
                #Calculate the minimum bounding geometry (MBG) for each polygon and adding the 
                #following field: MBG_Width, MBG_Length, MBG_Orientation
                
                result = processing.run("native:orientedminimumboundingbox", {'INPUT':DissolvedWithRow,'OUTPUT':'TEMPORARY_OUTPUT'})
                Boxed = result['OUTPUT']
                #arcpy.MinimumBoundingGeometry_management(tempWS + "b3pf.shp", tempWS + "c4b.shp",
                #                                         "RECTANGLE_BY_WIDTH", "NONE", "", "MBG_FIELDS")
        
                #Add and calculate the MBG Width/Length field
                #arcpy.AddField_management(tempWS + "c4b.shp", "MBG_W_L", "DOUBLE", "8", "2")
                # 
                #arcpy.CalculateField_management(tempWS + "c4b.shp", "MBG_W_L",
                #                                '!MBG_Width!/!MBG_Length!', "PYTHON3")
                RatioCalc = tempWS + "/c4b.shp"
                result = processing.run("native:fieldcalculator", {'INPUT':Boxed,'FIELD_NAME':'MBG_W_L','FIELD_TYPE':0,
                                                          'FIELD_LENGTH':8,'FIELD_PRECISION':2,
                                                          'FORMULA':' "width" / "height" ','OUTPUT':RatioCalc})
                result = processing.run("native:fieldcalculator", {'INPUT':RatioCalc,'FIELD_NAME':'row2','FIELD_TYPE':1,
                                                          'FIELD_LENGTH':0,'FIELD_PRECISION':0,'FORMULA':' $id ',
                                                          'OUTPUT':'TEMPORARY_OUTPUT'})
                RatioCalcWithRow = result['OUTPUT']
                
                #Join the MBG fields to the feature outline polygon shapefile
                #arcpy.SpatialJoin_analysis(tempWS + "b3pf.shp", tempWS + "c4b.shp", tempWS + "d5g.shp",
                #                           "JOIN_ONE_TO_MANY","KEEP_ALL", "","WITHIN")
                Joined = tempWS + "/d5g.shp"
                processing.run("native:joinattributesbylocation", {'INPUT':DissolvedWithRow,'PREDICATE':[0,1,5],'JOIN':RatioCalcWithRow,
                                                                   'JOIN_FIELDS':[],'METHOD':0,'DISCARD_NONMATCHING':False,
                                                                   'PREFIX':'','OUTPUT':Joined})
                
                JoinedClean = tempWS + "/d5g1.shp"
                processing.run("native:extractbyexpression", {'INPUT':Joined,'EXPRESSION':' "row" = "row_2" ',
                                                              'OUTPUT':JoinedClean})
                expression = ' "width"  >= '+minWidth+ ' and "MBG_W_L" >= '+ minRatio
                #minsRemoved = tempWS + "/e6d.shp"
                result = processing.run("native:extractbyexpression", {'INPUT':JoinedClean,'EXPRESSION':expression,
                                                                       'OUTPUT':'TEMPORARY_OUTPUT'})
                minsRemoved = result['OUTPUT']
                
                #Reshaping polygon containing holes
                if delHoles:
                    print("Deleting holes ...")
                    result = processing.run("qgis:exportaddgeometrycolumns", {'INPUT':minsRemoved,'CALC_METHOD':0,'OUTPUT':'TEMPORARY_OUTPUT'})
                    exportadd = result['OUTPUT']
                    processing.run("native:fieldcalculator", {'INPUT':exportadd,'FIELD_NAME':'holesize','FIELD_TYPE':0,'FIELD_LENGTH':0,'FIELD_PRECISION':0,'FORMULA':'("area") * 0.55 ','OUTPUT':'TEMPORARY_OUTPUT'})
                    fieldCalc = result['OUTPUT']
                    #processing.run("native:deleteholes", {'INPUT':fieldCalc,'MIN_AREA':float('"holesize"'),'OUTPUT':EliminatePolygonPart})
                    result = processing.run("native:deleteholes", {'INPUT':fieldCalc,'MIN_AREA':QgsProperty.fromExpression('"holesize"'),
                                                                   'OUTPUT':'TEMPORARY_OUTPUT'})
                    EliminatePolygonPart = result['OUTPUT']

                    #arcpy.EliminatePolygonPart_management(tempWS + "e6d.shp", tempWS + "e6h.shp","PERCENT",
                    #                                     "","55","CONTAINED_ONLY")
                else:
                    #os.rename(minsRemoved,EliminatePolygonPart)
                    EliminatePolygonPart = minsRemoved
                
                
                #Getting Maximum Relief values from Mosaic raster
                #Use the created shapefile as a mask for extraction of the Maximum Relief 
                #value from the Mosaic raster
                print("Extracting Relief ...")
                
                #Add Maximum Relief values
                ZonalHeight = tempWS + "/f7vt.shp"
                result = processing.run("native:zonalstatisticsfb", {'INPUT':EliminatePolygonPart,'INPUT_RASTER':minusDEM,
                                                            'RASTER_BAND':1,'COLUMN_PREFIX':'Zonal_',
                                                            'STATISTICS':[6],'OUTPUT':ZonalHeight})
                #ZonalHeight = result['OUTPUT'] #f7vt
                
                expression = ' "Zonal_max" >= '+minVR
                lowVRRemoved = tempWS + "/f7se.shp"
                result = processing.run("native:extractbyexpression", {'INPUT':ZonalHeight,'EXPRESSION':expression,
                                                                       'OUTPUT':lowVRRemoved})
                result = processing.run("native:renametablefield", {'INPUT':lowVRRemoved,'FIELD':'Zonal_max',
                                                                    'NEW_NAME':'Relief','OUTPUT':'TEMPORARY_OUTPUT'})
                renamedField = result['OUTPUT']
                
                #Delete useless fields
                keepFields = ["Relief","width","height","MBG_W_L","angle"]
                result = processing.run("native:retainfields", {'INPUT':renamedField,'FIELDS':keepFields,
                                                       'OUTPUT':'TEMPORARY_OUTPUT'})
                retainedFields = result['OUTPUT']
                
                #Add area and perimeter fields
                result = processing.run("qgis:exportaddgeometrycolumns", {'INPUT':retainedFields,'CALC_METHOD':0,
                                                                          'OUTPUT':'TEMPORARY_OUTPUT'})
                addedArea = result['OUTPUT']

                ######################################################################
                # START GEOMORPHONS###################################################
                
                if WBTyn:
                    
                    print("Processing geomorphons ...")
                    
                    if fillDirec == "Positive":
                        #remapString  = "2 5 1; 1 0; 6 10 0"
                        remapString  = ['0','1','0',
                                        '1','5','1',
                                        '5','10','0']
                    else:
                        #remapString  = "1 6 0; 7 1; 8 0; 9 10 1"
                        remapString  = ['0','6','0',
                                        '6','7','1',
                                        '7','8','0',
                                        '8','10','1']
                    
                    #Reclassify(inputGeo, "Value", remapString, "NODATA")
                    result = processing.run("native:reclassifybytable", {'INPUT_RASTER':inputGeo,
                                                                               'RASTER_BAND':1,
                                                                               'TABLE':remapString,'NO_DATA':0,
                                                                               'RANGE_BOUNDARIES':0,
                                                                               'NODATA_FOR_MISSING':True,
                                                                               'DATA_TYPE':0,
                                                                               'CREATE_OPTIONS':None,
                                                                               'OUTPUT':'TEMPORARY_OUTPUT'})
                    OutRecla_Rid = result['OUTPUT']
                    
                    # MajorityFilter(OutRecla_Rid, "EIGHT", "HALF")
                    result = processing.run("grass:r.neighbors", {'input':OutRecla_Rid,'selection':None,
                                                                  'method':2,'size':1,'gauss':None,'quantile':'',
                                                                  '-c':False,'-a':False,'weight':'','output':'TEMPORARY_OUTPUT'})
                    outMajFilt_Rid = result['output']
                    
                    # ZonalStatisticsAsTable(tempWS + "f7se.shp", "FID",outMajFilt_Rid, tempWS + "f7ge","NODATA", "SUM")
                    result = processing.run("native:zonalstatisticsfb", {'INPUT':addedArea,'INPUT_RASTER':outMajFilt_Rid,
                                                                         'RASTER_BAND':1,'COLUMN_PREFIX':'_',
                                                                         'STATISTICS':[1],'OUTPUT':'TEMPORARY_OUTPUT'})
                    zonalStatsSum = result['OUTPUT']
                    
                    #arcpy.MakeFeatureLayer_management(tempWS + "f7se.shp", tempWS + "rcalc5_lyr")
                    # arcpy.AddJoin_management(tempWS + "rcalc5_lyr", "FID",tempWS + "f7ge", "FID", "KEEP_ALL")
                     
                    Pixel_area = Cellsize_X*Cellsize_X
                    relExpression =' "_sum" * ' + str(Pixel_area) + ' / "area" '
                    print(relExpression)
                    result = processing.run("native:fieldcalculator", {'INPUT':zonalStatsSum,'FIELD_NAME':'Geo_rat',
                                                                       'FIELD_TYPE':0,'FIELD_LENGTH':8,'FIELD_PRECISION':3,
                                                                       'FORMULA':relExpression,'OUTPUT':'TEMPORARY_OUTPUT'})
                    GeoRatAll = result['OUTPUT']
                    expression = ' "Geo_rat" >= '+minGeoRatio
                    result = processing.run("native:extractbyexpression", {'INPUT':GeoRatAll,'EXPRESSION':expression,
                                                                       'OUTPUT':'TEMPORARY_OUTPUT'})
                    poly2 = result['OUTPUT']
                    #Delete useless fields
                    keepFields = ["Relief","width","height","MBG_W_L","angle","area","perimeter","Geo_Rat"]
                    result = processing.run("native:retainfields", {'INPUT':poly2,'FIELDS':keepFields,
                                                           'OUTPUT':'TEMPORARY_OUTPUT'})
                    poly1 = result['OUTPUT']
                    #processing.run("native:deletecolumn", {'INPUT':poly2,'COLUMN':['_sum'],
                    #                                       'OUTPUT':'TEMPORARY_OUTPUT'})
                    #poly1 = result['OUTPUT']
                    
                else:
                    poly1 = addedArea
        
                ######################################################################
                ######################################################################
                # END GEOMORPHONS#####################################################                  

                #Reshaping the delineated polygons
                print("Reshaping the delineated polygons ...")
                result = processing.run("native:basicstatisticsforfields", {'INPUT_LAYER':poly1,'FIELD_NAME':'width',
                                                                   'OUTPUT':'TEMPORARY_OUTPUT',
                                                                   'OUTPUT_HTML_FILE':'TEMPORARY_OUTPUT'})
                min_width = result['MIN']
                
                if in_smoosimp == "No simplification or smoothing":
                    print("No simplification or smoothing")
                    smoothGeometry = poly1
                
                elif in_smoosimp == "Light simplification and smoothing":
                    
                    print("Light simplification and smoothing...")
            
                    # Simplifying
                    result = processing.run("native:coveragesimplify", {'INPUT':poly1,'TOLERANCE':float(min_width/10),
                                                                        'PRESERVE_BOUNDARY':False,
                                                                        'OUTPUT':'TEMPORARY_OUTPUT'})
                    simplify1 = result['OUTPUT']
                    result = processing.run("native:smoothgeometry", {'INPUT':simplify1,'ITERATIONS':1,'OFFSET':0.5,
                                                                      'MAX_ANGLE':180,'OUTPUT':'TEMPORARY_OUTPUT'})
                    smoothGeometry = result['OUTPUT']

                   
                elif in_smoosimp == "Heavy simplification and smoothing":
                    
                    print("Heavy simplification and smoothing...")
            
                    # Simplifying
                    result = processing.run("native:coveragesimplify", {'INPUT':poly1,'TOLERANCE':float(min_width/5),
                                                                        'PRESERVE_BOUNDARY':False,
                                                                        'OUTPUT':'TEMPORARY_OUTPUT'})
                    simplify1 = result['OUTPUT']
                    # Smoothing
                    result = processing.run("native:smoothgeometry", {'INPUT':simplify1,'ITERATIONS':1,'OFFSET':0.5,
                                                                      'MAX_ANGLE':180,'OUTPUT':'TEMPORARY_OUTPUT'})
                    smoothGeometry = result['OUTPUT']
               
                #Final operations
                print('Final cleaning ...')    
           
                result = processing.run("native:extractbyattribute", {'INPUT':smoothGeometry,'FIELD':'width',
                                                                      'OPERATOR':2,'VALUE':'0.0',
                                                                      'OUTPUT':outFeat})
                #Save file
                print('Saving official shapefile ...')
                fname = os.path.dirname(outFeat)
                vlayer = QgsVectorLayer(outFeat, str(outFeat[len(fname)+1:-4]), "ogr")
                QgsProject.instance().addMapLayer(vlayer)
                
                #Delete temporary files if requested/possible
                if delTemp:
                    print('Deleting temporary files ...')
                    remove_tif_files(filledDEM) #a1f
                    remove_shp_files(Ras2Poly) #b3p
                    remove_shp_files(Dissolved2) #b3pf
                    remove_shp_files(RatioCalc)  #c4b
                    remove_shp_files(Joined) #d5g
                    remove_shp_files(JoinedClean) #d5g1
                    remove_shp_files(ZonalHeight) #f7vt
                    remove_shp_files(lowVRRemoved) #f7se
                    remove_tif_files(featuresDEMout2) #featuresDEM
                #Print file report
                txtFile = open(str(outFeat[:-4]) + "_Info.txt", "w")
                txtFile.write("Script: CoMMa Delineation ToolBox v1.2" "\n")
                txtFile.write("\n")
                txtFile.write("File Name: " + str(outFeat) + "\n")
                txtFile.write("Input DEM: " + os.path.basename(inputDEM) + "\n")

                if inputDer != "":
                    txtFile.write("Delineation type: Fill algorithm on derivative" + "\n")
                    txtFile.write("Input derivative: " + os.path.basename(inputDer) + "\n")
                else:
                    txtFile.write("Delineation type: Fill algorithm on bathymetry" + "\n")
                    txtFile.write("Input derivative: none" + "\n")
                    
                txtFile.write("\n")
                txtFile.write("Relief type: " + fillDirec + "\n")
                txtFile.write("\n")
                txtFile.write("Input cell size: " + str(Cellsize_X) + "\n")
                txtFile.write("Cutoff value: " + cutoffVR + "\n")
                txtFile.write("Minimum relief value: " + minVR + "\n")
                txtFile.write("Minimum Width (m): " + minWidth + "\n")
                txtFile.write("Minimum W/L Ratio: " + minRatio + "\n")
                txtFile.write("Buffer extent (m): " + delBuffer + "\n")
                txtFile.write("Simplification or smoothing?: " + in_smoosimp + "\n")
                txtFile.write("\n")
                if WBTyn:
                    txtFile.write("Input geomorphons file: " + os.path.basename(inputGeo) + "\n")
                    txtFile.write("Geomorphons positives vs area ratio: " + minGeoRatio + "\n")
                else:
                    txtFile.write("The geomorphons add-on was not implemented")

                txtFile.close()
                self.loading_screen.stop_animation()



            return
