# -*- coding: utf-8 -*-
"""
/***************************************************************************
 BathyMorph
                                 A QGIS plugin
 Calculates Slope Contours and fills small gaps in bathymetry 
 Generated by Plugin Builder: http://g-sherman.github.io/Qgis-Plugin-Builder/
                              -------------------
        begin                : 2023-03-10
        git sha              : $Format:%H$
        copyright            : (C) 2023 by Dr 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 
from qgis.core import QgsProject 

# Import the code for the dialog
from .BathyMorph_dialog import BathyMorphDialog
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 QgsSingleBandGrayRenderer,QgsContrastEnhancement

import processing
import sys
import traceback
import os, glob

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

    def select_input_file1(self): 
        filename, _filter = QFileDialog.getOpenFileName(selfMT.dlg, "Select input bathymetry file ","", '*.img *.tif*') 
        # Add layer to frame, find last in list, add to end of list, create all new lists
        selfMT.dlg.lineCombo1.clear() 
        selfMT.dlg.lineCombo1.insertItem(0,filename)
        selfMT.dlg.lineCombo1.setCurrentIndex(0)
        #autofill
        autoPoly = filename[:-4]+"_interp.img"
        selfMT.dlg.lineEdit_2.setText(autoPoly)
        if os.path.exists(autoPoly):
            selfMT.dlg.exists1.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists1.setText("")
        autoPoly = filename[:-4]+"_slope.img"
        selfMT.dlg.lineEdit_3.setText(autoPoly)
        if os.path.exists(autoPoly):
            selfMT.dlg.exists2.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists2.setText("")
        autoPoly = filename[:-4]+"_contours.shp"
        selfMT.dlg.lineEdit_4.setText(autoPoly)
        if os.path.exists(autoPoly):
            selfMT.dlg.exists5.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists5.setText("")
        autoPoly = filename[:-4]+"_roughness.img"
        selfMT.dlg.lineEdit_5.setText(autoPoly)
        if os.path.exists(autoPoly):
            selfMT.dlg.exists3.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists3.setText("")
        autoPoly = filename[:-4]+"_hires.img"
        selfMT.dlg.lineEdit_6.setText(autoPoly)
        if os.path.exists(autoPoly):
            selfMT.dlg.exists4.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists4.setText("")
        autoPoly = filename[:-4]+"_BPI.img"
        selfMT.dlg.lineEdit_7.setText(autoPoly)
        if os.path.exists(autoPoly):
            selfMT.dlg.exists7.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists7.setText("")
    def select_input_file2(self): 
        filename, _filter = QFileDialog.getSaveFileName(selfMT.dlg, "Select output interpolated file ","", '*.img') 
        selfMT.dlg.lineEdit_2.setText(filename) 
        if os.path.exists(filename):
            selfMT.dlg.exists1.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists1.setText("")
    def select_input_file3(self): 
        filename, _filter = QFileDialog.getSaveFileName(selfMT.dlg, "Select output slope file","", '*.img') 
        selfMT.dlg.lineEdit_3.setText(filename) 
        if os.path.exists(filename):
            selfMT.dlg.exists2.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists2.setText("")
    def select_input_file4(self): 
        filename, _filter = QFileDialog.getSaveFileName(selfMT.dlg, "Select output contours shapefile","", '*.shp') 
        selfMT.dlg.lineEdit_4.setText(filename) 
        if os.path.exists(filename):
            selfMT.dlg.exists5.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists5.setText("")
    def select_input_file5(self): 
        filename, _filter = QFileDialog.getSaveFileName(selfMT.dlg, "Select output roughness file","", '*.img') 
        selfMT.dlg.lineEdit_5.setText(filename) 
        if os.path.exists(filename):
            selfMT.dlg.exists3.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists3.setText("")         
    def select_input_file6(self): 
        filename, _filter = QFileDialog.getSaveFileName(selfMT.dlg, "Select output high resolution file","", '*.img') 
        selfMT.dlg.lineEdit_6.setText(filename) 
        if os.path.exists(filename):
            selfMT.dlg.exists4.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists4.setText("")
    def select_input_file7(self): 
        filename, _filter = QFileDialog.getSaveFileName(selfMT.dlg, "Select output BPI file","", '*.img') 
        selfMT.dlg.lineEdit_7.setText(filename) 
        if os.path.exists(filename):
            selfMT.dlg.exists7.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists7.setText("")
    def indexChanged(self): 
        selectedLayerIndex = selfMT.dlg.lineCombo1.currentIndex()
        currentText = selfMT.dlg.lineCombo1.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.lineEdit_2.text()[0:len(currentText[:-4])]
        if filename1[0:3] == "_in" or currentText[:-4] == filename1[0:len(currentText[:-4])]:
            filename = currentText
        autoPoly = filename[:-4]+"_interp.img"
        selfMT.dlg.lineEdit_2.setText(autoPoly)
        if os.path.exists(autoPoly):
            selfMT.dlg.exists1.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists1.setText("")
        autoPoly = filename[:-4]+"_slope.img"
        selfMT.dlg.lineEdit_3.setText(autoPoly)
        if os.path.exists(autoPoly):
            selfMT.dlg.exists2.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists2.setText("")
        autoPoly = filename[:-4]+"_contours.shp"
        selfMT.dlg.lineEdit_4.setText(autoPoly)
        if os.path.exists(autoPoly):
            selfMT.dlg.exists5.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists5.setText("")
        autoPoly = filename[:-4]+"_roughness.img"
        selfMT.dlg.lineEdit_5.setText(autoPoly)
        if os.path.exists(autoPoly):
            selfMT.dlg.exists3.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists3.setText("")
        autoPoly = filename[:-4]+"_hires.img"
        selfMT.dlg.lineEdit_6.setText(autoPoly)
        if os.path.exists(autoPoly):
            selfMT.dlg.exists4.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists4.setText("")
        autoPoly = filename[:-4]+"_BPI.img"
        selfMT.dlg.lineEdit_7.setText(autoPoly)
        if os.path.exists(autoPoly):
            selfMT.dlg.exists7.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists7.setText("")

    def help(self): 
        import webbrowser
        import marinetools
        MThelp = os.path.dirname(marinetools.__file__) + "\\bathymorph\\BathyMorph.pdf"
        webbrowser.open(MThelp)  

    def run(self):
        import random
        import marinetools
        from marinetools.bathymorph.BathyMorph_dialog import BathyMorphDialog
        global selfMT
        """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
        if self.first_start == True:
            self.first_start = False
            self.dlg = BathyMorphDialog()
            selfMT = self
            self.dlg.pushButton_1.clicked.connect(BathyMorph.select_input_file1) 
            self.dlg.pushButton_2.clicked.connect(BathyMorph.select_input_file2) 
            self.dlg.pushButton_3.clicked.connect(BathyMorph.select_input_file3) 
            self.dlg.pushButton_4.clicked.connect(BathyMorph.select_input_file4) 
            self.dlg.pushButton_5.clicked.connect(BathyMorph.select_input_file5) 
            self.dlg.pushButton_6.clicked.connect(BathyMorph.select_input_file6) 
            self.dlg.pushButton_7.clicked.connect(BathyMorph.select_input_file7) 
            self.dlg.lineCombo1.currentIndexChanged.connect(BathyMorph.indexChanged)
            self.dlg.helpButton.clicked.connect(BathyMorph.help) 
            
        # Fetch the currently loaded layers
        layers = QgsProject.instance().layerTreeRoot().children() 
        layers = QgsProject.instance().mapLayers().values()
        # Clear the contents of the comboBox from previous runs
        # self.dlg = MBES_SegmentationDialog()
        self.dlg.lineCombo1.clear() 
        # Populate the comboBox with names of all the loaded layer   
        self.dlg.lineCombo1.addItems([layer.name() for layer in layers if str(layer.type())== "1" or str(layer.type())== "LayerType.Raster"])
        BathyMorph.indexChanged(self) 
        
        # show the dialog
        self.dlg.show()
        # Run the dialog event loop
        result = self.dlg.exec_()
        
        # See if OK was pressed

        if result:
            import glob
            selectedLayerIndex = self.dlg.lineCombo1.currentIndex()
            currentText = selfMT.dlg.lineCombo1.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
            filename1 = filename
            
            # Test input file to find if supported 
            from marinetools.fileTest import fileTest
            filename1 = fileTest.main(filename1,"raster")
            if filename1[0:5] == "Error":
                print(filename1)
                return
            
            filename2 = self.dlg.lineEdit_2.text()  
            filename3 = self.dlg.lineEdit_3.text()  
            filename4 = self.dlg.lineEdit_4.text()  
            filename5 = self.dlg.lineEdit_5.text()  
            filename6 = self.dlg.lineEdit_6.text()
            filename7 = self.dlg.lineEdit_7.text()

            Smooth = self.dlg.SmoothingFactor.text()
            Interval = self.dlg.ContourInterval.text()
            Resolution = self.dlg.ResolutionOut.text()
            inner = self.dlg.inner.text()
            outer = self.dlg.outer.text()
            if Smooth == "":
                Smooth = "7"
            if Interval == "":
                Interval = "10"
            if Resolution == "":
                Resolution = "100"
            if inner == "":
                inner = "2"
            if outer == "":
                outer = "10"
            if self.dlg.DoInterpolate.isChecked() == True:
                DoInterp=1
            else:
                DoInterp=0
            if self.dlg.DoSlope.isChecked() == True:
                DoSlope=1
            else:
                DoSlope=0
            if self.dlg.DoContours.isChecked() == True:
                DoContours=1
            else:
                DoContours=0
            if self.dlg.DoRoughness.isChecked() == True:
                DoRoughness=1
            else:
                DoRoughness=0
            if self.dlg.DoHiRes.isChecked() == True:
                DoHiRes=1
            else:
                DoHiRes=0
            if self.dlg.DoBPI.isChecked() == True:
                DoBPI=1
            else:
                DoBPI=0
            if self.dlg.DoStandardise.isChecked() == True:
                DoStandardise=1
            else:
                DoStandardise=0

            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
   
            newdir = str(os.path.dirname(filename1) + "/tempMT")
            if not os.path.exists(newdir):
                os.mkdir(newdir)

            alphabet = 'ZYXWVUTSRQPONMLKJIHGFEDCBA0987654321'
            rand = alphabet[random.randint(1,25)] + alphabet[random.randint(1,35)] + alphabet[random.randint(1,35)] + alphabet[random.randint(1,35)]
            #Print file report
            txtFile = open(str(filename1[:-4]) + "_BathyMorph_Info.txt", "w")
            txtFile.write("Script: Bathynetry Morphometry" "\n")
            txtFile.write("\n")
            txtFile.write("Input DEM: " + os.path.basename(filename1) + "\n")                
            txtFile.write("\n")
            
            if DoInterp:
                # Interpolate output to filename2
                if os.path.exists(filename2):
                    os.remove(filename2)
                newInterp = filename2
                
                #take the file and make it usable by raster calculator
                tempName = os.path.split(filename1)[1]
                if tempName.count('.') > 1:
                    name = tempName[:-4]
                else:
                    name = tempName.split('.')[0]
                layerRef = name + '@1'
                # change to 0.0 and 1.0
                expr = '"'+ layerRef + '" * 0.0 + 1.0'
                infile = QgsRasterLayer(filename1)
                crs = infile.crs()
                extent = infile.extent()
                print("Expression = "+str(expr))

                result = processing.run("native:rastercalc", {'LAYERS':filename1,'EXPRESSION':expr,'EXTENT':extent,
                                                     'CELL_SIZE':None,'CRS':crs,'OUTPUT':'TEMPORARY_OUTPUT'})
                polyCov = result['OUTPUT']
                result = processing.run("gdal:polygonize", {'INPUT':polyCov,'BAND':1,'FIELD':'DN',
                                                            'EIGHT_CONNECTEDNESS':False,'EXTRA':'',
                                                            'OUTPUT':'TEMPORARY_OUTPUT'})
                delHole = result['OUTPUT']
                result = processing.run("native:deleteholes", {'INPUT':delHole,'MIN_AREA':0,'OUTPUT':'TEMPORARY_OUTPUT'})
                coverage1 = result['OUTPUT']
                tempfile1 = newdir + "/tempfile1.tif"
                tempfile1a = newdir + "/tempfile1.tif.aux.xml"

                processing.run("gdal:fillnodata", {'INPUT':filename1,'BAND':1,'DISTANCE':Smooth,'ITERATIONS':0,'NO_MASK':False,
                                                   'MASK_LAYER':None,'OPTIONS':'','EXTRA':'','OUTPUT':tempfile1})
                #fillnodata1 = result['OUTPUT']
                processing.run("gdal:cliprasterbymasklayer", {'INPUT':tempfile1,'MASK':coverage1,'SOURCE_CRS':None,'TARGET_CRS':None,'TARGET_EXTENT':None,'NODATA':None,'ALPHA_BAND':False,'CROP_TO_CUTLINE':True,
                                                              'KEEP_RESOLUTION':False,'SET_RESOLUTION':False,'X_RESOLUTION':None,'Y_RESOLUTION':None,
                                                              'MULTITHREADING':False,'OPTIONS':None,'DATA_TYPE':0,'EXTRA':'','OUTPUT':newInterp})
                txtFile.write("Interpolated file     : " + str(newInterp) + "\n")
                txtFile.write("Interpolation distance: " + str(Smooth) + "\n")
                txtFile.write("\n")
                #if os.path.exists(tempfile1):
                #    os.remove(tempfile1)
                #if os.path.exists(tempfile1a):
                #   os.remove(tempfile1a)
            else:
                newInterp = filename1

            if DoSlope:
                # Calculate Slope output to filename3 but use newInterp
                if os.path.exists(filename3):
                    os.remove(filename3)
                processing.run("native:slope", {'INPUT':newInterp,'Z_FACTOR':1,'OUTPUT':filename3})
                txtFile.write("Slope file created   : " + str(filename3) + "\n")
                txtFile.write("\n")
                fname = os.path.dirname(str(filename3))
                vlayer = QgsRasterLayer(str(filename3), str(filename3[len(fname)+1:-4]))
                QgsProject.instance().addMapLayer(vlayer)
                
                renderer=vlayer.renderer()
                provider=vlayer.dataProvider()
                GrayRenderer = QgsSingleBandGrayRenderer(provider,1) #create a new renderer
                myType = renderer.dataType(1)
                myEnhancement = QgsContrastEnhancement(myType)
                contrast_enhancement = QgsContrastEnhancement.StretchToMinimumMaximum
                myEnhancement.setContrastEnhancementAlgorithm(contrast_enhancement,True)
                myEnhancement.setMinimumValue(0)   #Set the minimum 0
                vlayerStats = vlayer.dataProvider().bandStatistics(1,QgsRasterBandStats.All)
                vlayerMax = vlayerStats.maximumValue
                myEnhancement.setMaximumValue(vlayerMax)
                vlayer.setRenderer(GrayRenderer)
                vlayer.renderer().setContrastEnhancement(myEnhancement)
                vlayer.renderer().setGradient(1) # WhiteToBlack
                vlayer.triggerRepaint()

            
            if DoInterp:
                fname = os.path.dirname(str(newInterp))
                vlayer = QgsRasterLayer(str(newInterp), str(filename2[len(fname)+1:-4]))
                QgsProject.instance().addMapLayer(vlayer)

            if DoContours:
                # make contours output to filename4 but use newInterp
                if os.path.exists(filename4):
                    os.remove(filename4)
                dict = processing.run("native:rasterlayerproperties", {'INPUT':newInterp,'BAND':1})
                cellsizeX = float(dict['PIXEL_WIDTH'])
                cutoff = cellsizeX * 10.0
                result = processing.run("gdal:contour", {'INPUT':newInterp,'BAND':1,'INTERVAL':Interval,'FIELD_NAME':'ELEV','CREATE_3D':False,'IGNORE_NODATA':False,'NODATA':None,'OFFSET':0,'EXTRA':'','OUTPUT':'TEMPORARY_OUTPUT'})
                newContour1 = result['OUTPUT']
                result = processing.run("qgis:exportaddgeometrycolumns", {'INPUT':newContour1,'CALC_METHOD':0,'OUTPUT':'TEMPORARY_OUTPUT'})
                newContour2 = result['OUTPUT']
                processing.run("native:extractbyattribute", {'INPUT':newContour2,'FIELD':'length','OPERATOR':2,'VALUE':cutoff,'OUTPUT':filename4})
                txtFile.write("Contour file created: " + str(filename4) + "\n")
                txtFile.write("Contour interval    : " + str(Interval) + "m \n")
                txtFile.write("\n")
                fname = os.path.dirname(str(filename4))
                vlayer = QgsVectorLayer(str(filename4), str(filename4[len(fname)+1:-4]), "ogr")
                QgsProject.instance().addMapLayer(vlayer)

            if DoRoughness:
                # Calculate Roughness output to filename5 but use newInterp
                if os.path.exists(filename5):
                    os.remove(filename5)
                newRoughness = newdir + "/newRoughness"+rand+".img"
                processing.run("gdal:roughness", {'INPUT':newInterp,'BAND':1,'COMPUTE_EDGES':False,'OPTIONS':'','OUTPUT':filename5})
                txtFile.write("Roughness file created: " + str(filename5) + "\n")
                txtFile.write("\n")
                fname = os.path.dirname(str(filename5))
                vlayer = QgsRasterLayer(str(filename5), str(filename5[len(fname)+1:-4]))
                QgsProject.instance().addMapLayer(vlayer)

            if DoHiRes:
                # Make interpolated output to filename6 but use newInterp
                if os.path.exists(filename6):
                    os.remove(filename6)
                newHiRes = newdir + "/newHiRes"+rand+".gpkg"
                outTin = newdir + "/outTin"+rand+".img"
                processing.run("native:pixelstopoints", {'INPUT_RASTER':newInterp,'RASTER_BAND':1,'FIELD_NAME':'VALUE','OUTPUT':newHiRes})
                #result = processing.run("native:pixelstopoints", {'INPUT_RASTER':newInterp,'RASTER_BAND':1,'FIELD_NAME':'VALUE','OUTPUT':'TEMPORARY_OUTPUT'})
                #newHiRes = result['OUTPUT']
                expr = newHiRes + '|layername=newHiRes'+rand+'::~::0::~::1::~::0'
                extent = QgsVectorLayer(newdir + "/newHiRes"+rand+".gpkg").extent()
                crs = QgsVectorLayer(newdir + "/newHiRes"+rand+".gpkg").crs()
                processing.run("qgis:tininterpolation", {'INTERPOLATION_DATA':expr,'METHOD':0,'EXTENT':extent,'PIXEL_SIZE':Resolution,'OUTPUT':outTin})
                tempName = os.path.split(newInterp)[1]
                if tempName.count('.') > 1:
                    name = tempName[:-4] # remove file extension as filename has extra "." - hopefully just 3 letters
                else:
                    name = tempName.split('.')[0]
                InRef = name + '@1'
                tempName = os.path.split(outTin)[1]
                name = tempName.split('.')[0]
                TinRef = name + '@1'
                expression = '"' + TinRef + '" * "' + InRef + '" / "' + InRef + '"'
                processing.run("qgis:rastercalculator", {'EXPRESSION':expression,'LAYERS':[outTin],'CELLSIZE':None,'EXTENT':None,'CRS':crs,'OUTPUT':filename6})
                txtFile.write("High Res file created: " + str(filename6) + "\n")
                txtFile.write("New resolution       : " + str(Resolution) + "m \n")
                txtFile.write("\n")
                fname = os.path.dirname(str(filename6))
                vlayer = QgsRasterLayer(str(filename6), str(filename6[len(fname)+1:-4]))
                QgsProject.instance().addMapLayer(vlayer)
                
            if DoBPI:
                if os.path.exists(filename7):
                    os.remove(filename7)
                temp1 = newdir + r"\\outer.img"
                try:
                    remove_img_file(temp1)
                except:
                    a=1
                temp2 = newdir + r"\\inner.img"
                try:
                    remove_img_file(temp2)
                except:
                    a=1
                bathy_nonull = newdir + r"\\bathy_nonull.img"
                try:
                    remove_img_file(bathy_nonull)
                except:
                    a=1
                bathy_noStan = newdir + r"\\bathy_noStan.img"
                try:
                    remove_img_file(bathy_noStan)
                except:
                    a=1
                processing.run("native:fillnodata", {'INPUT':newInterp,'BAND':1,'FILL_VALUE':0,'OUTPUT':bathy_nonull})
                processing.run("otb:Smoothing", {'in':bathy_nonull,'out':temp1,'type':'mean','type.mean.radius':outer,'outputpixeltype':5})
                processing.run("otb:Smoothing", {'in':bathy_nonull,'out':temp2,'type':'mean','type.mean.radius':inner,'outputpixeltype':5})
                rado = float(outer) * float(outer) # area of smoothed area (less the pi)
                radi = float(inner) * float(inner)
                denom = rado-radi
                formula = "C - (((A*"+str(rado)+")-B*("+str(radi)+"))/"+str(denom)+")"
                if DoStandardise:
                    processing.run("gdal:rastercalculator", {'INPUT_A':temp1,'BAND_A':1,'INPUT_B':temp2,'BAND_B':1,'INPUT_C':newInterp,'BAND_C':1,'FORMULA':formula,'OUTPUT':bathy_noStan})
                    stats = processing.run("native:rasterlayerstatistics", {'INPUT':str(bathy_noStan),'BAND':1})
                    InputMean = stats["MEAN"]
                    InputStd = stats["STD_DEV"]
                    tempName = os.path.split(bathy_noStan)[1]
                    if tempName.count('.') > 1:
                        name = tempName[:-4]
                    else:
                        name = tempName.split('.')[0]
                    name = tempName
                    layerRef = name + '@1'
                    formula = '((("' + layerRef + '" - '+str(InputMean)+')/'+str(InputStd)+')*100.0)'
                    infile = QgsRasterLayer(newInterp)
                    crs = infile.crs()
                    extent = infile.extent()
                    processing.run("qgis:rastercalculator", {'EXPRESSION':formula,'LAYERS':[bathy_noStan],'CELLSIZE':None,'EXTENT':extent,'CRS':None,'OUTPUT':filename7})
                else:
                    processing.run("gdal:rastercalculator", {'INPUT_A':temp1,'BAND_A':1,'INPUT_B':temp2,'BAND_B':1,'INPUT_C':newInterp,'BAND_C':1,'FORMULA':formula,'OUTPUT':filename7})
                    
                txtFile.write("BPI file created: " + str(filename7) + "\n")
                txtFile.write("BPI outer radius: " + str(outer) + "\n") 
                txtFile.write("BPI inner radius: " + str(inner) + "\n")
                txtFile.write("\n")
                fname = os.path.dirname(str(filename7))
                vlayer = QgsRasterLayer(str(filename7), str(filename7[len(fname)+1:-4]))
                QgsProject.instance().addMapLayer(vlayer)
                
            self.loading_screen.stop_animation()
            txtFile.close()

            try:
                os.rmdir(newdir)
            except:
                print("Temporary directory is not empty to delete")
            pass

