# -*- coding: utf-8 -*-

"""
/***************************************************************************
 QProto
                                 A QGIS plugin
 QProto algorithm
 Generated by Plugin Builder: http://g-sherman.github.io/Qgis-Plugin-Builder/
                              -------------------
        begin                : 2018-04-19
        copyright            : (C) 2018 by Matteo Ghetta (Faunalia)
        email                : matteo.ghetta@faunalia.eu
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   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.                                   *
 *                                                                         *
 ***************************************************************************/
"""

__author__ = 'Matteo Ghetta (Faunalia)'
__date__ = '2018-04-19'
__copyright__ = '(C) 2018 by Matteo Ghetta (Faunalia)'

# This will get replaced with a git SHA1 when you do a git archive

__revision__ = '$Format:%H$'

from qgis.PyQt.QtCore import QCoreApplication, QVariant
from qgis.PyQt.QtGui import QIcon
from qgis.core import (QgsProcessing,
                       QgsProcessingAlgorithm,
                       QgsRasterLayer,
                       QgsPointXY,
                       QgsGeometry,
                       QgsProcessingParameterFeatureSource,
                       QgsProcessingParameterRasterLayer,
                       QgsProcessingParameterField,
                       QgsProcessingParameterNumber,
                       QgsProcessingParameterFolderDestination,
                       QgsVectorLayer,
                       QgsVectorFileWriter,
                       QgsField,
                       QgsFields,
                       QgsFeature,
                       QgsWkbTypes,
                       QgsDistanceArea,
                       QgsCoordinateTransformContext,
                       QgsSpatialIndex,
                       QgsProcessingException,
                       QgsProcessingUtils)

import processing
import os
import math
from osgeo import gdal
import numpy as np
from datetime import datetime
import time
import tempfile

from .utils import (map2pixel,
                    triangolo,
                    pixel2map,
                    OutStyle)

pluginPath = os.path.dirname(__file__)

