# -*- coding: utf-8 -*-
"""
/***************************************************************************
 Synthetic
                                 A QGIS plugin
 Removes Spurious pixels from bathymetry data that can look like trees
 Generated by Plugin Builder: http://g-sherman.github.io/Qgis-Plugin-Builder/
                              -------------------
        begin                : 2025-03-18
        git sha              : $Format:%H$
        copyright            : (C) 2025 by Dr Tim Le Bas
        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.                                   *
 *                                                                         *
 ***************************************************************************/
"""
import sys
import traceback
import os, glob
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 

# Import the code for the dialog
from .Synthetic_dialog import SyntheticDialog
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
from qgis.core import QgsProperty

import processing



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

    def __init__(self, iface):
        """Constructor.

        :param iface: An interface instance that will be passed to this class
            which provides the hook by which you can manipulate the QGIS
            application at run time.
        :type iface: QgsInterface
        """
        # 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',
            'Synthetic_{}.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'&Tree Removal')

        # 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):
        """Get the translation for a string using Qt translation API.

        We implement this ourselves since we do not inherit QObject.

        :param message: String for translation.
        :type message: str, QString

        :returns: Translated version of message.
        :rtype: QString
        """
        # noinspection PyTypeChecker,PyArgumentList,PyCallByClass
        return QCoreApplication.translate('Synthetic', 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/Synthetic/icon.png'
        self.add_action(
            icon_path,
            text=self.tr(u'Synthetic'),
            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'&Sythetic'),
                action)
            self.iface.removeToolBarIcon(action)
    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.BathyCombo1.clear() 
        selfMT.dlg.BathyCombo1.insertItem(0,filename)
        selfMT.dlg.BathyCombo1.setCurrentIndex(0)
        #autofill
        Range = selfMT.dlg.Range.text()
        if Range == "":
            Range = "3"
        Resolution = selfMT.dlg.Resolution.text()
        if Resolution == "":
            Resolution = "10"
        autoPoly = filename[:-4]+"_synthetic_"+Range+".img"
        selfMT.dlg.Output_lineEdit.setText(autoPoly)
        if os.path.exists(autoPoly):
            selfMT.dlg.exists2.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists2.setText("")
    def select_input_file2(self): 
        filename, _filter = QFileDialog.getSaveFileName(selfMT.dlg, "Select output interpolated file ","", '*.img') 
        selfMT.dlg.Output_lineEdit.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.getOpenFileName(selfMT.dlg, "Select input trackline file ","", '*.shp') 
        # Add layer to frame, find last in list, add to end of list, create all new lists
        selfMT.dlg.TracklineCombo1.clear() 
        selfMT.dlg.TracklineCombo1.insertItem(0,filename)
        selfMT.dlg.TracklineCombo1.setCurrentIndex(0)
        #autofill
        Range = selfMT.dlg.Range.text()
        if Range == "":
            Range = "100"
        autoPoly = filename[:-4]+"_synthetic_"+Range+".img"
        selfMT.dlg.Output_lineEdit.setText(autoPoly)
        if os.path.exists(autoPoly):
            selfMT.dlg.exists2.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists2.setText("")
        layer = QgsVectorLayer(filename, "ogr")
        fields = layer.fields()
        selfMT.dlg.comboBox.clear() 
        # Populate the comboBox with names of all the loaded fields   
        selfMT.dlg.comboBox.addItems([field.name() for field in fields]) 
    def VarDepChanged(self): 
        selectedLayerIndex = selfMT.dlg.TracklineCombo1.currentIndex()
        currentText = selfMT.dlg.TracklineCombo1.currentText()
        layers = QgsProject.instance().mapLayers().values()
        a=0
        filename="NULL"
        for layer in (layer1 for layer1 in layers if str(layer1.type())== "0" or str(layer1.type())== "LayerType.Vector"):
            if a == selectedLayerIndex:
                filename = str(layer.source())
            a=a+1
        if currentText not in filename:
            filename = currentText

        surveyLine = filename 
        layer = QgsVectorLayer(surveyLine, "ogr")
        fields = layer.fields()
        selfMT.dlg.comboBox.clear() 
        selfMT.dlg.comboBox.addItems([field.name() for field in fields]) 
        selfMT.dlg.ConstDep.clear() 
        selfMT.dlg.ConstAlt.clear() 
    def ConstDepChanged(self): 
        selfMT.dlg.comboBox.clear() 
        selfMT.dlg.ConstAlt.clear() 
    def ConstAltChanged(self): 
        selfMT.dlg.comboBox.clear() 
        selfMT.dlg.ConstDep.clear() 
    def indexChanged(self): 
        selectedLayerIndex = selfMT.dlg.BathyCombo1.currentIndex()
        currentText = selfMT.dlg.BathyCombo1.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.Output_lineEdit.text()[0:len(currentText[:-4])]
        if filename1[0:3] == "_sy" or currentText[:-4] == filename1[0:len(currentText[:-4])]:
            filename = currentText
        #autofill
        Range = selfMT.dlg.Range.text()
        if Range == "":
            Range = "100"
        autoPoly = filename[:-4]+"_synthetic_"+Range+".img"
        selfMT.dlg.Output_lineEdit.setText(autoPoly)
        if os.path.exists(autoPoly):
            selfMT.dlg.exists2.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists2.setText("")

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

    def run(self):
        import random
        import marinetools
        from marinetools.Synthetic.Synthetic_dialog import SyntheticDialog
        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 = SyntheticDialog()
            selfMT = self
            self.dlg.pushButton_1.clicked.connect(Synthetic.select_input_file1) 
            self.dlg.pushButton_2.clicked.connect(Synthetic.select_input_file2) 
            self.dlg.pushButton_3.clicked.connect(Synthetic.select_input_file3) 
            self.dlg.BathyCombo1.currentIndexChanged.connect(Synthetic.indexChanged)
            self.dlg.TracklineCombo1.currentIndexChanged.connect(Synthetic.VarDepChanged)
            self.dlg.Range.textChanged.connect(Synthetic.indexChanged) 
            self.dlg.ConstDep.textChanged.connect(Synthetic.ConstDepChanged) 
            self.dlg.ConstAlt.textChanged.connect(Synthetic.ConstAltChanged) 
            self.dlg.helpButton.clicked.connect(Synthetic.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.BathyCombo1.clear() 
        self.dlg.TracklineCombo1.clear() 
        # Populate the comboBox with names of all the loaded layer   
        self.dlg.BathyCombo1.addItems([layer.name() for layer in layers if str(layer.type())== "1" or str(layer.type())== "LayerType.Raster"])
        self.dlg.TracklineCombo1.addItems([layer.name() for layer in layers if str(layer.type())== "0" or str(layer.type())== "LayerType.Vector"])
        Synthetic.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.BathyCombo1.currentIndex()
            currentText = selfMT.dlg.BathyCombo1.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
            inputBathy = filename
            
            selectedLayerIndex = self.dlg.TracklineCombo1.currentIndex()
            currentText = selfMT.dlg.TracklineCombo1.currentText()
            layers = QgsProject.instance().mapLayers().values()
            a=0
            filename="NULL"
            for layer in (layer1 for layer1 in layers if str(layer1.type())== "0" or str(layer1.type())== "LayerType.Vector"):
                if a == selectedLayerIndex:
                    filename = str(layer.source())
                a=a+1
            if currentText not in filename:
                filename = currentText

            surveyLine = filename 

            # Test input file to find if supported 
            from marinetools.fileTest import fileTest
            #inputBathy = fileTest.main(inputBathy,"raster")
            #if inputBathy[0:5] == "Error":
            #    print(inputBathy)
            #    return
            
            OutputFile2 = self.dlg.Output_lineEdit.text()  
  
            Range = self.dlg.Range.text()
            if Range == "":
                Range = "100"
            MaxRange = 10000000
            dict = processing.run("native:rasterlayerproperties", {'INPUT':inputBathy,'BAND':1})
            Resolution = float(dict['PIXEL_WIDTH'])

            newdir = str(os.path.dirname(inputBathy) + "/tempMT")
            if not os.path.exists(newdir):
                os.mkdir(newdir)

            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
            
            alphabet = 'ZYXWVUTSRQPONMLKJIHGFEDCBA0987654321'
            rand = alphabet[random.randint(1,25)] + alphabet[random.randint(1,35)] + alphabet[random.randint(1,35)] + alphabet[random.randint(1,35)]

            VarDepth = self.dlg.comboBox.currentText()
            ConstAlt = self.dlg.ConstAlt.text()
            ConstDep = self.dlg.ConstDep.text()
            
            #Just find the crs
            result = processing.run("native:buffer", {'INPUT':surveyLine,'DISTANCE':Range,'SEGMENTS':5,
                                                      'END_CAP_STYLE':0,'JOIN_STYLE':0,'MITER_LIMIT':2,'DISSOLVE':True,
                                                      'SEPARATE_DISJOINT':False,'OUTPUT':'TEMPORARY_OUTPUT'})
            bufferVector = result['OUTPUT']
            crs = bufferVector.crs()
            RangeDeg = Range
            ResolutionDeg = Resolution
            if str(crs) == "<QgsCoordinateReferenceSystem: EPSG:4326>":
                RangeDeg = str(float(Range)/111000.0)
                ResolutionDeg = str(float(Resolution)/111000.0)

            # Find if points or line
            result = processing.run("gdal:ogrinfo", {'INPUT':surveyLine,'ALL_LAYERS':False,'SUMMARY_ONLY':True,
                                                    'NO_METADATA':False,'EXTRA':'','OUTPUT':'TEMPORARY_OUTPUT'})
            f= open(result['OUTPUT'])
            oneliner = f.read()
            f.close()
            lines = oneliner.splitlines()
            for line in lines:
                vectorTypePos = line.find("Geometry:")
                if vectorTypePos != -1:
                    vectorType = line[vectorTypePos+10:]
            if vectorType == "Line String":
                result = processing.run("native:pointsalonglines", {'INPUT':surveyLine,'DISTANCE':ResolutionDeg,'START_OFFSET':0,
                                                                    'END_OFFSET':0,'OUTPUT':'TEMPORARY_OUTPUT'})
                linePoints = result['OUTPUT']
                surveyLine = linePoints

            #convert bathymetry to points (VECTOR POINTS)
            result = processing.run("native:pixelstopoints", {'INPUT_RASTER':inputBathy,'RASTER_BAND':1,'FIELD_NAME':'depth',
                                                              'OUTPUT':'TEMPORARY_OUTPUT'})
            allPoints = result['OUTPUT']

            # Get slope and azimuth (aspect) of bathymetry
            result = processing.run("grass:r.slope.aspect", {'elevation':inputBathy,'format':0,'precision':2,'-a':True,
                                                             '-e':False,'-n':True,'zscale':1,'min_slope':0,
                                                             'slope':'TEMPORARY_OUTPUT','aspect':'TEMPORARY_OUTPUT',
                                                             'GRASS_REGION_PARAMETER':None,'GRASS_REGION_CELLSIZE_PARAMETER':0,
                                                             'GRASS_RASTER_FORMAT_OPT':'','GRASS_RASTER_FORMAT_META':''})
            slopeMap = result['slope']
            aspectMap = result['aspect']
            
            tempfile1 = newdir + "/tempfile1_"+rand+".shp"
            tempfile2 = newdir + "/tempfile2_"+rand+".shp"
            tempfile3 = newdir + "/tempfile3_"+rand+".shp"
            tempfile4 = newdir + "/tempfile4_"+rand+".shp"
 
            stats = processing.run("native:rasterlayerstatistics", {'INPUT':str(inputBathy),'BAND':1})
            InputMean = stats["MEAN"]
            InputStd = stats["STD_DEV"]
            pointsGap = abs(InputMean)- abs(InputStd)
            RangeSquared = str(float(Range)*float(Range))
            result = processing.run("native:rastersampling", {'INPUT':surveyLine,'RASTERCOPY':inputBathy,
                                                              'COLUMN_PREFIX':'depth','OUTPUT':'TEMPORARY_OUTPUT'})
            addedDepths = result['OUTPUT']

            # Points given for line VarDepth has Fish depth values
            if VarDepth != "":
                GroundRange = 'sqrt('+RangeSquared+'-((-"depth1"-'+VarDepth+')*(-"depth1"-'+VarDepth+')))'
                expression = '-"depth1" > ' + VarDepth + ' AND -"depth1" - '+ VarDepth +' < '+ Range
                if InputMean > 0.0:
                    GroundRange = 'sqrt('+RangeSquared+'-(("depth1"-'+VarDepth+'")*("depth1"-'+VarDepth+')))'
                    expression = '"depth1" > ' + VarDepth + ' AND "depth1" - '+ VarDepth +' < '+ Range
                result = processing.run("native:extractbyexpression", {'INPUT':addedDepths,'EXPRESSION':expression,
                                                                       'OUTPUT':'TEMPORARY_OUTPUT'})
                selectedAbove = result['OUTPUT']
                addedDepths = selectedAbove
            # Points given for line Constant Fish depth values
            if ConstDep != "":
                GroundRange = 'sqrt('+RangeSquared+'-((-"depth1"-'+ConstDep+')*(-"depth1"-'+ConstDep+')))'
                expression = '-"depth1" > ' + ConstDep + ' AND -"depth1" - '+ ConstDep +' < '+ Range
                if InputMean > 0.0:
                    GroundRange = 'sqrt('+RangeSquared+'-(("depth1"-'+ConstDep+')*("depth1"-"'+ConstDep+'")))'
                    expression = '"depth1" > ' + ConstDep + ' AND "depth1" - '+ ConstDep +' < '+ Range
                result = processing.run("native:extractbyexpression", {'INPUT':addedDepths,'EXPRESSION':expression,
                                                                       'OUTPUT':'TEMPORARY_OUTPUT'})
                selectedAbove = result['OUTPUT']
                addedDepths = selectedAbove
            # Points given for constant  Fish altitude
            if ConstAlt != "":
                AltSquared = str(float(ConstAlt)*float(ConstAlt))
                GroundRange = str(sqrt(float(RangeSquared)-float(AltSquared)))
            #print(GroundRange)
            if str(crs) == "<QgsCoordinateReferenceSystem: EPSG:4326>":
                GroundRange = GroundRange + '/111000.0'
            #print(GroundRange)
            result = processing.run("native:fieldcalculator", {'INPUT':addedDepths,'FIELD_NAME':'coverage','FIELD_TYPE':0,
                                                               'FIELD_LENGTH':0,'FIELD_PRECISION':0,
                                                               'FORMULA':GroundRange,'OUTPUT':tempfile1})
            #addedGroundRange = result['OUTPUT']

            # Check we haev some points!
            result = processing.run("gdal:ogrinfo", {'INPUT':tempfile1,'ALL_LAYERS':False,'SUMMARY_ONLY':True,
                                                    'NO_METADATA':False,'EXTRA':'','OUTPUT':'TEMPORARY_OUTPUT'})

            f= open(result['OUTPUT'])
            oneliner = f.read()
            f.close()
            lines = oneliner.splitlines()
            for line in lines:
                featureCountPos = line.find("Feature Count:")
                if featureCountPos != -1:
                    featureCount = line[featureCountPos+15:]
            #print(str(featureCount))
            if int(featureCount) < 2:
                self.loading_screen.stop_animation()
                # Create an informational message box
                msgbox = QMessageBox(QMessageBox.Information,"Synthetic Imagery Output",
                                     "There is not any imagery from this trackline",QMessageBox.Ok)
                msgbox.exec_() # Display the message box
                return
                            

            addedGroundRange = tempfile1
            result = processing.run("native:buffer", {'INPUT':addedGroundRange,'DISTANCE':QgsProperty.fromExpression('"coverage"'),
                                             'SEGMENTS':5,'END_CAP_STYLE':0,'JOIN_STYLE':0,'MITER_LIMIT':2,
                                             'DISSOLVE':True,'OUTPUT':'TEMPORARY_OUTPUT'})
            stencilClip = result['OUTPUT']
            
            #clip points by buffer (CLIPPED)
            result = processing.run("native:clip", {'INPUT':allPoints,'OVERLAY':stencilClip,'OUTPUT':'TEMPORARY_OUTPUT'})
            clippedPoints = result['OUTPUT']
            
            #create lines of distance to line (SHORTEST LINES)
            result = processing.run("native:shortestline", {'SOURCE':clippedPoints,'DESTINATION':addedGroundRange,'METHOD':0,
                                                            'NEIGHBORS':1,'DISTANCE':None,'OUTPUT':'TEMPORARY_OUTPUT'})
            shortestLine = result['OUTPUT']
            #shortestLine = tempfile4
            

            #create points at beginning and end distance set to more than range e.g. 10000000 (INTERPOLATED POINTS)
            result = processing.run("native:pointsalonglines", {'INPUT':shortestLine,'DISTANCE':MaxRange,'START_OFFSET':0,
                                                                'END_OFFSET':0,'OUTPUT':'TEMPORARY_OUTPUT'})
            distancePoints = result['OUTPUT']
            result = processing.run("native:rastersampling", {'INPUT':distancePoints,'RASTERCOPY':aspectMap,
                                                              'COLUMN_PREFIX':'aspectXX','OUTPUT':'TEMPORARY_OUTPUT'})
            withAspectPoints = result['OUTPUT']
            result = processing.run("native:rastersampling", {'INPUT':withAspectPoints,'RASTERCOPY':slopeMap,
                                                              'COLUMN_PREFIX':'slopeXXX','OUTPUT':'TEMPORARY_OUTPUT'})
            withSlopePoints = result['OUTPUT']
            
            # Lambert's Law?
            
            #expression = 'atan( (- "depth" - "FishDepth" )/ "distance" ) * 57.29577' # alpha - This is the insonification angle
            expression = ' "slopeXXX1"  *  cos( ( "aspectXX1"  -  "angle" ) * 0.0174533 )' # theta - This is the slope relative to the nearest trackline
            #expression = ' cos( ( ("slopeXXX1"  *  cos( ( "aspectXX1"  -  "angle" ) * 0.0174533 ))) * 0.0174533)' # lambert value - This is the cos(90 - theta)
            #expression = ' ( ((atan( (- "depth" - "FishDepth" )/ "distance" ) * 57.29577) + (("slopeXXX1"  *  cos( ( "aspectXX1"  -  "angle" )) * 0.0174533 ))) * 0.0174533)' # lambert value - This is the cos(90 - theta)
            result = processing.run("native:fieldcalculator", {'INPUT':withSlopePoints,'FIELD_NAME':'theta','FIELD_TYPE':0,
                                                       'FIELD_LENGTH':10,'FIELD_PRECISION':4,'FORMULA':expression,
                                                       'OUTPUT':'TEMPORARY_OUTPUT'})
            StartlinePoints = result['OUTPUT']
            
            #StartlinePoints = tempfile3
            
            # Rasterize points to raster (RASTERIZED)
            # distance to trackline and angle to trackline or any fieldname in original trackline points
            #result = processing.run("gdal:rasterize", {'INPUT':StartlinePoints,'FIELD':'distance','BURN':0,'USE_Z':False,
            #result = processing.run("gdal:rasterize", {'INPUT':StartlinePoints,'FIELD':'angle','BURN':0,'USE_Z':False,
            #angleRaster = result['OUTPUT']
            result = processing.run("gdal:rasterize", {'INPUT':StartlinePoints,'FIELD':'theta','BURN':0,'USE_Z':False,
                                                       'UNITS':1,'WIDTH':Resolution,'HEIGHT':Resolution,'EXTENT':None,'NODATA':0,'OPTIONS':None,
                                                       'DATA_TYPE':5,'INIT':None,'INVERT':False,'EXTRA':'','OUTPUT':'TEMPORARY_OUTPUT'})
            thetaValues = result['OUTPUT']

            # Clip and put to OutputFile2
            rseult = processing.run("gdal:cliprasterbymasklayer", {'INPUT':thetaValues,'MASK':stencilClip,'SOURCE_CRS':None,
                                                                   'TARGET_CRS':None,'TARGET_EXTENT':None,'NODATA':None,
                                                                   'ALPHA_BAND':False,'CROP_TO_CUTLINE':True,
                                                                   'KEEP_RESOLUTION':True,'SET_RESOLUTION':False,
                                                                   'X_RESOLUTION':None,'Y_RESOLUTION':None,
                                                                   'MULTITHREADING':False,'OPTIONS':None,'DATA_TYPE':0,
                                                                   'EXTRA':'','OUTPUT':OutputFile2})
            



            self.loading_screen.stop_animation()

            fname = os.path.dirname(str(OutputFile2))
            vlayer = QgsRasterLayer(str(OutputFile2), str(OutputFile2[len(fname)+1:-4]))
            QgsProject.instance().addMapLayer(vlayer)
            
            pass


































