# -*- coding: utf-8 -*-
"""
/***************************************************************************
 BTM
                                 A QGIS plugin
  Analyzes benthic terrain for the purposes of classifying surficial seafloor characteristics that may be used in studies of benthic habitat, geomorphology, prediction of benthic fish species distribution, marine protected area design, and more
 Generated by Plugin Builder: http://g-sherman.github.io/Qgis-Plugin-Builder/
                              -------------------
        begin                : 2023-03-04
        git sha              : $Format:%H$
        copyright            : (C) 2023 by Tim Le Bas, National Oceanography Centre, Southampton. UK.
        email                : tim.lebas@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 # added
from qgis.core import QgsProject # added

# Import the code for the dialog
#from .BTM_dialog import BTMDialog
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

import processing
import sys
import traceback
import os
import csv
from xml.dom.minidom import parse
import shutil

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 BTM:
    """QGIS Plugin Implementation."""

    def select_input_file(self): 
        filename, _filter = QFileDialog.getOpenFileName(selfMT.dlg, "Select input raster Bathymetry file ","", '*.img *.tif') # added
        # Add layer to frame, find last in list, add to end of list, create all new lists
        selfMT.dlg.BathyInput.clear() 
        selfMT.dlg.BathyInput.insertItem(0,filename)
        selfMT.dlg.BathyInput.setCurrentIndex(0)
        #autofill
        BroadInner = selfMT.dlg.BroadInner.text()
        if BroadInner == "":
            BroadInner = "5"
        BroadOuter = selfMT.dlg.BroadOuter.text()
        if BroadOuter == "":
            BroadOuter = "50"
        FineInner = selfMT.dlg.FineInner.text()
        if FineInner == "":
            FineInner = "2"
        FineOuter = selfMT.dlg.FineOuter.text()
        if FineOuter == "":
            FineOuter = "10"
        autoPoly = filename[:-4]+"_BTM_"+BroadInner+"_"+BroadOuter+"_"+FineInner+"_"+FineOuter+".shp"
        autorast = filename[:-4]+"_zones_"+BroadInner+"_"+BroadOuter+"_"+FineInner+"_"+FineOuter+".img"
        selfMT.dlg.OutputPoly.setText(autoPoly)
        if os.path.exists(autoPoly):
            selfMT.dlg.exists1.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists1.setText("")
        selfMT.dlg.OutputRaster.setText(autorast)
        if os.path.exists(autorast):
            selfMT.dlg.exists2.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists2.setText("")
       
    def indexChanged(self): 
        selectedLayerIndex = selfMT.dlg.BathyInput.currentIndex()
        currentText = selfMT.dlg.BathyInput.currentText()
        layers = QgsProject.instance().mapLayers().values()
        a=0
        filename="NULL"
        for layer in (layer1 for layer1 in layers if str(layer1.type())== "1" or str(layer1.type())== "LayerType.Raster"):
            if a == selectedLayerIndex:
                filename = str(layer.source())
            a=a+1
        filename1= selfMT.dlg.OutputPoly.text()[0:len(currentText[:-4])]
        if filename1[0:3] == "_BT" or currentText[:-4] == filename1[0:len(currentText[:-4])]:
            filename = currentText
        #autofill
        BroadInner = selfMT.dlg.BroadInner.text()
        if BroadInner == "":
            BroadInner = "5"
        BroadOuter = selfMT.dlg.BroadOuter.text()
        if BroadOuter == "":
            BroadOuter = "50"
        FineInner = selfMT.dlg.FineInner.text()
        if FineInner == "":
            FineInner = "2"
        FineOuter = selfMT.dlg.FineOuter.text()
        if FineOuter == "":
            FineOuter = "10"
        autoPoly = filename[:-4]+"_BTM_"+BroadInner+"_"+BroadOuter+"_"+FineInner+"_"+FineOuter+".shp"
        autorast = filename[:-4]+"_zones_"+BroadInner+"_"+BroadOuter+"_"+FineInner+"_"+FineOuter+".img"
        selfMT.dlg.OutputPoly.setText(autoPoly)
        if os.path.exists(autoPoly):
            selfMT.dlg.exists1.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists1.setText("")
        selfMT.dlg.OutputRaster.setText(autorast)
        if os.path.exists(autorast):
            selfMT.dlg.exists2.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists2.setText("")

    def select_dictionary_file(self): # added
        filename, _filter = QFileDialog.getOpenFileName(selfMT.dlg, "Select classification dictionary file ","", '*.csv') # added
        selfMT.dlg.Dictionary.setText(filename) # added
        
    def select_outputRaster_file(self): # added
        filename, _filter = QFileDialog.getSaveFileName(selfMT.dlg, "Select output raster zones file ","", '*.img') # added
        selfMT.dlg.OutputRaster.setText(filename) # added
        if os.path.exists(filename):
            selfMT.dlg.exists2.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists2.setText("")
        
    def select_outputPoly_file(self): # added
        filename, _filter = QFileDialog.getSaveFileName(selfMT.dlg, "Select output polygon file ","", '*.shp') # added
        selfMT.dlg.OutputPoly.setText(filename) # added
        if os.path.exists(filename):
            selfMT.dlg.exists1.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists1.setText("")
        
    def help(self): 
        import webbrowser
        import marinetools
        MThelp = os.path.dirname(marinetools.__file__) + "\\btm\\BTM.pdf"
        webbrowser.open(MThelp)
        
    def otbInstall(self):
        import webbrowser
        import marinetools
        MThelp = os.path.dirname(marinetools.__file__) + "\\obia\\OTBinstall.pdf"
        webbrowser.open(MThelp)

    def run(self):
        """Run method that performs all the real work"""
        import tempfile,glob,os
        import random
        import shutil
        import math
        import marinetools
        from marinetools.btm.BTM_dialog import BTMDialog
        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
        self.first_start = True
        if self.first_start == True:
            self.first_start = False
            self.dlg = BTMDialog()
            selfMT = self
            self.dlg.InputFile.clicked.connect(BTM.select_input_file) # added
            self.dlg.DictFile.clicked.connect(BTM.select_dictionary_file) # added
            self.dlg.OutFileRaster.clicked.connect(BTM.select_outputRaster_file) # added
            self.dlg.OutFilePoly.clicked.connect(BTM.select_outputPoly_file) # added
            self.dlg.BroadInner.textChanged.connect(BTM.indexChanged) 
            self.dlg.BroadOuter.textChanged.connect(BTM.indexChanged) 
            self.dlg.FineInner.textChanged.connect(BTM.indexChanged) 
            self.dlg.FineOuter.textChanged.connect(BTM.indexChanged) 
            self.dlg.helpButton.clicked.connect(BTM.help) 
            self.dlg.BathyInput.currentIndexChanged.connect(BTM.indexChanged)
        # Fetch the currently loaded layers
        layers = QgsProject.instance().mapLayers().values()
        self.dlg.BathyInput.clear() 
        # Populate the comboBox with names of all the loaded layer   
        self.dlg.BathyInput.addItems([layer.name() for layer in layers if str(layer.type())== "1" or str(layer.type())== "LayerType.Raster"])
        BTM.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.BathyInput.currentIndex()
            currentText = selfMT.dlg.BathyInput.currentText()
            layers = QgsProject.instance().mapLayers().values()
            a=0
            filename="NULL"
            for layer in (layer1 for layer1 in layers if str(layer1.type())== "1" or str(layer1.type())== "LayerType.Raster"):
                if a == selectedLayerIndex:
                    filename = str(layer.source())
                a=a+1
            if currentText not in filename:
                filename = currentText
            bathy = filename

            # Test input file to find if supported 
            from marinetools.fileTest import fileTest
            bathy = fileTest.main(bathy,"raster")
            if bathy[0:5] == "Error":
                print(bathy)
                return
            
            out_dictionary        = self.dlg.Dictionary.text()  
            output_zones          = self.dlg.OutputRaster.text()  
            output_zones_Poly     = self.dlg.OutputPoly.text()  
            broad_bpi_inner_radius=self.dlg.BroadInner.text() 
            broad_bpi_outer_radius=self.dlg.BroadOuter.text() 
            fine_bpi_inner_radius =self.dlg.FineInner.text()  
            fine_bpi_outer_radius =self.dlg.FineOuter.text()
            try:
                processing.run("otb:ImageEnvelope", {'in':bathy,'out':'TEMPORARY_OUTPUT','sr':0,'elev.dem':'','elev.geoid':'','elev.default':0,'proj':''})
            except:
                QMessageBox.information(None, "Information:", "This tool requires (OrfeoToolBox (OTB) to be installed\nDownload from https://www.orfeo-toolbox.org/download/\n and install locally.")
                BTM.otbInstall(self)
                return
            if os.path.exists(output_zones):
                os.remove(output_zones)
            if os.path.exists(output_zones_Poly):
                os.remove(output_zones_Poly)
            if self.dlg.DeleteInter.isChecked() == True:
                DeleteInter=1
            else:
                DeleteInter=0
            if self.dlg.InvertDepths.isChecked() == True:
                InvertDepths=1
            else:
                InvertDepths=0

            if broad_bpi_inner_radius == "":
                broad_bpi_inner_radius = "5"
            if broad_bpi_outer_radius == "":
                broad_bpi_outer_radius = "50"
            if fine_bpi_inner_radius == "":
                fine_bpi_inner_radius = "2"
            if fine_bpi_outer_radius == "":
                fine_bpi_outer_radius = "10"
            # Do something useful here - delete the line containing pass and
            # substitute with your code.
            print(str(bathy))
            print(str(broad_bpi_inner_radius))
            print(str(broad_bpi_outer_radius))
            print(str(fine_bpi_inner_radius))
            print(str(fine_bpi_outer_radius))
            print(str(out_dictionary))
            print(str(output_zones))
            print(str(DeleteInter))
            print(str(InvertDepths))
            print(str(output_zones_Poly))
            #Print file report
            txtFile = open(str(output_zones_Poly[:-4]) + "_Info.txt", "w")
            txtFile.write("Script: Benthic Terrain Modeller" "\n")
            txtFile.write("\n")
            txtFile.write("File Name:     " + str(output_zones_Poly) + "\n")
            txtFile.write("Output Raster: " + str(output_zones) + "\n")
            txtFile.write("Input DEM:     " + os.path.basename(bathy) + "\n")                
            txtFile.write("\n")
            txtFile.write("Broad-scale BPI inner radius: " + str(broad_bpi_inner_radius) + "\n")
            txtFile.write("Broad-scale BPI outer radius: " + str(broad_bpi_outer_radius) + "\n")
            txtFile.write("Fine-scale BPI inner radius:  " + str(fine_bpi_inner_radius) + "\n")
            txtFile.write("Fine-scale BPI outer radius:  " + str(fine_bpi_outer_radius) + "\n")
            txtFile.write("\n")
            if DeleteInter == 1:
                txtFile.write("Intermediate files have been deleted where possible\n")
            else:
                txtFile.write("Intermediate files have not been deleted\n")

            if InvertDepths == 1:
                txtFile.write("Depths have been inverted\n")
            else:
                txtFile.write("Depths have not been inverted\n")

            
            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
            
            # if no dictionary provided make one's own
            basePath = str(os.path.dirname(bathy))
            if str(out_dictionary) == "":
                
                if not os.path.exists(str(basePath) + r"\tempMT"):
                    os.makedirs(str(basePath) + r"\tempMT")
                out_dictionary = str(basePath) + r"\tempMT\temporary_dictionary.csv"
                tempfile = str(basePath) + r"\tempMT\tmpoutputpoints.img"
                
                stats = processing.run("native:rasterlayerstatistics", {'INPUT':str(bathy),'BAND':1})
                InputMean = stats["MEAN"]
                InputStd = stats["STD_DEV"]
                if float(str(InputMean)) < 0.0:
                    InputMean = str(0.0 - float(str(InputMean)))
                processing.run("native:slope", {'INPUT':str(bathy),'Z_FACTOR':1,'OUTPUT':tempfile})
                stats = processing.run("native:rasterlayerstatistics", {'INPUT':str(tempfile),'BAND':1})
                SlopeMean = stats["MEAN"]
                SlopeStd = stats["STD_DEV"]

                self.iface.messageBar().pushMessage("Bathymetry Average = "+str(InputMean)+"\n")
                self.iface.messageBar().pushMessage("Bathymetry StanDev = "+str(InputStd)+"\n")
                self.iface.messageBar().pushMessage("     Slope Average = "+str(SlopeMean)+"\n")
                self.iface.messageBar().pushMessage("     Slope StanDev = "+str(SlopeStd)+"\n")

                file = open(out_dictionary,"w")

                file.write("Class,Zone,BroadBPI_Lower,BroadBPI_Upper,FineBPI_Lower,FineBPI_Upper,Slope_Lower,Slope_Upper,Depth_Lower,Depth_Upper\n") 
                file.write("1,Peak,-10000,-100,-10000,-100,0,90,0,12000\n") 
                file.write("2,Ridge,-10000,-100,-100,100,0,90,0,12000\n") 
                file.write("3,Ridge Trough,-10000,-100,100,10000,0,90,0,12000\n") 
                file.write("4,Flat with Ridge,-100,100,-10000,-100,0,90,0,12000\n") 
                file.write("5,Flat with Trough,-100,100,100,10000,0,90,0,12000\n") 
                file.write("6,Flat (Deeper),-100,100,-100,100,0,"+str(SlopeMean)+","+str(InputMean)+",12000\n") 
                file.write("7,Flat (Shallower),-100,100,-100,100,0,"+str(SlopeMean)+",0,"+str(InputMean)+"\n") 
                file.write("8,Steep (Deeper),-100,100,-100,100,"+str(SlopeMean)+",90,"+str(InputMean)+",12000\n") 
                file.write("9,Steep (Shallower),-100,100,-100,100,"+str(SlopeMean)+",90,0,"+str(InputMean)+"\n") 
                file.write("10,Depression with Ridge,100,10000,-10000,-100,0,90,0,12000\n") 
                file.write("11,Depression,100,10000,-100,100,0,90,0,12000\n") 
                file.write("12,Depression with Trough,100,10000,100,10000,0,90,0,12000\n") 
                file.close() 
                rememberDict = basePath + r"\Dictionary.csv"
                shutil.copyfile(out_dictionary,rememberDict)
                BTM.remove_img_file(tempfile)
            file = open(out_dictionary,"r") 
            lines = file.read()
            file.close()
            txtFile.write("Classification Dictionary:    " + str(out_dictionary) + "\n")
            txtFile.write(lines)
            txtFile.close()
                
            BTM.runMod(
                input_bathymetry      =str(bathy),
                broad_bpi_inner_radius=str(broad_bpi_inner_radius),
                broad_bpi_outer_radius=str(broad_bpi_outer_radius),
                fine_bpi_inner_radius =str(fine_bpi_inner_radius),
                fine_bpi_outer_radius =str(fine_bpi_outer_radius),
                classification_dict   =out_dictionary,
                output_zones          =str(output_zones),
                keepInter             =DeleteInter,
                invert                =InvertDepths)
            
            alphabet = 'ZYXWVUTSRQPONMLKJIHGFEDCBA0987654321'
            rand = alphabet[random.randint(1,25)] + alphabet[random.randint(1,35)] + alphabet[random.randint(1,35)] + alphabet[random.randint(1,35)]
            temp3 =basePath + "/tempMT/RAsieve.tif"
            temp4 =basePath + "/tempMT/VCpolygons"+rand+".shp"
            temp5 =basePath + "/tempMT/VGaggregate"+rand+".shp"
            temp6 =basePath + "/tempMT/VTaddfield"+rand+".shp"
            MinSize=10
            power = math.log(int(MinSize),2)+1.0 # Increments of power of two sizes
            first = str(output_zones)
            for x in range(1,int(power)):
                next=basePath + "/tempMT/RAsieve"+str(x)+".tif"
                Size = 2**x
                processing.run("gdal:sieve", {'INPUT':first,'THRESHOLD':Size,'EIGHT_CONNECTEDNESS':False,
                                              'NO_MASK':False,'MASK_LAYER':None,'EXTRA':'','OUTPUT':next})
                #self.remove_tif_files(first)
                first = next
            processing.run("gdal:sieve", {'INPUT':first,'THRESHOLD':MinSize,'EIGHT_CONNECTEDNESS':False,
                                          'NO_MASK':False,'MASK_LAYER':None,'EXTRA':'','OUTPUT':temp3})
            BTM.remove_tif_files(first)

            processing.run("native:pixelstopolygons", {'INPUT_RASTER':temp3,'RASTER_BAND':1,'FIELD_NAME':'VALUE','OUTPUT':temp4})
            self.iface.messageBar().pushMessage("Make polygons")
            processing.run("native:aggregate", {'INPUT':temp4,'GROUP_BY':'"VALUE"',
                                                'AGGREGATES':[{'aggregate': 'mean','delimiter': ',','input': '"VALUE"','length': 20,'name': 'VALUE','precision': 8,'sub_type': 0,'type': 6,'type_name': 'double precision'}],
                                                'OUTPUT':temp5})

            temp7 =basePath + "/tempMT/VTaddfield"+rand+".xlsx"
            temp8 =basePath + "/tempMT/VTaddfield2"+rand+".shp"
            temp9 =basePath + "/tempMT/VTaddfield3"+rand+".shp"
            processing.run("native:savefeatures", {'INPUT':out_dictionary,'OUTPUT':temp7,'LAYER_NAME':'','DATASOURCE_OPTIONS':'','LAYER_OPTIONS':''})
            processing.run("native:fieldcalculator", {'INPUT':temp5,'FIELD_NAME':'name2','FIELD_TYPE':2,'FIELD_LENGTH':0,'FIELD_PRECISION':0,'FORMULA':' to_string( VALUE )','OUTPUT':temp8})
            processing.run("native:joinattributestable", {'INPUT':temp8,'FIELD':'name2','INPUT_2':temp7,'FIELD_2':'Field1','FIELDS_TO_COPY':[],'METHOD':1,'DISCARD_NONMATCHING':False,'PREFIX':'Dict_','OUTPUT':temp6})
            processing.run("native:renametablefield", {'INPUT':temp6,'FIELD':'Dict_Fie_1','NEW_NAME':'Geomorphol','OUTPUT':temp9})
            processing.run("native:multiparttosingleparts", {'INPUT':temp9,'OUTPUT':output_zones_Poly})
        
            self.loading_screen.stop_animation()

            fname = os.path.dirname(str(output_zones_Poly))
            vlayer = QgsVectorLayer(str(output_zones_Poly), str(output_zones_Poly[len(fname)+1:-4]), "ogr")
            QgsProject.instance().addMapLayer(vlayer)
            BTM.colour_polygons_random(vlayer)
            if DeleteInter:
                import glob
                delClasses = basePath + "/tempMT/*.*"
                for f in glob.glob(delClasses):
                    try:
                        os.remove(f)
                    except:
                        print(f + " not removed")
                """delinnner = basePath + "/tempMT/inner*.*"
                for f in glob.glob(delinnner):
                    os.remove(f)
                delouter = basePath + "/tempMT/outer*.*"
                for f in glob.glob(delouter):
                    os.remove(f)
                """
        return

    def runMod(input_bathymetry, broad_bpi_inner_radius,
             broad_bpi_outer_radius, fine_bpi_inner_radius,
             fine_bpi_outer_radius, classification_dict, output_zones, keepInter, invert):

        """
        Compute complete model. The crux of this computation maps ranges
        of values provided in the classification dictionary (a CSV or Excel
        spreadsheet) to bathymetry derivatives: standardized
        fine- and broad- scale BPI and slope.
        """
        # Local workspace variables:
        basePath = str(os.path.dirname(input_bathymetry))
        out_workspace = basePath + r"\tempMT"        

        # local variables:
        broad_bpi = os.path.join(out_workspace, "broad_bpi")
        fine_bpi = os.path.join(out_workspace, "fine_bpi")
        slope_rast = os.path.join(out_workspace, "slope")
        broad_std = os.path.join(out_workspace, "broad_std")
        fine_std = os.path.join(out_workspace, "fine_std")

        MainRaster = basePath + r"\\tempMT\\MainRaster.img"
        try:
            remove_tif_file(MainRaster)
        except:
            a=1
        slope_rast = basePath + r"\\tempMT\\SlopeRast.img"
        try:
            remove_img_file(slope_rast)
        except:
            a=1
            
        if invert == 1:
            processing.run("gdal:rastercalculator", {'INPUT_A':input_bathymetry,'BAND_A':1,'FORMULA':'A*(-1.0)','OUTPUT':MainRaster})
        else:
            processing.run("gdal:rastercalculator", {'INPUT_A':input_bathymetry,'BAND_A':1,'FORMULA':'A*(1.0)','OUTPUT':MainRaster})

        # Process: Build Broad Scale BPI
        broad_bpi = basePath + r"\tempMT\\BroadBPI.img"
        try:
            remove_img_file(broad_bpi)
        except:
            a=1
        fine_bpi = basePath + r"\\tempMT\\FineBPI.img"
        try:
            remove_img_file(fine_bpi)
        except:
            a=1
        BTM.bpi3(MainRaster, broad_bpi_inner_radius,
                 broad_bpi_outer_radius, broad_bpi)

        # Process: Build Fine Scale BPI
        BTM.bpi3(MainRaster, fine_bpi_inner_radius,
                 fine_bpi_outer_radius, fine_bpi)
        
        broad_std = basePath + r"\\tempMT\\BroadStd.img"
        try:
            remove_img_file(broad_std)
        except:
            a=1
        fine_std = basePath + r"\\tempMT\\FineStd.img"
        try:
            remove_img_file(fine_std)
        except:
            a=1
 
        # Process: Standardize BPIs
        BTM.standardize_bpi(broad_bpi, broad_std)
        BTM.standardize_bpi(fine_bpi, fine_std)

        # Process: Calculate Slope
        processing.run("native:slope", {'INPUT':MainRaster,'Z_FACTOR':1,'OUTPUT':slope_rast})

        # Process: Zone Classification Builder
        BTM.classify(classification_dict, broad_std, fine_std, slope_rast, MainRaster, output_zones)

        return

    def bpi3(bathy=None, inner_radius=None, outer_radius=None,
             out_raster=None):
        """
        Create a bathymetric position index (BPI) raster, which
        measures the average value in a 'donut' of locations, excluding
        cells too close to the origin point, and outside a set distance.
        """
        basePath = str(os.path.dirname(bathy))
        temp1 = basePath + r"\\outertemp.img"
        try:
            remove_img_file(temp1)
        except:
            a=1
        temp2 = basePath + r"\\innertemp.img"
        try:
            remove_img_file(temp2)
        except:
            a=1
        bathy_nonull = basePath + r"\\bathy_nonull.img"
        try:
            remove_img_file(bathy_nonull)
        except:
            a=1
        processing.run("native:fillnodata", {'INPUT':bathy,'BAND':1,'FILL_VALUE':0,'OUTPUT':bathy_nonull})
        processing.run("otb:Smoothing", {'in':bathy_nonull,'out':temp1,'type':'mean','type.mean.radius':outer_radius,'outputpixeltype':5})
        processing.run("otb:Smoothing", {'in':bathy_nonull,'out':temp2,'type':'mean','type.mean.radius':inner_radius,'outputpixeltype':5})
        rado = float(outer_radius) * float(outer_radius) # area of smoothed area (less the pi)
        radi = float(inner_radius) * float(inner_radius)
        denom = rado-radi
        formula = "C - (((A*"+str(rado)+")-B*("+str(radi)+"))/"+str(denom)+")"
        #print(formula)
        processing.run("gdal:rastercalculator", {'INPUT_A':temp1,'BAND_A':1,'INPUT_B':temp2,'BAND_B':1,'INPUT_C':bathy,'BAND_C':1,'FORMULA':formula,'OUTPUT':out_raster})
        return

    def standardize_bpi(bpi_raster=None, out_raster=None):
        stats = processing.run("native:rasterlayerstatistics", {'INPUT':str(bpi_raster),'BAND':1})
        InputMean = stats["MEAN"]
        InputStd = stats["STD_DEV"]
        tempName = os.path.split(bpi_raster)[1]
        if tempName.count('.') > 1:
            name = tempName[:-4]
        else:
            name = tempName.split('.')[0]
        layerRef = name + '@1'

        formula = "(((A - "+str(InputMean)+")/"+str(InputStd)+")*100.0)"
        formula = '((("' + layerRef + '" - '+str(InputMean)+')/'+str(InputStd)+')*100.0)'
        print(formula)
        #processing.run("gdal:rastercalculator", {'INPUT_A':bpi_raster,'BAND_A':1,'FORMULA':formula,'OUTPUT':out_raster})
        infile = QgsRasterLayer(bpi_raster)
        crs = infile.crs()
        extent = infile.extent()
        processing.run("qgis:rastercalculator", {'EXPRESSION':formula,'LAYERS':[bpi_raster],'CELLSIZE':None,'EXTENT':extent,'CRS':None,'OUTPUT':out_raster})

        return
    
    def classify(classification_file, bpi_broad_std, bpi_fine_std, slope, bathy, out_raster):
        """
        Perform raster classification, based on classification mappings
        and provided raster derivatives (fine- and broad- scale BPI,
        slope, and the original raster). Outputs a classified raster.
        """

        basePath = os.path.dirname(bathy)

        # Read in the BTM Document; the class handles parsing a variety of inputs.
        btm_doc = BtmDocument(classification_file)
        classes = btm_doc.classification()
        tmpClass = basePath + r"\\tmpClass.img"
        try:
            remove_img_file(fine_std)
        except:
            a=1
        tmpClassIn = basePath + r"\\tmpClass1.img"
        processing.run("gdal:rastercalculator", {'INPUT_A':bathy,'BAND_A':1,'FORMULA':'A*0.0','OUTPUT':tmpClassIn})
        bro = QgsRasterLayer(bpi_broad_std,"bro")        
        fin = QgsRasterLayer(bpi_fine_std,"fin")        
        slo = QgsRasterLayer(slope,"slo")        
        dep = QgsRasterLayer(bathy,"dep")        
        QgsProject.instance().addMapLayer(bro, False)
        QgsProject.instance().addMapLayer(fin, False)
        QgsProject.instance().addMapLayer(slo, False)
        QgsProject.instance().addMapLayer(dep, False)
        key = {'0': 'None'}
        for item in classes:
            cur_class = str(item["Class"])
            cur_name = str(item["Zone"])
            tmpClassIn = basePath + r"\\tmpClass"+str(int(item["Class"]))+".img"
            tmpClassOut = basePath + r"\\tmpClass"+str(int(item["Class"])+1)+".img"
            cla = QgsRasterLayer(tmpClassIn,"cla")        
            QgsProject.instance().addMapLayer(cla, False)

            DL = str(float(str(item["Depth_LowerBounds"])))
            DU = str(float(str(item["Depth_UpperBounds"])))
            SL = str(float(str(item["Slope_LowerBounds"])))
            SU = str(float(str(item["Slope_UpperBounds"])))
            FL = str(float(str(item["LSB_LowerBounds"])))
            FU = str(float(str(item["LSB_UpperBounds"])))
            BL = str(float(str(item["SSB_LowerBounds"])))
            BU = str(float(str(item["SSB_UpperBounds"])))

            express = ' ( "bro@1" >'+BL+') * ( "bro@1" <'+BU+') * ( "fin@1" >'+FL+') * ( "fin@1" <'+FU+') * ( "slo@1" >'+SL+') * ( "slo@1" <'+SU+') * ( "dep@1" >'+DL+') * ( "dep@1" <'+DU+') * '+ str(item["Class"])+ ' + "cla@1"' 
            processing.run("qgis:rastercalculator", {'EXPRESSION':express,'LAYERS':bathy,'CELLSIZE':None,'EXTENT':None,'CRS':None,'OUTPUT':tmpClassOut})
            QgsProject.instance().removeMapLayer(cla.id())
        processing.run("gdal:rastercalculator", {'INPUT_A':tmpClassOut,'BAND_A':1,'FORMULA':'A*1.0', 'OUTPUT':out_raster})
        QgsProject.instance().removeMapLayer(bro.id())
        QgsProject.instance().removeMapLayer(fin.id())
        QgsProject.instance().removeMapLayer(slo.id())
        QgsProject.instance().removeMapLayer(dep.id())

        return
    
    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_img_file(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(".img",".*")):
            os.remove(f)
        return
         
    def colour_polygons_random(layer):
        # provide file name index and field's unique values
        from random import randrange

        fni = layer.fields().indexFromName('Geomorphol')
        unique_values = layer.uniqueValues(fni)

        # fill categories
        categories = []
        for unique_value in unique_values:
            # initialize the default symbol for this geometry type
            symbol = QgsSymbol.defaultSymbol(layer.geometryType())
            # configure a symbol layer
            layer_style = {}
            layer_style['color'] = '%d, %d, %d' % (randrange(0, 256), randrange(0, 256), randrange(0, 256))
            layer_style['outline'] = '#000000'
            symbol_layer = QgsSimpleFillSymbolLayer.create(layer_style)
            # replace default symbol layer with the configured one
            if symbol_layer is not None:
                symbol.changeSymbolLayer(0, symbol_layer)
            # create renderer object
            category = QgsRendererCategory(unique_value, symbol, str(unique_value))
            # entry for the list of category items
            categories.append(category)
        # create renderer object
        renderer = QgsCategorizedSymbolRenderer('Geomorphol', categories)

        # assign the created renderer to the layer
        if renderer is not None:
            layer.setRenderer(renderer)

        layer.triggerRepaint()
  
class BtmDocument(object):
    """ A wrapper class for handling any kind of BTM Classification file.
    The actual parsing happens in other classes, but this way we can keep
    a consistent method to access the data regardless of the input file type.
    """
    def __init__(self, filename):
        self.filename = filename
        self.doctype = self._doctype()
        self.schema = self._get_schema()

    def _doctype(self):
        # map of 'known' extensions to filetypes. don't bother with
        # a more formal mimetyping for now.
        known_types = {'.csv' : 'CSV','.xls' : 'Excel','.xlsx': 'Excel','.xml': 'XML'}
        ext = os.path.splitext(self.filename)[1].lower()
        #if known_types.has_key(ext):
        if (ext in known_types):
            dtype = known_types[ext]
        else:
            raise TypeError("Invalid document type for {}".format(self.filename))
        return dtype

    def _get_schema(self):
        """ map file types to their respective classes. """
        schema_map = {
            'CSV' : 'BtmCsvDocument',
            'Excel' : 'BtmExcelDocument',
            'XML' : 'BtmXmlDocument'
        }
        if (self.doctype in schema_map):
            schema_name = schema_map[self.doctype]
        # do some introspection to pull in the class, based on the document
        # type provided.
        schema = globals()[schema_name]
        return schema(self.filename)

    # push up methods contained within our subclasses
    def name(self):
        return self.schema.name()

    def description(self):
        return self.schema.description()

    def classification(self):
        return self.schema.classification()

class BtmXmlDocument(BtmDocument):
    def __init__(self, filename):
        self.dom = parse(filename)
        self.node_dict = self.node_to_dict(self.dom)

    def name(self):
        return self.node_dict['ClassDict']['PrjName']

    def description(self):
        return self.node_dict['ClassDict']['PRJDescription']

    def classification(self):
        # XXX: improve this to provide more direct access to the classes
        # currently returns a list of dictionaries, each one with a key:val pair;
        # would be useful to just have a header row, and map the names
        # from the header
        return self.node_dict['ClassDict']['Classifications']['ClassRec']

    def get_text_from_node(self, node):
        """
        scans through all children of node and gathers the
        text.
        """
        t = ""
        emptyNode = node.hasChildNodes()
        if emptyNode:
            for n in node.childNodes:
                if n.nodeType == n.TEXT_NODE:
                    t += n.nodeValue
                else:
                    raise NotTextNodeError
        else:
            t = None
        return t

    def node_to_dict(self, node):
        dic = {}
        multlist = {} # holds temporary lists where there are multiple children
        multiple = 0
        for n in node.childNodes:
            if n.nodeType != n.ELEMENT_NODE:
                continue
            # find out if there are multiple records
            if len(node.getElementsByTagName(n.nodeName)) > 1:
                multiple = 1
                # and set up the list to hold the values
                if not multlist.has_key(n.nodeName):
                    multlist[n.nodeName] = []
            try:
                # text node
                text = self.get_text_from_node(n)
            except NotTextNodeError:
                if multiple:
                    # append to our list
                    multlist[n.nodeName].append(self.node_to_dict(n))
                    dic.update({n.nodeName:multlist[n.nodeName]})
                    continue
                else:
                    # 'normal' node
                    dic.update({n.nodeName:self.node_to_dict(n)})
                    continue
            # text node
            if multiple:
                multlist[n.nodeName].append(text)
                dic.update({n.nodeName:multlist[n.nodeName]})
            else:
                dic.update({n.nodeName:text})
        return dic



class BtmExcelDocument(BtmDocument):
    # TODO: FORCE A TEST FOR ARCGIS 10.2 INSTALLATION

    def __init__(self, filename):
        self.filename = filename
        self.header = None # filled in by parse_workbook
        self.workbook = self.parse_workbook(self.filename)

    def name(self):
        return os.path.basename(self.filename)

    def description(self):
        # TODO: implement metadata in Excel?
        return None

    def classification(self):
        in_workbook = self.workbook
        result_rows = []
        for row in in_workbook:
            # replace empty strings with Nones
            row_clean = [None if x == '' else x for x in row]

            if len(row_clean) != 10:
                message = "Encountered malformed row which requires correction:" + \
                        "{}\"".format(",".join(row))
                raise ValueError(message)
            # don't parse the header, assume columns are in expected order.
            (class_code, zone, broad_lower, broad_upper, fine_lower, fine_upper, \
            slope_lower, slope_upper, depth_lower, depth_upper) = row_clean

            # for now: fake the format used by the XML documents.
            res_row = {'Class': str(int(class_code)),
                       'Zone': zone,
                       'SSB_LowerBounds': broad_lower,
                       'SSB_UpperBounds': broad_upper,
                       'LSB_LowerBounds': fine_lower,
                       'LSB_UpperBounds': fine_upper,
                       'Slope_LowerBounds': slope_lower,
                       'Slope_UpperBounds': slope_upper,
                       'Depth_LowerBounds': depth_lower,
                       'Depth_UpperBounds': depth_upper}
            result_rows.append(res_row)
        return result_rows

    def parse_workbook(self, filename):
        result = []
        try:
            import xlrd
        except ImportError:
            msg = "Reading Excel files requires the `xlrd` library, which is "+ \
                    "included in ArcGIS 10.2+. If you'd like Excel support in"+\
                    " ArcGIS 10.1, please install `xlrd` manually."
            raise Exception(msg)

        with xlrd.open_workbook(filename) as workbook:
            # assume data is in the first sheet
            sheet = workbook.sheet_by_index(0)

            # FIXME: assume all column labels are fixed
            self.header = ['Class', 'Zone', 'bpi_broad_lower', 'bpi_broad_upper',
                    'bpi_fine_lower', 'bpi_fine_upper', 'slope_lower',
                    'slope_upper', 'depth_lower', 'depth_upper']
            for row in range(2, sheet.nrows):
                cell = sheet.cell(row, 0)
                # an empty row terminates the set
                if cell.value in ["", None] :
                    break
                # we have an expected row of classes.
                else:
                    result.append([sheet.cell(row, col).value for col in range(10)])
        return result


class BtmCsvDocument(BtmDocument):
    def __init__(self, filename):
        self.filename = filename
        self.header = None # filled in by parse_csv
        self.csv = self.parse_csv(self.filename)

    def name(self):
        return os.path.basename(self.filename)

    def description(self):
        # TODO: implement metadata in CSV dictionaries?
        return None

    def classification(self):
        in_csv = self.csv
        result_rows = []
        for row in in_csv:
            # replace empty strings with Nones
            row_clean = [None if x == '' else x for x in row]

            # don't parse the header, assume columns are in expected order.
            (class_code, zone, broad_lower, broad_upper, fine_lower, fine_upper, \
            slope_lower, slope_upper, depth_lower, depth_upper) = row_clean

            # for now: fake the format used by the XML documents.
            res_row = {'Class': class_code,
                       'Zone': zone,
                       'SSB_LowerBounds': broad_lower,
                       'SSB_UpperBounds': broad_upper,
                       'LSB_LowerBounds': fine_lower,
                       'LSB_UpperBounds': fine_upper,
                       'Slope_LowerBounds': slope_lower,
                       'Slope_UpperBounds': slope_upper,
                       'Depth_LowerBounds': depth_lower,
                       'Depth_UpperBounds': depth_upper}
            result_rows.append(res_row)
        return result_rows

    def parse_csv(self, filename):
        result = None
        with open(filename, 'r') as f:
            # Use the sniffer to figure out what kind of input we're getting
            sample = f.read(1024)
            f.seek(0)
            sniff_obj = csv.Sniffer()
            try:
                dialect = sniff_obj.sniff(sample)
                has_header = sniff_obj.has_header(sample)
            except csv.Error:
                # If the CSV is malformed (e.g. a missing ',' in a row),
                # this error can be raised. In that case, we shouldn't
                # give up, but instead just set the default dialect and
                # assume a header.
                dialect = "excel"
                has_header = True

            # read in CSV, respecting the detected dialect
            in_csv = csv.reader(f, dialect)

            if has_header:
                self.header = next(in_csv)

            # everything but the header
            result = [r for r in in_csv]
            f.close()
        return result