class QProtoAlgorithm(QgsProcessingAlgorithm):

    INPUT = 'INPUT'
    DTM = 'DTM'
    ID_FIELD = 'ID_FIELD'
    ELEV_FIELD = 'ELEV_FIELD'
    ASPECT_FIELD = 'ASPECT_FIELD'
    WEIGHT_FIELD = 'WEIGHT_FIELD'
    MASS_FIELD = 'MASS_FIELD'
    SHADOW_FIELD = 'SHADOW_FIELD'
    FIX_SHADOW = 'FIX_SHADOW'
    LATERAL_FIELD = 'LATERAL_FIELD'
    FIX_LATERAL = 'FIX_LATERAL'
    VISIBILITY = 'VISIBILITY'
    OUTPUT = 'OUTPUT'

    def name(self):
        return 'QProto'

    def displayName(self):
        return self.tr(self.name())

    def group(self):
        return self.tr(self.groupId())

    def groupId(self):
        return 'QProto'

    def tr(self, string):
        return QCoreApplication.translate('Processing', string)

    def createInstance(self):
        return QProtoAlgorithm()

    def tags(self):
        return self.tr("qproto, rock, slope").split(",")

    def icon(self):
        return QIcon(os.path.join(pluginPath, "icon.svg"))

    def shortHelpString(self):
        return self.tr(
        '''

<p>QPROTO, based on the cone method developed by Jaboyedoff and Labiouse (2011), identifies the areas that are most exposed to a rockfall phenomenon only considering the topography of the slope and the calibration of some empirical parameters. Furthermore, it allows the quantitative estimation of block velocity and kinetic energy and the evaluation of a time-independent rockfall relative hazard.</p>

<p><u>Digital Elevation Model:</u> Raster layer of investigated area [m]</p>

<p><u>Input Point Layer:</u> vector layer of the source points, containing the following attributes:</p>

<p style="margin-left: 40px;"><u>ID:</u> Identifying number of the source point;</p>

<p style="margin-left: 40px;"><u>Elevation:</u> Elevation of the source point [m];</p>

<p style="margin-left: 40px;"><u>Aspect:</u> Dip direction of the slope in the source point, ranging between 0-360 [&deg;]. It defines the orientation of the cone;</p>

<p style="margin-left: 40px;"><u>Detachment Propensity:</u>&nbsp;Propensity index to detachment associated to each source point, ranging between 0-1 [-];</p>

<p style="margin-left: 40px;"><u>Boulder mass:</u> Boulder mass [kg];</p>

<p style="margin-left: 40px;"><u>Energy Line Angle:</u> Energy line angle of the cone having apex in the source point, ranging between 0-90 [&deg;]. Recommended values between 20&deg;-50&deg;. Optionally, a fixed value can be used for quick analysis. In this case the value has to be indicated in <em>Fixed value of the Energy Line Angle</em>.</p>

<p style="margin-left: 40px;"><u>Lateral spreading angle:</u> Lateral dispersion angle of the cone having apex in the source point, ranging between 0-90 [&deg;]. Recommended values between 0&deg;-30&deg;. This angle is added and subtracted to the dip direction of the slope (aspect) to define the limits of the visibility cone in the horizontal plane. Optionally, a fixed value can be used for quick analysis. In this case the value has to be indicated in <em>Fixed value of the Lateral Spreading Angle</em>.</p>

<p><u>Fixed value of the Energy Line Angle:</u> Empty if Energy Line Angle is considered as an attribute of each source point [&deg;]</p>

<p><u>Fixed value of the Lateral Spreading Angle:</u> Empty if Lateral Dispersion Angle is considered as an attribute of the point [&deg;].</p>

<p><u>Visibility Distance:</u> Distance to which the analysis can be extended [m]</p>

<p>&nbsp;</p>

<p><u>Help:</u> <a href="https://gitlab.com/faunalia/QPROTO/blob/master/help/QPROTO_User_Manual_1.4.pdf">https://gitlab.com/faunalia/QPROTO/blob/master/help/QPROTO_User_Manual_1.4.pdf</a></p>

<p><u>Training dataset:</u> <a href="https://gitlab.com/faunalia/QPROTO/raw/master/help/trainig_dataset.zip">https://gitlab.com/faunalia/QPROTO/tree/master/help</a></p>

<p><u>info:</u> <a href="mailto:marta.castelli@polito.it">marta.castelli@polito.it</a></p>



        '''
        )

    def helpUrl(self):
        return "https://gitlab.com/faunalia/QPROTO/blob/master/help/QPROTO_User_Manual_1.4.pdf"

    def initAlgorithm(self, config):
        '''
        define all the input parameters
        '''

        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.INPUT,
                self.tr('Input point layer'),
                [QgsProcessing.TypeVectorPoint]
            )
        )

        self.addParameter(
            QgsProcessingParameterRasterLayer(
                self.DTM,
                self.tr('Digital Elevation Model')
            )
        )

        self.addParameter(
            QgsProcessingParameterField(
                self.ID_FIELD,
                self.tr('ID'),
                'ID_IN',
                self.INPUT
            )
        )

        self.addParameter(
            QgsProcessingParameterField(
                self.ELEV_FIELD,
                self.tr('Elevation [m]'),
                'QUOTA_IN',
                self.INPUT
            )
        )

        self.addParameter(
            QgsProcessingParameterField(
                self.ASPECT_FIELD,
                self.tr('Aspect [°]'),
                'ASPECT',
                self.INPUT
            )
        )

        self.addParameter(
            QgsProcessingParameterField(
                self.WEIGHT_FIELD,
                self.tr('Detachment propensity [-]'),
                'PESO',
                self.INPUT
            )
        )

        self.addParameter(
            QgsProcessingParameterField(
                self.MASS_FIELD,
                self.tr('Boulder mass [kg]'),
                'MASSA',
                self.INPUT
            )
        )

        self.addParameter(
            QgsProcessingParameterField(
                self.SHADOW_FIELD,
                self.tr('Energy line angle [°]'),
                None,
                self.INPUT,
                optional=True
            )
        )

        self.addParameter(
            QgsProcessingParameterNumber(
                self.FIX_SHADOW,
                self.tr('Fixed value for energy line angle [°]'),
                QgsProcessingParameterNumber.Double,
                30,
                optional=True
            )
        )

        self.addParameter(
            QgsProcessingParameterField(
                self.LATERAL_FIELD,
                self.tr('Lateral spreading angle [°]'),
                None,
                self.INPUT,
                optional=True
            )
        )

        self.addParameter(
            QgsProcessingParameterNumber(
                self.FIX_LATERAL,
                self.tr('Fixed value for lateral spreading angle [°]'),
                QgsProcessingParameterNumber.Double,
                10,
                optional=True
            )
        )

        self.addParameter(
            QgsProcessingParameterNumber(
                self.VISIBILITY,
                self.tr('Visibility distance [m]'),
                QgsProcessingParameterNumber.Double,
                800
            )
        )

        self.addParameter(
            QgsProcessingParameterFolderDestination(
                self.OUTPUT,
                self.tr('Output Folder')
            )
        )

    def processAlgorithm(self, parameters, context, feedback):
        """
        Here is where the processing itself takes place.
        """

        # start time of the algorithm
        start_time = datetime.now()

        # get all the inputs
        layer = self.parameterAsSource(parameters, self.INPUT, context)
        raster = self.parameterAsRasterLayer(parameters, self.DTM, context)
        ID_Field = self.parameterAsFields(parameters, self.ID_FIELD, context)
        Elevation_Field = self.parameterAsFields(parameters, self.ELEV_FIELD, context)
        Aspect_Field = self.parameterAsFields(parameters, self.ASPECT_FIELD, context)
        Weight_Field = self.parameterAsFields(parameters, self.WEIGHT_FIELD, context)
        Mass_Field = self.parameterAsFields(parameters, self.MASS_FIELD, context)
        Shadow_Angle_Field = self.parameterAsFields(parameters, self.SHADOW_FIELD, context)
        Fixed_value_for_the_Shadow_Angle = self.parameterAsDouble(parameters, self.FIX_SHADOW, context)
        Lateral_Dispersion_Field = self.parameterAsFields(parameters, self.LATERAL_FIELD, context)
        Fixed_value_for_the_Lateral_Dispersion = self.parameterAsDouble(parameters, self.FIX_LATERAL, context)
        Visibility_Distance = self.parameterAsDouble(parameters, self.VISIBILITY, context)

        # set the output folder
        Results = self.parameterAsString(parameters, self.OUTPUT, context)

        feedback.pushConsoleInfo(str(parameters['OUTPUT']))

        # set the output folder
        Results = self.parameterAsString(parameters, self.OUTPUT, context)

        # get the CRS from the source layer
        crs = layer.sourceCrs()
        # get the extent from the raster
        ex = raster.extent()
        extent = '{}, {}, {},{}'.format(ex.xMinimum(),
                                        ex.xMaximum(),
                                        ex.yMinimum(),
                                        ex.yMaximum())

        # read the dtm as a numpy array
        rast_ds = gdal.Open(raster.source())
        geoTransform = rast_ds.GetGeoTransform()
        Output_Resolution = geoTransform[1]
        rast_array = np.array(rast_ds.GetRasterBand(1).ReadAsArray())

        v_layer = QgsVectorLayer('Point?crs=' + crs.toWkt(), 'points', "memory") # memory layer
        pr = v_layer.dataProvider()
        pr.addAttributes([
            QgsField('ID_FIN', QVariant.Int),
            QgsField('ELEV_FIN', QVariant.Int),
            QgsField('ID_IN', QVariant.Int),
            QgsField('ELEV_IN', QVariant.Int),
            QgsField('SHAD_ANGLE', QVariant.Int),
            QgsField('LAT_DISP', QVariant.Int),
            QgsField('WEIGHT', QVariant.Double),
            QgsField('V', QVariant.Double),
            QgsField('E', QVariant.Double),
            QgsField('W_EN', QVariant.Double)
        ])

        v_layer.updateFields()

        # Creation of the final vector containing the points whitin the rockfall propagation area (this vector is saved in the output directory).
        fields = QgsFields()
        fields.append(QgsField('ID_FIN', QVariant.Int)) # adding of the 'ID_FIN' field
        fields.append(QgsField('ELEV_FIN', QVariant.Int)) # adding of the 'ELEV_FIN' field
        fields.append(QgsField('ID_IN', QVariant.Int)) # adding of the 'ID_IN' field
        fields.append(QgsField('ELEV_IN', QVariant.Int)) # adding of the 'ELEV_IN' field
        fields.append(QgsField('SHAD_ANGLE', QVariant.Int)) # adding of the 'SHAD_ANGLE' field
        fields.append(QgsField('LAT_DISP', QVariant.Int)) # adding of the 'LAT_DISP' field
        fields.append(QgsField('WEIGHT', QVariant.Double)) # adding of the 'WEIGHT' field
        fields.append(QgsField('V', QVariant.Double)) # adding of the 'V' field
        fields.append(QgsField('E', QVariant.Double)) # adding of the 'E' field
        fields.append(QgsField('W_EN', QVariant.Double)) # adding of the 'W_EN' field
        fields.append(QgsField('COUNT', QVariant.Int)) # adding of the 'COUNT' field
        fields.append(QgsField('SUM_WEIGHT', QVariant.Double)) # adding of the 'SUM_WEIGHT' field
        fields.append(QgsField('W_TOT_EN', QVariant.Double)) # adding of the 'W_TOT_EN' field
        pathh = os.path.join(Results, 'finalpoints.shp')
        # write the final shapefile using native method and not Processing one
        finalpoints = QgsVectorFileWriter(pathh, None, fields, QgsWkbTypes.Point, crs, "ESRI Shapefile")
        totalfeatures = [] # this list contains all the features within the rockfall propagation area


        # Counting of the source points (it is necessary to evaluate the progress of the analysis)
        feat_count = layer.featureCount()
        feedback.setProgressText('Total number of source points to process: {}'.format(str(feat_count)))

        i = 1
        coni_creati = 0
        id = 0

        dist_calc = QgsDistanceArea()
        dist_calc.setSourceCrs(crs, QgsCoordinateTransformContext())
        dist_calc.setEllipsoid(crs.ellipsoidAcronym())

        # Check for proper field types
        input_fields = [ID_Field, Elevation_Field, Aspect_Field, Weight_Field, Mass_Field]
        if Shadow_Angle_Field:
            input_fields.append(Shadow_Angle_Field)
        if Lateral_Dispersion_Field:
            input_fields.append(Lateral_Dispersion_Field)


        for f in layer.getFeatures():

            if not Shadow_Angle_Field:
                shad_ang = Fixed_value_for_the_Shadow_Angle
            else:
                shad_ang = f[Shadow_Angle_Field[0]]

            if not Lateral_Dispersion_Field:
                lat_disp = Fixed_value_for_the_Lateral_Dispersion
            else:
                lat_disp = f[Lateral_Dispersion_Field[0]]

            message = 'Processing souce point No. {} of {} ...'.format(str(i), str(feat_count))
            feedback.setProgressText(message)
            feedback.setProgress((i/float(feat_count)) * 89)
            i = i + 1
            if feedback.isCanceled():
                break



            # Creation of a triangle that limits the lateral dispersion of the phenomenon (the triangle is created as geometry)
            point1 = f.geometry().asPoint() # Point 1 (it coincides with the i-th source point)
            tr = triangolo(point1, f, Aspect_Field[0], Visibility_Distance, lat_disp)
            tr_geom = tr.geometry()
            tr_ex = tr_geom.boundingBox()
            tr_extent = '{},{},{},{}'.format(
                tr_ex.xMinimum() - Output_Resolution,
                tr_ex.xMaximum() + Output_Resolution,
                tr_ex.yMinimum() - Output_Resolution,
                tr_ex.yMaximum() + Output_Resolution
            )

            # Execution of "r.viewshed" for the i-th source point
            coordStr = '{},{}'.format(point1.x(), point1.y())
            # view_path = os.path.join(Results, 'viewshed.tif')
            view_path = os.path.join(QgsProcessingUtils.tempFolder(), 'viewshed.tif')
            viewshed = processing.run('grass7:r.viewshed', {
                'input': raster,
                'coordinates': coordStr,
                'observer_elevation': 0,
                'target_elevation': 0,
                'max_distance': Visibility_Distance,
                'refraction_coeff': 0.14286,
                'memory': 500,
                '-c':False,
                '-r':False,
                '-b':False,
                '-e':False,
                'output':'memory:',
                'GRASS_REGION_PARAMETER': tr_extent,
                'GRASS_REGION_CELLSIZE_PARAMETER':0,
                'GRASS_RASTER_FORMAT_OPT':'',
                'GRASS_RASTER_FORMAT_META':'',
                'output': view_path
            }, context=context, feedback=feedback)

            # Reading of the "r.viewshed" output as array
            view = QgsRasterLayer(viewshed['output'], 'viewshed','gdal')

            view_ds = gdal.Open(viewshed['output'])
            viewGeoTransform = view_ds.GetGeoTransform()

            view_rast_array = np.array(view_ds.GetRasterBand(1).ReadAsArray())
            # The "r.viewshed" output is filtered with the aim of saving all the points within the energy line
            v_ex=view.extent()

            (firstRow, firstColumn, lastRow, lastColumn) = map2pixel(viewGeoTransform, v_ex)
            below_enLine = []
            for row in range(firstRow, lastRow):
                for col in range(firstColumn, lastColumn):
                    pt_in=(view_rast_array[col, row])
                    if pt_in <= (90 - shad_ang): # if this condition is true, the point is within the energy line
                        (x, y) = pixel2map(viewGeoTransform, row, col)
                        vpoint = QgsPointXY(x,y)
                        below_enLine.append(vpoint)

            # Creation of the minimum convex hull which contains all the 'below_enLine' items
            outGeom = QgsGeometry()
            tmpGeom = QgsGeometry(outGeom.fromMultiPointXY(below_enLine))
            convex = tmpGeom.convexHull()

            # Clipping of 'convex' with the triangle
            intPol = convex.intersection(tr_geom)

            # Generation of the points within the previous result ('intPol')
            intPol_ex = intPol.boundingBox()
            (firstRow, firstColumn, lastRow, lastColumn) = map2pixel(geoTransform, intPol_ex)
            cone_creation = False
            for row in range(firstRow, lastRow + 1):
                for col in range(firstColumn, lastColumn + 1):
                    (x, y) = pixel2map(geoTransform, row, col)
                    jPoint = QgsPointXY()
                    jPoint.setX(x)
                    jPoint.setY(y)
                    if intPol.contains(jPoint):
                        quota=int(rast_array[col, row]) # elevation value of the j-th point of 'intPol'
                        distance = dist_calc.measureLine([jPoint,point1])
                        det = ((f[Elevation_Field[0]]) - quota - distance * math.tan(math.radians(shad_ang))) # from Jaboyedoff e Labiouse (2011)
                        if det >= 0: # if 'det' is positive, the current point is within the energy line
                            inAttr = [id] # adding of 'ID_FIN'
                            inAttr.append(quota) # adding of 'ELEV_FIN'
                            inAttr.append(f[ID_Field[0]]) # adding of 'ID_IN' from the i-th source point
                            inAttr.append(f[Elevation_Field[0]]) # adding of 'ELEV_IN' from the i-th source point
                            inAttr.append(shad_ang) # adding of 'SHAD_ANGLE' from the i-th source point
                            inAttr.append(lat_disp) # adding of 'LAT_DISP' from the i-th source point
                            inAttr.append(f[Weight_Field[0]]) # adding of 'WEIGHT' from the i-th source point
                            v = math.sqrt(2 * 9.81 * det) # calculation of the velocity
                            e = 0.5 * (f[Mass_Field[0]]) * (v ** 2) # calculation of the energy
                            w_en = f[Weight_Field[0]] * e
                            inAttr.append(v) # adding of 'V'
                            inAttr.append(e) # adding of 'E'
                            inAttr.append(w_en) # adding of 'W_EN'
                            outFeature = QgsFeature()
                            outFeature.setGeometry(QgsGeometry.fromPointXY(jPoint)) # output geometry
                            outFeature.setAttributes(inAttr) # output attributes
                            totalfeatures.append(outFeature)
                            id += 1
                            cone_creation = True

            if cone_creation:
                message = 'Some points added to "finalpoints.shp".'
                coni_creati += 1
                feedback.setProgressText(message)
            else:
                message = 'No points added to "finalpoints.shp".'
                feedback.setProgressText(message)

        # to remove eventual useless file created
        # for item in Results:
        #     if item.startswith("finalpoints"):
        #         os.remove(os.path.join( path, item))

        pr.addFeatures(totalfeatures)

        end_time_pr = datetime.now()
        message = '{} cones have been generated starting from {} source points'.format(coni_creati, feat_count)
        feedback.setProgressText(message)
        message = '{} points processed in (hh:mm:ss): {}'.format(feat_count, end_time_pr - start_time)
        feedback.setProgressText(message)
        feedback.setProgress(90)

        w_en = {} # this dictionary contains the field 'W_EN' of each point within the rockfall propagation area
        weights = {} # this dictionary contains the field 'WEIGHT' of each point within the rockfall propagation area
        allfeatures={} # this dictionary collects all the points within the rockfall propagation area
        for elem in v_layer.getFeatures():
            allfeatures[elem.id()] = QgsFeature(elem)
            weights[elem.id()] = elem[6]
            w_en[elem.id()] = elem[9]


        index = QgsSpatialIndex() # this index contains all the features of 'allfeatures'
        for w in allfeatures.values():
            index.insertFeature(w)

        for feat in allfeatures.values():
            inGeom = feat.geometry() # input geometry
            inAttr = feat.attributes() # input attributes
            idsList = index.intersects(inGeom.boundingBox())
            inAttr.append(len(idsList))
            inAttr.append(sum(weights[x] for x in idsList)) # sum of the 'WEIGHT' fields for each duplicate geometry
            inAttr.append(sum(w_en[y] for y in idsList)) # sum of the 'W_EN_TOT' fields for each duplicate geometry
            outFeat = QgsFeature()
            outFeat.setGeometry(inGeom) # output geometry
            outFeat.setAttributes(inAttr) # output attributes
            finalpoints.addFeature(outFeat)

        del finalpoints

        finalpoints = QgsVectorLayer(os.path.join(Results , 'finalpoints.shp'), 'finalpoints', 'ogr')
        # add the layers created from the algorithm without explicitly declared an outputs
        context.addLayerToLoadOnCompletion(finalpoints.source(), context.LayerDetails(
            name='finalpoints',
            project=context.project()
        ))

        feedback.setProgressText('Creating raster outputs ...')


        # Creation of "susceptibility.sdat" + style
        feedback.setProgressText('Creating susceptibility output ...')
        feedback.setProgress(91)
        visualName = 'susceptibility.sdat'
        susceptibility = os.path.join(Results, visualName)

        # SAGA rasterize
        saga_susceptibility = processing.run("saga:rasterize", {
            'INPUT': finalpoints.source(),
            'FIELD': "SUM_WEIGHT",
            'OUTPUT': 2,
            # 0 first, 1 last, 2 min, 3 max, 4 mean. Default is 0
            'MULTIPLE': 0,
            'LINE_TYPE': 0,
            'POLY_TYPE': 0,
            'GRID_TYPE': 4,
            'TARGET_USER_XMIN TARGET_USER_XMAX TARGET_USER_YMIN TARGET_USER_YMAX': extent,
            'TARGET_USER_SIZE': Output_Resolution,
            'TARGET_USER_FITS': 1,
            'GRID': susceptibility
        }, context=context, feedback=feedback)

        # create the raster layer from the SAGA output
        susceptibility = QgsRasterLayer(susceptibility, visualName)
        # call the OutStyle function to create the style (qml) to the layer
        sty = OutStyle(susceptibility, Results, visualName)

        # call Processing algorithm to set the final style
        processing.run("qgis:setstyleforrasterlayer",{
        'INPUT': susceptibility.source(),
        'STYLE':sty
        }, context=context, feedback=feedback)

        # add the layers created from the algorithm without explicitly declared an outputs
        context.addLayerToLoadOnCompletion(susceptibility.source(), context.LayerDetails(
            name='susceptibility',
            project=context.project()
        ))

        # Creation of "count.sdat" + style
        feedback.setProgressText('Creating count output ...')
        feedback.setProgress(92)
        visualName = 'count.sdat'
        count = os.path.join(Results, visualName)

        # SAGA rasterize
        saga_count = processing.run("saga:rasterize", {
            'INPUT': finalpoints.source(),
            'FIELD': "COUNT",
            'OUTPUT': 2,
            # 0 first, 1 last, 2 min, 3 max, 4 mean. Default is 0
            'MULTIPLE': 0,
            'LINE_TYPE': 0,
            'POLY_TYPE': 0,
            'GRID_TYPE': 4,
            'TARGET_USER_XMIN TARGET_USER_XMAX TARGET_USER_YMIN TARGET_USER_YMAX': extent,
            'TARGET_USER_SIZE': Output_Resolution,
            'TARGET_USER_FITS': 1,
            'GRID': count
        }, context=context, feedback=feedback)

        # create the raster layer from the SAGA output
        count = QgsRasterLayer(count, visualName)
        # call the OutStyle function to create the style (qml) to the layer
        sty = OutStyle(count, Results, visualName)

        # call Processing algorithm to set the final style
        processing.run("qgis:setstyleforrasterlayer",{
        'INPUT': count.source(),
        'STYLE':sty
        }, context=context, feedback=feedback)

        # add the layers created from the algorithm without explicitly declared an outputs
        context.addLayerToLoadOnCompletion(count.source(), context.LayerDetails(
            name='count',
            project=context.project()
        ))


        # Creation of "v_min.sdat" + style
        feedback.setProgressText('Creating v_min output ...')
        feedback.setProgress(93)
        visualName = 'v_min.sdat'
        v_min = os.path.join(Results, visualName)

        # SAGA rasterize
        saga_v_min = processing.run("saga:rasterize", {
            'INPUT': finalpoints.source(),
            'FIELD': "V",
            'OUTPUT': 2,
            # 0 first, 1 last, 2 min, 3 max, 4 mean. Default is 0
            'MULTIPLE': 2,
            'LINE_TYPE': 0,
            'POLY_TYPE': 0,
            'GRID_TYPE': 4,
            'TARGET_USER_XMIN TARGET_USER_XMAX TARGET_USER_YMIN TARGET_USER_YMAX': extent,
            'TARGET_USER_SIZE': Output_Resolution,
            'TARGET_USER_FITS': 1,
            'GRID': v_min
        }, context=context, feedback=feedback)

        # create the raster layer from the SAGA output
        v_min = QgsRasterLayer(v_min, visualName)
        # call the OutStyle function to create the style (qml) to the layer
        sty = OutStyle(v_min, Results, visualName)

        # call Processing algorithm to set the final style
        processing.run("qgis:setstyleforrasterlayer",{
        'INPUT': v_min.source(),
        'STYLE':sty
        }, context=context, feedback=feedback)

        # add the layers created from the algorithm without explicitly declared an outputs
        context.addLayerToLoadOnCompletion(v_min.source(), context.LayerDetails(
            name='v_min',
            project=context.project()
        ))

        # Creation of "v_mean.sdat" + style
        feedback.setProgressText('Creating v_mean output ...')
        feedback.setProgress(94)
        visualName = 'v_mean.sdat'
        v_mean = os.path.join(Results, visualName)

        # SAGA rasterize
        saga_v_mean = processing.run("saga:rasterize", {
            'INPUT': finalpoints.source(),
            'FIELD': "V",
            'OUTPUT': 2,
            # 0 first, 1 last, 2 min, 3 max, 4 mean. Default is 0
            'MULTIPLE': 4,
            'LINE_TYPE': 0,
            'POLY_TYPE': 0,
            'GRID_TYPE': 4,
            'TARGET_USER_XMIN TARGET_USER_XMAX TARGET_USER_YMIN TARGET_USER_YMAX': extent,
            'TARGET_USER_SIZE': Output_Resolution,
            'TARGET_USER_FITS': 1,
            'GRID': v_mean
        }, context=context, feedback=feedback)


        v_mean = QgsRasterLayer(v_mean, visualName)
        # call the OutStyle function to create the style (qml) to the layer
        sty = OutStyle(v_mean, Results, visualName)

        # call Processing algorithm to set the final style
        processing.run("qgis:setstyleforrasterlayer",{
        'INPUT': v_mean.source(),
        'STYLE':sty
        }, context=context, feedback=feedback)

        # add the layers created from the algorithm without explicitly declared an outputs
        context.addLayerToLoadOnCompletion(v_mean.source(), context.LayerDetails(
            name='v_mean',
            project=context.project()
        ))


        # Creation of "v_max.sdat" + style
        feedback.setProgressText('Creating v_max output ...')
        feedback.setProgress(95)
        visualName = 'v_max.sdat'
        v_max = os.path.join(Results, visualName)

        # SAGA rasterize
        saga_v_max = processing.run("saga:rasterize", {
            'INPUT': finalpoints.source(),
            'FIELD': "V",
            'OUTPUT': 2,
            # 0 first, 1 last, 2 min, 3 max, 4 mean. Default is 0
            'MULTIPLE': 3,
            'LINE_TYPE': 0,
            'POLY_TYPE': 0,
            'GRID_TYPE': 4,
            'TARGET_USER_XMIN TARGET_USER_XMAX TARGET_USER_YMIN TARGET_USER_YMAX': extent,
            'TARGET_USER_SIZE': Output_Resolution,
            'TARGET_USER_FITS': 1,
            'GRID': v_max
        }, context=context, feedback=feedback)

        # create the raster layer from the SAGA output
        v_max = QgsRasterLayer(v_max, visualName)
        # call the OutStyle function to create the style (qml) to the layer
        sty = OutStyle(v_max, Results, visualName)

        # call Processing algorithm to set the final style
        processing.run("qgis:setstyleforrasterlayer",{
        'INPUT': v_max.source(),
        'STYLE':sty
        }, context=context, feedback=feedback)

        # add the layers created from the algorithm without explicitly declared an outputs
        context.addLayerToLoadOnCompletion(v_max.source(), context.LayerDetails(
            name='v_max',
            project=context.project()
        ))


        # Creation of "e_min.sdat" + style
        feedback.setProgressText('Creating e_min output ...')
        feedback.setProgress(96)
        visualName = 'e_min.sdat'
        e_min = os.path.join(Results, visualName)

        # SAGA rasterize
        saga_e_min = processing.run("saga:rasterize", {
            'INPUT': finalpoints.source(),
            'FIELD': "E",
            'OUTPUT': 2,
            # 0 first, 1 last, 2 min, 3 max, 4 mean. Default is 0
            'MULTIPLE': 2,
            'LINE_TYPE': 0,
            'POLY_TYPE': 0,
            'GRID_TYPE': 4,
            'TARGET_USER_XMIN TARGET_USER_XMAX TARGET_USER_YMIN TARGET_USER_YMAX': extent,
            'TARGET_USER_SIZE': Output_Resolution,
            'TARGET_USER_FITS': 1,
            'GRID': e_min
        }, context=context, feedback=feedback)

        # create the raster layer from the SAGA output
        e_min = QgsRasterLayer(e_min, visualName)
        # call the OutStyle function to create the style (qml) to the layer
        sty = OutStyle(e_min, Results, visualName)

        # call Processing algorithm to set the final style
        processing.run("qgis:setstyleforrasterlayer",{
        'INPUT': e_min.source(),
        'STYLE':sty
        }, context=context, feedback=feedback)

        # add the layers created from the algorithm without explicitly declared an outputs
        context.addLayerToLoadOnCompletion(e_min.source(), context.LayerDetails(
            name='e_min',
            project=context.project()
        ))


        # Creation of "e_mean.sdat" + style
        feedback.setProgressText('Creating e_mean output ...')
        feedback.setProgress(97)
        visualName = 'e_mean.sdat'
        e_mean = os.path.join(Results, visualName)

        # SAGA rasterize
        saga_e_mean = processing.run("saga:rasterize", {
            'INPUT': finalpoints.source(),
            'FIELD': "E",
            'OUTPUT': 2,
            # 0 first, 1 last, 2 min, 3 max, 4 mean. Default is 0
            'MULTIPLE': 4,
            'LINE_TYPE': 0,
            'POLY_TYPE': 0,
            'GRID_TYPE': 4,
            'TARGET_USER_XMIN TARGET_USER_XMAX TARGET_USER_YMIN TARGET_USER_YMAX': extent,
            'TARGET_USER_SIZE': Output_Resolution,
            'TARGET_USER_FITS': 1,
            'GRID': e_mean
        }, context=context, feedback=feedback)

        # create the raster layer from the SAGA output
        e_mean = QgsRasterLayer(e_mean, visualName)
        # call the OutStyle function to create the style (qml) to the layer
        sty = OutStyle(e_mean, Results, visualName)

        # call Processing algorithm to set the final style
        processing.run("qgis:setstyleforrasterlayer",{
        'INPUT': e_mean.source(),
        'STYLE':sty
        }, context=context, feedback=feedback)

        # add the layers created from the algorithm without explicitly declared an outputs
        context.addLayerToLoadOnCompletion(e_mean.source(), context.LayerDetails(
            name='e_mean',
            project=context.project()
        ))


        # Creation of "e_max.sdat" + style
        feedback.setProgressText('Creating e_max output ...')
        feedback.setProgress(98)
        visualName = 'e_max.sdat'
        e_max = os.path.join(Results, visualName)

        # SAGA rasterize
        saga_e_max = processing.run("saga:rasterize", {
            'INPUT': finalpoints.source(),
            'FIELD': "E",
            'OUTPUT': 2,
            # 0 first, 1 last, 2 min, 3 max, 4 mean. Default is 0
            'MULTIPLE': 3,
            'LINE_TYPE': 0,
            'POLY_TYPE': 0,
            'GRID_TYPE': 4,
            'TARGET_USER_XMIN TARGET_USER_XMAX TARGET_USER_YMIN TARGET_USER_YMAX': extent,
            'TARGET_USER_SIZE': Output_Resolution,
            'TARGET_USER_FITS': 1,
            'GRID': e_max
        }, context=context, feedback=feedback)

        # create the raster layer from the SAGA output
        e_max = QgsRasterLayer(e_max, visualName)
        # call the OutStyle function to create the style (qml) to the layer
        sty = OutStyle(e_max, Results, visualName)

        # call Processing algorithm to set the final style
        processing.run("qgis:setstyleforrasterlayer",{
        'INPUT': e_max.source(),
        'STYLE':sty
        }, context=context, feedback=feedback)

        # add the layers created from the algorithm without explicitly declared an outputs
        context.addLayerToLoadOnCompletion(e_max.source(), context.LayerDetails(
            name='e_max',
            project=context.project()
        ))


        # Creation of "w_en.sdat" + style
        feedback.setProgressText('Creating w_en output ...')
        feedback.setProgress(99)
        visualName = 'w_en.sdat'
        w_en = os.path.join(Results, visualName)

        # SAGA rasterize
        saga_w_en = processing.run("saga:rasterize", {
            'INPUT': finalpoints.source(),
            'FIELD': "W_EN",
            'OUTPUT': 2,
            # 0 first, 1 last, 2 min, 3 max, 4 mean. Default is 0
            'MULTIPLE': 3,
            'LINE_TYPE': 0,
            'POLY_TYPE': 0,
            'GRID_TYPE': 4,
            'TARGET_USER_XMIN TARGET_USER_XMAX TARGET_USER_YMIN TARGET_USER_YMAX': extent,
            'TARGET_USER_SIZE': Output_Resolution,
            'TARGET_USER_FITS': 1,
            'GRID': w_en
        }, context=context, feedback=feedback)

        # create the raster layer from the SAGA output
        w_en = QgsRasterLayer(w_en, visualName)
        # call the OutStyle function to create the style (qml) to the layer
        sty = OutStyle(w_en, Results, visualName)

        # call Processing algorithm to set the final style
        processing.run("qgis:setstyleforrasterlayer",{
        'INPUT': w_en.source(),
        'STYLE':sty
        }, context=context, feedback=feedback)

        # add the layers created from the algorithm without explicitly declared an outputs
        context.addLayerToLoadOnCompletion(w_en.source(), context.LayerDetails(
            name='w_en',
            project=context.project()
        ))


        # Creation of "w_tot_en.sdat" + style
        feedback.setProgressText('Creating w_tot_en output ...')
        feedback.setProgress(100)
        visualName = 'w_tot_en.sdat'
        w_tot_en =os.path.join(Results, visualName)

        # SAGA rasterize
        saga_w_tot_en = processing.run("saga:rasterize", {
            'INPUT': finalpoints.source(),
            'FIELD': "W_TOT_EN",
            'OUTPUT': 2,
            # 0 first, 1 last, 2 min, 3 max, 4 mean. Default is 0
            'MULTIPLE': 3,
            'LINE_TYPE': 0,
            'POLY_TYPE': 0,
            'GRID_TYPE': 4,
            'TARGET_USER_XMIN TARGET_USER_XMAX TARGET_USER_YMIN TARGET_USER_YMAX': extent,
            'TARGET_USER_SIZE': Output_Resolution,
            'TARGET_USER_FITS': 1,
            'GRID': w_tot_en
        }, context=context, feedback=feedback)

        # create the raster layer from the SAGA output
        w_tot_en = QgsRasterLayer(w_tot_en, visualName)
        # call the OutStyle function to create the style (qml) to the layer
        sty = OutStyle(w_tot_en, Results, visualName)

        # call Processing algorithm to set the final style
        processing.run("qgis:setstyleforrasterlayer",{
        'INPUT': w_tot_en.source(),
        'STYLE':sty
        }, context=context, feedback=feedback)

        # add the layers created from the algorithm without explicitly declared an outputs
        context.addLayerToLoadOnCompletion(w_tot_en.source(), context.LayerDetails(
            name='w_tot_en',
            project=context.project()
        ))

        end_time = datetime.now()
        total_time = end_time - start_time
        feedback.setProgressText('Total duration of the analysis (hh:mm:ss) {}'.format(total_time))

        feedback.pushConsoleInfo('Total duration of the analysis (hh:mm:ss) {}'.format(total_time))

        return {self.OUTPUT: Results}


    def checkParameterValues(self, parameters, context):
        '''
        overwrite the checking of default parameters of Processing and throw
        an exception if the Result folder if the temporary folder
        due to incompatibility of SAGA to store information there
        '''

        if parameters['OUTPUT'] == 'TEMPORARY_OUTPUT':
            return False, self.tr("Temporary folder not supported \nPlease choose another folder destination")
        return super(QProtoAlgorithm, self).checkParameterValues(parameters, context)
