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

"""
/***************************************************************************
 ProcessingUMEP
                                 A QGIS plugin
 UMEP for processing toolbox
 Generated by Plugin Builder: http://g-sherman.github.io/Qgis-Plugin-Builder/
                              -------------------
        begin                : 2020-04-02
        copyright            : (C) 2020 by Fredrik Lindberg
        email                : fredrikl@gvc.gu.se
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   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__ = 'Fredrik Lindberg'
__date__ = '2020-04-02'
__copyright__ = '(C) 2020 by Fredrik Lindberg'

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

__revision__ = '$Format:%H$'

from qgis.PyQt.QtCore import QCoreApplication, QDate, QTime, Qt, QVariant
from qgis.core import (QgsProcessing,
                       QgsProcessingAlgorithm,
                       QgsProcessingParameterString,
                       QgsProcessingParameterBoolean,
                       QgsProcessingParameterNumber,
                       QgsProcessingParameterFolderDestination,
                       QgsProcessingParameterRasterDestination,
                       QgsProcessingParameterFileDestination,
                       QgsProcessingParameterFile,
                       QgsProcessingException,
                       QgsProcessingParameterDefinition,
                       QgsProcessingParameterEnum,
                       QgsProcessingParameterFeatureSource,
                       QgsProcessingParameterField,
                       QgsProcessingParameterRasterLayer,
                       QgsVectorLayer)
from processing.gui.wrappers import WidgetWrapper
from qgis.PyQt.QtWidgets import QDateEdit, QTimeEdit
import numpy as np
from osgeo import gdal, osr
from osgeo.gdalconst import *
import os
from qgis.PyQt.QtGui import QIcon
import inspect
from pathlib import Path, PurePath
from ..util.misc import get_ders, saveraster
import zipfile
from osgeo.gdalconst import *
from ..util.SEBESOLWEIGCommonFiles.Solweig_v2015_metdata_noload import Solweig_2015a_metdata_noload
from ..util.SEBESOLWEIGCommonFiles import Solweig_v2015_metdata_noload as metload
from ..util.SEBESOLWEIGCommonFiles.clearnessindex_2013b import clearnessindex_2013b
from ..functions.SOLWEIGpython.Tgmaps_v1 import Tgmaps_v1
from ..functions.SOLWEIGpython import Solweig_2022a_calc_forprocessing as so
from ..functions.SOLWEIGpython import WriteMetadataSOLWEIG
from ..functions.SOLWEIGpython import PET_calculations as p
from ..functions.SOLWEIGpython import UTCI_calculations as utci
from ..functions.SOLWEIGpython.CirclePlotBar import PolarBarPlot

import matplotlib.pyplot as plt

# For "Save necessary rasters for TreePlanter tool"
from shutil import copyfile

class ProcessingSOLWEIGAlgorithm(QgsProcessingAlgorithm):
    """
    This algorithm is a processing version of SOLWEIG
    """

    # Constants used to refer to parameters and outputs. They will be
    # used when calling the algorithm from another algorithm, or when
    # calling from the QGIS console.
    
    #Spatial data
    INPUT_DSM = 'INPUT_DSM'
    INPUT_SVF = 'INPUT_SVF'
    INPUT_CDSM = 'INPUT_CDSM'
    INPUT_TDSM = 'INPUT_TDSM'
    INPUT_HEIGHT = 'INPUT_HEIGHT'
    INPUT_ASPECT = 'INPUT_ASPECT'
    TRANS_VEG = 'TRANS_VEG'
    LEAF_START = 'LEAF_START'
    LEAF_END = 'LEAF_END'
    CONIFER_TREES = 'CONIFER_TREES'
    INPUT_THEIGHT = 'INPUT_THEIGHT'
    INPUT_LC = 'INPUT_LC'
    USE_LC_BUILD = 'USE_LC_BUILD'
    INPUT_DEM = 'INPUT_DEM'
    SAVE_BUILD = 'SAVE_BUILD'
    INPUT_ANISO = 'INPUT_ANISO'

    #Enivornmental parameters
    ALBEDO_WALLS = 'ALBEDO_WALLS'
    ALBEDO_GROUND = 'ALBEDO_GROUND'
    EMIS_WALLS = 'EMIS_WALLS'
    EMIS_GROUND = 'EMIS_GROUND'

    #Tmrt parameters
    ABS_S = 'ABS_S'
    ABS_L = 'ABS_L'
    POSTURE = 'POSTURE'

    #Meteorology
    INPUT_MET = 'INPUTMET'
    ONLYGLOBAL = 'ONLYGLOBAL'
    UTC = 'UTC'

    #PET parameters
    AGE = 'AGE'
    ACTIVITY = 'ACTIVITY'
    CLO = 'CLO'
    WEIGHT = 'WEIGHT'
    HEIGHT = 'HEIGHT'
    SEX = 'SEX'
    SENSOR_HEIGHT = 'SENSOR_HEIGHT'

    #Optional settings
    # POI = 'POI'
    POI_FILE = 'POI_FILE'
    POI_FIELD = 'POI_FIELD'
    CYL = 'CYL'

    #Output
    OUTPUT_DIR = 'OUTPUT_DIR'
    OUTPUT_TMRT = 'OUTPUT_TMRT'
    OUTPUT_LUP = 'OUTPUT_LUP'
    OUTPUT_KUP = 'OUTPUT_KUP'
    OUTPUT_KDOWN = 'OUTPUT_KDOWN'
    OUTPUT_LDOWN = 'OUTPUT_LDOWN'
    OUTPUT_SH = 'OUTPUT_SH'
    OUTPUT_TREEPLANTER = 'OUTPUT_TREEPLANTER'


    def initAlgorithm(self, config):
        #spatial
        self.addParameter(QgsProcessingParameterRasterLayer(self.INPUT_DSM,
            self.tr('Building and ground Digital Surface Model (DSM)'), None, optional=False))
        self.addParameter(QgsProcessingParameterFile(self.INPUT_SVF,
            self.tr('Sky View Factor grids (.zip)'), extension='zip'))
        self.addParameter(QgsProcessingParameterRasterLayer(self.INPUT_HEIGHT,
            self.tr('Wall height raster'), '', optional=False))
        self.addParameter(QgsProcessingParameterRasterLayer(self.INPUT_ASPECT,
            self.tr('Wall aspect raster'), '', optional=False))
        self.addParameter(QgsProcessingParameterRasterLayer(self.INPUT_CDSM,
            self.tr('Vegetation Canopy DSM'), '', optional=True))
        self.addParameter(QgsProcessingParameterNumber(self.TRANS_VEG,
            self.tr('Transmissivity of light through vegetation (%):'),
            QgsProcessingParameterNumber.Integer,
            QVariant(3), True, minValue=0, maxValue=100))
        
        self.addParameter(QgsProcessingParameterNumber(self.LEAF_START,
            self.tr('First day of year with leaves on trees (if deciduous)'), QgsProcessingParameterNumber.Integer,
            QVariant(97), False, minValue=0, maxValue=366))

        self.addParameter(QgsProcessingParameterNumber(self.LEAF_END,
            self.tr('Last day of year with leaves on trees (if deciduous)'), QgsProcessingParameterNumber.Integer,
            QVariant(300), False, minValue=0, maxValue=366))

        self.addParameter(QgsProcessingParameterBoolean(self.CONIFER_TREES,
            self.tr("Coniferous trees (deciduous default)"), defaultValue=False))
        
        self.addParameter(QgsProcessingParameterRasterLayer(self.INPUT_TDSM,
            self.tr('Vegetation Trunk-zone DSM'), '', optional=True))
        self.addParameter(QgsProcessingParameterNumber(self.INPUT_THEIGHT,
            self.tr("Trunk zone height (percent of Canopy Height). Used if no Vegetation Trunk-zone DSM is loaded"),
            QgsProcessingParameterNumber.Double,
            QVariant(25.0), optional=True, minValue=0.1, maxValue=99.9))
        self.addParameter(QgsProcessingParameterRasterLayer(self.INPUT_LC,
            self.tr('UMEP land cover grid'), '', optional=True))
        self.addParameter(QgsProcessingParameterBoolean(self.USE_LC_BUILD,
            self.tr("Use land cover grid to derive building grid"), defaultValue=False, optional=True))
        self.addParameter(QgsProcessingParameterRasterLayer(self.INPUT_DEM,
            self.tr('Digital Elevation Model (DEM)'), '', optional=True))
        self.addParameter(QgsProcessingParameterBoolean(self.SAVE_BUILD,
            self.tr("Save generated building grid"), defaultValue=False, optional=True))
        self.addParameter(QgsProcessingParameterFile(self.INPUT_ANISO,
            self.tr('Shadow maps used for anisotropic model for sky diffuse and longwave radiation (.npz)'), extension='npz', optional=True))

        #Environmental parameters
        self.addParameter(QgsProcessingParameterNumber(self.ALBEDO_WALLS,
            self.tr('Albedo (walls)'), QgsProcessingParameterNumber.Double,
            QVariant(0.20), False, minValue=0, maxValue=1))
        self.addParameter(QgsProcessingParameterNumber(self.ALBEDO_GROUND,
            self.tr('Albedo (ground)'), QgsProcessingParameterNumber.Double,
            QVariant(0.15), False, minValue=0, maxValue=1))        
        self.addParameter(QgsProcessingParameterNumber(self.EMIS_WALLS,
            self.tr('Emissivity (walls)'), QgsProcessingParameterNumber.Double,
            QVariant(0.90), False, minValue=0, maxValue=1)) 
        self.addParameter(QgsProcessingParameterNumber(self.EMIS_GROUND,
            self.tr('Emissivity (ground)'), QgsProcessingParameterNumber.Double,
            QVariant(0.95), False, minValue=0, maxValue=1)) 

        #Tmrt parameters
        self.addParameter(QgsProcessingParameterNumber(self.ABS_S,
            self.tr('Absorption of shortwave radiation of human body'), QgsProcessingParameterNumber.Double,
            QVariant(0.70), False, minValue=0, maxValue=1))
        self.addParameter(QgsProcessingParameterNumber(self.ABS_L,
            self.tr('Absorption of longwave radiation of human body'), QgsProcessingParameterNumber.Double,
            QVariant(0.95), False, minValue=0, maxValue=1))
        self.addParameter(QgsProcessingParameterEnum(
            self.POSTURE, self.tr('Posture of human body'), ['Standing', 'Sitting'], defaultValue=0))
        self.addParameter(QgsProcessingParameterBoolean(self.CYL,
            self.tr("Consider human as cylinder instead of box"), defaultValue=True))

        #Meteorology
        self.addParameter(QgsProcessingParameterFile(self.INPUT_MET,
            self.tr('Input meteorological file (.txt)'), extension='txt'))
        self.addParameter(QgsProcessingParameterBoolean(self.ONLYGLOBAL,
            self.tr("Estimate diffuse and direct shortwave radiation from global radiation"), defaultValue=False))
        self.addParameter(QgsProcessingParameterNumber(self.UTC,
            self.tr('Coordinated Universal Time (UTC) '),
            QgsProcessingParameterNumber.Integer,
            QVariant(0), False, minValue=-12, maxValue=12)) 
        
        
        #ADVANCED PARAMETERS
        #POIs for thermal comfort estimations
        # poi = QgsProcessingParameterBoolean(self.POI,
        #     self.tr("Include Point of Interest(s) for thermal comfort calculations (PET and UTCI)"), defaultValue=False)
        # poi.setFlags(poi.flags() | QgsProcessingParameterDefinition.FlagAdvanced)
        # self.addParameter(poi)
        poifile = QgsProcessingParameterFeatureSource(self.POI_FILE,
            self.tr('Vector point file including Point of Interest(s) for thermal comfort calculations (PET and UTCI)'), [QgsProcessing.TypeVectorPoint], optional=True)
        poifile.setFlags(poifile.flags() | QgsProcessingParameterDefinition.FlagAdvanced)
        self.addParameter(poifile)
        poi_field = QgsProcessingParameterField(self.POI_FIELD,
            self.tr('ID field'),'', self.POI_FILE, QgsProcessingParameterField.Numeric, optional=True)
        poi_field.setFlags(poi_field.flags() | QgsProcessingParameterDefinition.FlagAdvanced)
        self.addParameter(poi_field)

        #PET parameters
        age = QgsProcessingParameterNumber(self.AGE, self.tr('Age (yy)'),
                QgsProcessingParameterNumber.Integer,
                QVariant(35), optional=True, minValue=0, maxValue=120)
        age.setFlags(age.flags() | QgsProcessingParameterDefinition.FlagAdvanced)
        self.addParameter(age)
        act = QgsProcessingParameterNumber(self.ACTIVITY, self.tr('Activity (W)'),
                QgsProcessingParameterNumber.Double,
                QVariant(80), optional=True, minValue=0, maxValue=1000)
        act.setFlags(act.flags() | QgsProcessingParameterDefinition.FlagAdvanced)
        self.addParameter(act)
        clo = QgsProcessingParameterNumber(self.CLO, self.tr('Clothing (clo)'),
                QgsProcessingParameterNumber.Double,
                QVariant(0.9), optional=True, minValue=0, maxValue=10)
        clo.setFlags(clo.flags() | QgsProcessingParameterDefinition.FlagAdvanced)
        self.addParameter(clo)
        wei = QgsProcessingParameterNumber(self.WEIGHT, self.tr('Weight (kg)'),
                QgsProcessingParameterNumber.Integer,
                QVariant(75), optional=True, minValue=0, maxValue=500) 
        wei.setFlags(wei.flags() | QgsProcessingParameterDefinition.FlagAdvanced)
        self.addParameter(wei)
        hei = QgsProcessingParameterNumber(self.HEIGHT, self.tr('Height (cm)'),
                QgsProcessingParameterNumber.Integer,
                QVariant(180), optional=True, minValue=0, maxValue=250) 
        hei.setFlags(hei.flags() | QgsProcessingParameterDefinition.FlagAdvanced)
        self.addParameter(hei)
        sex = QgsProcessingParameterEnum(
            self.SEX, self.tr('Sex'), ['Male', 'Female'], optional=True, defaultValue=0)
        sex.setFlags(sex.flags() | QgsProcessingParameterDefinition.FlagAdvanced)
        self.addParameter(sex)
        shei = QgsProcessingParameterNumber(self.SENSOR_HEIGHT, self.tr('Height of wind sensor (m agl)'),
                QgsProcessingParameterNumber.Double,
                QVariant(10), optional=True, minValue=0, maxValue=250) 
        shei.setFlags(shei.flags() | QgsProcessingParameterDefinition.FlagAdvanced)
        self.addParameter(shei)

        #OUTPUT
        self.addParameter(QgsProcessingParameterBoolean(self.OUTPUT_TMRT,
            self.tr("Save Mean Radiant Temperature raster(s)"), defaultValue=True))
        self.addParameter(QgsProcessingParameterBoolean(self.OUTPUT_KDOWN,
            self.tr("Save Incoming shortwave radiation raster(s)"), defaultValue=False))
        self.addParameter(QgsProcessingParameterBoolean(self.OUTPUT_KUP,
            self.tr("Save Outgoing shortwave radiation raster(s)"), defaultValue=False))
        self.addParameter(QgsProcessingParameterBoolean(self.OUTPUT_LDOWN,
            self.tr("Save Incoming longwave radiation raster(s)"), defaultValue=False))
        self.addParameter(QgsProcessingParameterBoolean(self.OUTPUT_LUP,
            self.tr("Save Outgoing longwave radiation raster(s)"), defaultValue=False))
        self.addParameter(QgsProcessingParameterBoolean(self.OUTPUT_SH,
            self.tr("Save shadow raster(s)"), defaultValue=False))
        self.addParameter(QgsProcessingParameterBoolean(self.OUTPUT_TREEPLANTER,
            self.tr("Save necessary raster(s) for the TreePlanter and Spatial TC tools"), defaultValue=False))
        self.addParameter(QgsProcessingParameterFolderDestination(self.OUTPUT_DIR,
                                                     'Output folder'))

        self.plugin_dir = os.path.dirname(__file__)
        self.temp_dir = os.path.dirname(self.plugin_dir) + '/temp'

    def processAlgorithm(self, parameters, context, feedback):
        np.seterr(divide='ignore', invalid='ignore')
 
        # InputParameters
        dsmlayer = self.parameterAsRasterLayer(parameters, self.INPUT_DSM, context)
        transVeg = self.parameterAsDouble(parameters, self.TRANS_VEG, context) / 100.
        firstdayleaf = self.parameterAsInt(parameters, self.LEAF_START, context)
        lastdayleaf = self.parameterAsInt(parameters, self.LEAF_END, context)
        conifer_bool = self.parameterAsBool(parameters, self.CONIFER_TREES, context)
        vegdsm = self.parameterAsRasterLayer(parameters, self.INPUT_CDSM, context)
        vegdsm2 = self.parameterAsRasterLayer(parameters, self.INPUT_TDSM, context)
        lcgrid = self.parameterAsRasterLayer(parameters, self.INPUT_LC, context)
        useLcBuild = self.parameterAsBool(parameters, self.USE_LC_BUILD, context)
        dem = None
        inputSVF = self.parameterAsString(parameters, self.INPUT_SVF, context)
        whlayer = self.parameterAsRasterLayer(parameters, self.INPUT_HEIGHT, context)
        walayer = self.parameterAsRasterLayer(parameters, self.INPUT_ASPECT, context)
        trunkr = self.parameterAsDouble(parameters, self.INPUT_THEIGHT, context)
        onlyglobal = self.parameterAsBool(parameters, self.ONLYGLOBAL, context)
        utc = self.parameterAsDouble(parameters, self.UTC, context)
        inputMet = self.parameterAsString(parameters, self.INPUT_MET, context)
        # usePOI = self.parameterAsBool(parameters, self.POI, context)
        poilyr = self.parameterAsVectorLayer(parameters, self.POI_FILE, context)
        poi_field = None
        mbody = None
        ht = None
        clo = None
        age = None
        activity = None
        sex = None
        sensorheight = None
        saveBuild = self.parameterAsBool(parameters, self.SAVE_BUILD, context)
        demforbuild = 0
        folderPathPerez = self.parameterAsString(parameters, self.INPUT_ANISO, context)
        poisxy = None
        poiname = None

        # Other parameters #
        absK = self.parameterAsDouble(parameters, self.ABS_S, context)
        absL = self.parameterAsDouble(parameters, self.ABS_L, context)
        pos = self.parameterAsInt(parameters, self.POSTURE, context)
        
        if self.parameterAsBool(parameters, self.CYL, context):
            cyl = 1
        else:
            cyl = 0

        if pos == 0:
            Fside = 0.22
            Fup = 0.06
            height = 1.1
            Fcyl = 0.28
        else:
            Fside = 0.166666
            Fup = 0.166666
            height = 0.75
            Fcyl = 0.2

        albedo_b = self.parameterAsDouble(parameters, self.ALBEDO_WALLS, context)
        albedo_g = self.parameterAsDouble(parameters, self.ALBEDO_GROUND, context)
        ewall = self.parameterAsDouble(parameters, self.EMIS_WALLS, context)
        eground = self.parameterAsDouble(parameters, self.EMIS_GROUND, context)
        elvis = 0 # option removed 20200907 in processing UMEP

        outputDir = self.parameterAsString(parameters, self.OUTPUT_DIR, context)
        outputTmrt = self.parameterAsBool(parameters, self.OUTPUT_TMRT, context)
        outputSh = self.parameterAsBool(parameters, self.OUTPUT_SH, context)
        outputKup = self.parameterAsBool(parameters, self.OUTPUT_KUP, context)
        outputKdown = self.parameterAsBool(parameters, self.OUTPUT_KDOWN, context)
        outputLup = self.parameterAsBool(parameters, self.OUTPUT_LUP, context)
        outputLdown = self.parameterAsBool(parameters, self.OUTPUT_LDOWN, context)
        outputTreeplanter = self.parameterAsBool(parameters, self.OUTPUT_TREEPLANTER, context)
        outputKdiff = False
        #outputSstr = False

        # If "Save necessary rasters for TreePlanter tool" is ticked, save the following raster for TreePlanter or Spatial TC
        if outputTreeplanter:
            outputTmrt = True
            outputKup = True
            outputKdown = True
            outputLup = True
            outputLdown = True
            outputSh = True
            saveBuild = True
            outputKdiff = True
            #outputSstr = True

        if parameters['OUTPUT_DIR'] == 'TEMPORARY_OUTPUT':
            if not (os.path.isdir(outputDir)):
                os.mkdir(outputDir)

        # Code from old plugin
        provider = dsmlayer.dataProvider()
        filepath_dsm = str(provider.dataSourceUri())
        gdal_dsm = gdal.Open(filepath_dsm)
        dsm = gdal_dsm.ReadAsArray().astype(float)
        sizex = dsm.shape[0]
        sizey = dsm.shape[1]
        rows = dsm.shape[0]
        cols = dsm.shape[1]

        # response to issue #85
        nd = gdal_dsm.GetRasterBand(1).GetNoDataValue()
        dsm[dsm == nd] = 0.
        # dsmcopy = np.copy(dsm)
        if dsm.min() < 0:
            dsmraise = np.abs(dsm.min())
            dsm = dsm + dsmraise
            feedback.setProgressText('Digital Surface Model (DSM) included negative values. DSM raised with ' + str(dsmraise) + 'm.')
        else:
            dsmraise = 0

        # Get latlon from grid coordinate system
        old_cs = osr.SpatialReference()
        dsm_ref = dsmlayer.crs().toWkt()
        old_cs.ImportFromWkt(dsm_ref)

        wgs84_wkt = """
        GEOGCS["WGS 84",
            DATUM["WGS_1984",
                SPHEROID["WGS 84",6378137,298.257223563,
                    AUTHORITY["EPSG","7030"]],
                AUTHORITY["EPSG","6326"]],
            PRIMEM["Greenwich",0,
                AUTHORITY["EPSG","8901"]],
            UNIT["degree",0.01745329251994328,
                AUTHORITY["EPSG","9122"]],
            AUTHORITY["EPSG","4326"]]"""

        new_cs = osr.SpatialReference()
        new_cs.ImportFromWkt(wgs84_wkt)

        transform = osr.CoordinateTransformation(old_cs, new_cs)
        widthx = gdal_dsm.RasterXSize
        heightx = gdal_dsm.RasterYSize
        geotransform = gdal_dsm.GetGeoTransform()
        minx = geotransform[0]
        miny = geotransform[3] + widthx * geotransform[4] + heightx * geotransform[5]
        lonlat = transform.TransformPoint(minx, miny)
        gdalver = float(gdal.__version__[0])
        if gdalver == 3.:
            lon = lonlat[1] #changed to gdal 3
            lat = lonlat[0] #changed to gdal 3
        else:
            lon = lonlat[0] #changed to gdal 2
            lat = lonlat[1] #changed to gdal 2
        scale = 1 / geotransform[1]

        alt = np.median(dsm)
        if alt < 0:
            alt = 3
        feedback.setProgressText('Longitude derived from DSM: ' + str(lon))
        feedback.setProgressText('Latitude derived from DSM: ' + str(lat))

        trunkfile = 0
        trunkratio = 0
        # psi = transVeg / 100.0

        # if useVegdem:
        if vegdsm is not None:
            usevegdem = 1
            feedback.setProgressText('Vegetation scheme activated')

            # load raster
            gdal.AllRegister()
            provider = vegdsm.dataProvider()
            filePathOld = str(provider.dataSourceUri())
            dataSet = gdal.Open(filePathOld)
            vegdsm = dataSet.ReadAsArray().astype(float)
            filePath_cdsm = filePathOld
            vegsizex = vegdsm.shape[0]
            vegsizey = vegdsm.shape[1]

            if not (vegsizex == sizex) & (vegsizey == sizey):
                raise QgsProcessingException("Error in Vegetation Canopy DSM: All rasters must be of same extent and resolution")

            if vegdsm2 is not None:
                gdal.AllRegister()
                provider = vegdsm2.dataProvider()
                filePathOld = str(provider.dataSourceUri())
                filePath_tdsm = filePathOld
                dataSet = gdal.Open(filePathOld)
                vegdsm2 = dataSet.ReadAsArray().astype(float)
            else:
                trunkratio = trunkr / 100.0
                vegdsm2 = vegdsm * trunkratio
                filePath_tdsm = None

            vegsizex = vegdsm2.shape[0]
            vegsizey = vegdsm2.shape[1]

            if not (vegsizex == sizex) & (vegsizey == sizey):  # &
                raise QgsProcessingException("Error in Trunk Zone DSM: All rasters must be of same extent and resolution")
        else:
            vegdsm = np.zeros([rows, cols])
            vegdsm2 = np.zeros([rows, cols])
            usevegdem = 0
            filePath_cdsm = None
            filePath_tdsm = None

        # Land cover
        if lcgrid is not None:
            landcover = 1
            feedback.setProgressText('Land cover scheme activated')

            # load raster
            gdal.AllRegister()
            provider = lcgrid.dataProvider()
            filePath_lc = str(provider.dataSourceUri())
            dataSet = gdal.Open(filePath_lc)
            lcgrid = dataSet.ReadAsArray().astype(float)

            lcsizex = lcgrid.shape[0]
            lcsizey = lcgrid.shape[1]

            if not (lcsizex == sizex) & (lcsizey == sizey):
                raise QgsProcessingException("Error in land cover grid: All grids must be of same extent and resolution")

            baddataConifer = (lcgrid == 3)
            baddataDecid = (lcgrid == 4)
            if baddataConifer.any():
                raise QgsProcessingException("Error in land cover grid: Land cover grid includes Confier land cover class. Ground cover information (underneath canopy) is required.")
            if baddataDecid.any():
                raise QgsProcessingException("Error in land cover grid: Land cover grid includes Decidiuous land cover class. Ground cover information (underneath canopy) is required.")
            if np.isnan(lcgrid).any():
                raise QgsProcessingException("Error in land cover grid: Land cover grid includes NaN values. Use the QGIS Fill NoData cells tool to remove NaN values.")
        else:
            filePath_lc = None
            landcover = 0

        # DEM #
        if not useLcBuild:
            demforbuild = 1
            dem = self.parameterAsRasterLayer(parameters, self.INPUT_DEM, context)

            if dem is None:
                raise QgsProcessingException("Error: No valid DEM selected")

            # load raster
            gdal.AllRegister()
            provider = dem.dataProvider()
            filePathOld = str(provider.dataSourceUri())
            dataSet = gdal.Open(filePathOld)
            dem = dataSet.ReadAsArray().astype(float)

            demsizex = dem.shape[0]
            demsizey = dem.shape[1]

            if not (demsizex == sizex) & (demsizey == sizey):
                raise QgsProcessingException( "Error in DEM: All grids must be of same extent and resolution")

            # response to issue and #230
            nd = dataSet.GetRasterBand(1).GetNoDataValue()
            dem[dem == nd] = 0.
            if dem.min() < 0:
                demraise = np.abs(dem.min())
                dem = dem + demraise
                feedback.setProgressText('Digital Evevation Model (DEM) included negative values. DEM raised with ' + str(demraise) + 'm.')
            else:
                demraise = 0

            alt = np.median(dem)
            if alt > 0:
                alt = 3.

            if (dsmraise != demraise) and (dsmraise - demraise > 0.5):
                feedback.setProgressText('WARNiNG! DEM and DSM was raised unequally (difference > 0.5 m). Check your input data!')

        #SVFs
        zip = zipfile.ZipFile(inputSVF, 'r')
        zip.extractall(self.temp_dir)
        zip.close()

        try:
            dataSet = gdal.Open(self.temp_dir + "/svf.tif")
            svf = dataSet.ReadAsArray().astype(float)
            dataSet = gdal.Open(self.temp_dir + "/svfN.tif")
            svfN = dataSet.ReadAsArray().astype(float)
            dataSet = gdal.Open(self.temp_dir + "/svfS.tif")
            svfS = dataSet.ReadAsArray().astype(float)
            dataSet = gdal.Open(self.temp_dir + "/svfE.tif")
            svfE = dataSet.ReadAsArray().astype(float)
            dataSet = gdal.Open(self.temp_dir + "/svfW.tif")
            svfW = dataSet.ReadAsArray().astype(float)

            if usevegdem == 1:
                dataSet = gdal.Open(self.temp_dir + "/svfveg.tif")
                svfveg = dataSet.ReadAsArray().astype(float)
                dataSet = gdal.Open(self.temp_dir + "/svfNveg.tif")
                svfNveg = dataSet.ReadAsArray().astype(float)
                dataSet = gdal.Open(self.temp_dir + "/svfSveg.tif")
                svfSveg = dataSet.ReadAsArray().astype(float)
                dataSet = gdal.Open(self.temp_dir + "/svfEveg.tif")
                svfEveg = dataSet.ReadAsArray().astype(float)
                dataSet = gdal.Open(self.temp_dir + "/svfWveg.tif")
                svfWveg = dataSet.ReadAsArray().astype(float)

                dataSet = gdal.Open(self.temp_dir + "/svfaveg.tif")
                svfaveg = dataSet.ReadAsArray().astype(float)
                dataSet = gdal.Open(self.temp_dir + "/svfNaveg.tif")
                svfNaveg = dataSet.ReadAsArray().astype(float)
                dataSet = gdal.Open(self.temp_dir + "/svfSaveg.tif")
                svfSaveg = dataSet.ReadAsArray().astype(float)
                dataSet = gdal.Open(self.temp_dir + "/svfEaveg.tif")
                svfEaveg = dataSet.ReadAsArray().astype(float)
                dataSet = gdal.Open(self.temp_dir + "/svfWaveg.tif")
                svfWaveg = dataSet.ReadAsArray().astype(float)
            else:
                svfveg = np.ones((rows, cols))
                svfNveg = np.ones((rows, cols))
                svfSveg = np.ones((rows, cols))
                svfEveg = np.ones((rows, cols))
                svfWveg = np.ones((rows, cols))
                svfaveg = np.ones((rows, cols))
                svfNaveg = np.ones((rows, cols))
                svfSaveg = np.ones((rows, cols))
                svfEaveg = np.ones((rows, cols))
                svfWaveg = np.ones((rows, cols))
        except:
            raise QgsProcessingException("SVF import error: The zipfile including the SVFs seems corrupt. Retry calcualting the SVFs in the Pre-processor or choose another file.")

        svfsizex = svf.shape[0]
        svfsizey = svf.shape[1]

        if not (svfsizex == sizex) & (svfsizey == sizey):  # &
            raise QgsProcessingException("Error in svf rasters: All grids must be of same extent and resolution")

        tmp = svf + svfveg - 1.
        tmp[tmp < 0.] = 0.
        # %matlab crazyness around 0
        svfalfa = np.arcsin(np.exp((np.log((1. - tmp)) / 2.)))

        feedback.setProgressText('Sky View Factor rasters loaded')

        # wall height layer
        if whlayer is None:
            raise QgsProcessingException("Error: No valid wall height raster layer is selected")
        provider = whlayer.dataProvider()
        filepath_wh = str(provider.dataSourceUri())
        self.gdal_wh = gdal.Open(filepath_wh)
        wallheight = self.gdal_wh.ReadAsArray().astype(float)
        vhsizex = wallheight.shape[0]
        vhsizey = wallheight.shape[1]
        if not (vhsizex == sizex) & (vhsizey == sizey):
            raise QgsProcessingException("Error in Wall height raster: All rasters must be of same extent and resolution")

        # wall aspectlayer
        if walayer is None:
            raise QgsProcessingException("Error: No valid wall aspect raster layer is selected")
        provider = walayer.dataProvider()
        filepath_wa = str(provider.dataSourceUri())
        self.gdal_wa = gdal.Open(filepath_wa)
        wallaspect = self.gdal_wa.ReadAsArray().astype(float)
        vasizex = wallaspect.shape[0]
        vasizey = wallaspect.shape[1]
        if not (vasizex == sizex) & (vasizey == sizey):
            raise QgsProcessingException("Error in Wall aspect raster: All rasters must be of same extent and resolution")

        voxelheight = geotransform[1]  # float

        # Metdata
        headernum = 1
        delim = ' '
        Twater = []

        try:
            self.metdata = np.loadtxt(inputMet,skiprows=headernum, delimiter=delim)
            metfileexist = 1
        except:
            raise QgsProcessingException("Error: Make sure format of meteorological file is correct. You can"
                                                        "prepare your data by using 'Prepare Existing Data' in "
                                                        "the Pre-processor")

        testwhere = np.where((self.metdata[:, 14] < 0.0) | (self.metdata[:, 14] > 1300.0))
        if testwhere[0].__len__() > 0:
             raise QgsProcessingException("Error: Kdown - beyond what is expected at line: " + str(testwhere[0] + 1))

        if self.metdata.shape[1] == 24:
            feedback.setProgressText("Meteorological data successfully loaded")
        else:
            raise QgsProcessingException("Error: Wrong number of columns in meteorological data. You can "
                                                        "prepare your data by using 'Prepare Existing Data' in "
                                                        "the Pre-processor")

        feedback.setProgressText("Calculating sun positions for each time step")
        location = {'longitude': lon, 'latitude': lat, 'altitude': alt}
        YYYY, altitude, azimuth, zen, jday, leafon, dectime, altmax = \
            Solweig_2015a_metdata_noload(self.metdata,location, utc)

        # Creating vectors from meteorological input
        DOY = self.metdata[:, 1]
        hours = self.metdata[:, 2]
        minu = self.metdata[:, 3]
        Ta = self.metdata[:, 11]
        RH = self.metdata[:, 10]
        radG = self.metdata[:, 14]
        radD = self.metdata[:, 21]
        radI = self.metdata[:, 22]
        P = self.metdata[:, 12]
        Ws = self.metdata[:, 9]
        
        # Check if diffuse and direct radiation exist
        if onlyglobal == 0:
            if np.min(radD) == -999:
                raise QgsProcessingException("Diffuse radiation include NoData values",
                                        'Tick in the box "Estimate diffuse and direct shortwave..." or aqcuire '
                                        'observed values from external data sources.')
            if np.min(radI) == -999:
                raise QgsProcessingException("Direct radiation include NoData values",
                                        'Tick in the box "Estimate diffuse and direct shortwave..." or aqcuire '
                                        'observed values from external data sources.')

        # POIs check
        if poilyr is not None: # usePOI:
            #header = 'yyyy id   it imin dectime altitude azimuth kdir kdiff kglobal kdown   kup    keast ksouth ' \
            #            'kwest knorth ldown   lup    least lsouth lwest  lnorth   Ta      Tg     RH    Esky   Tmrt    ' \
            #            'I0     CI   Shadow  SVF_b  SVF_bv KsideI PET UTCI'

            header = 'yyyy id   it imin dectime altitude azimuth kdir kdiff kglobal kdown   kup    keast ksouth ' \
                        'kwest knorth ldown   lup    least lsouth lwest  lnorth   Ta      Tg     RH    Esky   Tmrt    ' \
                        'I0     CI   Shadow  SVF_b  SVF_bv KsideI PET UTCI  CI_Tg   CI_TgG  KsideD  Lside   diffDown    Kside'                            

            # poilyr = self.parameterAsVectorLayer(parameters, self.POI_FILE, context) 
            # if poilyr is None:
                # raise QgsProcessingException("No valid point layer is selected")

            poi_field = self.parameterAsFields(parameters, self.POI_FIELD, context)
            # if poi_field[0] is None:
            #     raise QgsProcessingException("An attribute field with unique values must be selected when using a POI vector file")
            vlayer = poilyr
            prov = vlayer.dataProvider()
            fields = prov.fields()
            idx = vlayer.fields().indexFromName(poi_field[0])
            numfeat = vlayer.featureCount()
            poiname = []
            poisxy = np.zeros((numfeat, 3)) - 999
            ind = 0
            for f in vlayer.getFeatures():  # looping through each POI
                y = f.geometry().centroid().asPoint().y()
                x = f.geometry().centroid().asPoint().x()

                poiname.append(f.attributes()[idx])
                poisxy[ind, 0] = ind
                poisxy[ind, 1] = np.round((x - minx) * scale)
                if miny >= 0:
                    poisxy[ind, 2] = np.round((miny + rows * (1. / scale) - y) * scale)
                else:
                    poisxy[ind, 2] = np.round((miny + rows * (1. / scale) - y) * scale)

                ind += 1

            for k in range(0, poisxy.shape[0]):
                poi_save = []  # np.zeros((1, 33))
                data_out = outputDir + '/POI_' + str(poiname[k]) + '.txt'
                np.savetxt(data_out, poi_save,  delimiter=' ', header=header, comments='')
            
            # Other PET variables
            mbody = self.parameterAsDouble(parameters, self.WEIGHT, context)
            ht = self.parameterAsDouble(parameters, self.HEIGHT, context) / 100.
            clo = self.parameterAsDouble(parameters, self.CLO, context)
            age = self.parameterAsDouble(parameters, self.AGE, context)
            activity = self.parameterAsDouble(parameters, self.WEIGHT, context)
            sex = self.parameterAsInt(parameters, self.SEX, context) + 1
            sensorheight = self.parameterAsDouble(parameters, self.SENSOR_HEIGHT, context)
            
            feedback.setProgressText("Point of interest (POI) vector data successfully loaded")

        # %Parameterisarion for Lup
        if not height:
            height = 1.1

        # %Radiative surface influence, Rule of thumb by Schmid et al. (1990).
        first = np.round(height)
        if first == 0.:
            first = 1.
        second = np.round((height * 20.))

        if usevegdem == 1:
            # Conifer or deciduous
            if conifer_bool:
                leafon = np.ones((1, DOY.shape[0]))
            else:
                leafon = np.zeros((1, DOY.shape[0]))
                if firstdayleaf > lastdayleaf:
                    leaf_bool = ((DOY > firstdayleaf) | (DOY < lastdayleaf))
                else:
                    leaf_bool = ((DOY > firstdayleaf) & (DOY < lastdayleaf))
                leafon[0, leaf_bool] = 1

            # % Vegetation transmittivity of shortwave radiation
            psi = leafon * transVeg
            psi[leafon == 0] = 0.5
            # amaxvalue
            vegmax = vegdsm.max()
            amaxvalue = dsm.max() - dsm.min()
            amaxvalue = np.maximum(amaxvalue, vegmax)

            # Elevation vegdsms if buildingDEM includes ground heights
            vegdsm = vegdsm + dsm
            vegdsm[vegdsm == dsm] = 0
            vegdsm2 = vegdsm2 + dsm
            vegdsm2[vegdsm2 == dsm] = 0

            # % Bush separation
            bush = np.logical_not((vegdsm2 * vegdsm)) * vegdsm

            svfbuveg = (svf - (1. - svfveg) * (1. - transVeg))  # % major bug fixed 20141203
        else:
            psi = leafon * 0. + 1.
            svfbuveg = svf
            bush = np.zeros([rows, cols])
            amaxvalue = 0

        # %Initialization of maps
        Knight = np.zeros((rows, cols))
        Tgmap1 = np.zeros((rows, cols))
        Tgmap1E = np.zeros((rows, cols))
        Tgmap1S = np.zeros((rows, cols))
        Tgmap1W = np.zeros((rows, cols))
        Tgmap1N = np.zeros((rows, cols))
        
        # building grid and land cover preparation
        sitein = self.plugin_dir + "/landcoverclasses_2016a.txt"
        f = open(sitein)
        lin = f.readlines()
        lc_class = np.zeros((lin.__len__() - 1, 6))
        for i in range(1, lin.__len__()):
            lines = lin[i].split()
            for j in np.arange(1, 7):
                lc_class[i - 1, j - 1] = float(lines[j])

        if demforbuild == 0:
            buildings = np.copy(lcgrid)
            buildings[buildings == 7] = 1
            buildings[buildings == 6] = 1
            buildings[buildings == 5] = 1
            buildings[buildings == 4] = 1
            buildings[buildings == 3] = 1
            buildings[buildings == 2] = 0
        else:
            buildings = dsm - dem
            buildings[buildings < 2.] = 1.
            buildings[buildings >= 2.] = 0.

        if saveBuild:
            saveraster(gdal_dsm, outputDir + '/buildings.tif', buildings)

        # Import shadow matrices (Anisotropic sky)
        if folderPathPerez:  #UseAniso
            anisotropic_sky = 1
            data = np.load(folderPathPerez)
            shmat = data['shadowmat']
            vegshmat = data['vegshadowmat']
            vbshvegshmat = data['vbshmat']
            if usevegdem == 1:
                diffsh = np.zeros((rows, cols, shmat.shape[2]))
                for i in range(0, shmat.shape[2]):
                    diffsh[:, :, i] = shmat[:, :, i] - (1 - vegshmat[:, :, i]) * (1 - transVeg) # changes in psi not implemented yet
            else:
                diffsh = shmat
                vegshmat += 1
                vbshvegshmat += 1

            # Estimate number of patches based on shadow matrices
            if shmat.shape[2] == 145:
                patch_option = 1 # patch_option = 1 # 145 patches
            elif shmat.shape[2] == 153:
                patch_option = 2 # patch_option = 2 # 153 patches
            elif shmat.shape[2] == 306:
                patch_option = 3 # patch_option = 3 # 306 patches
            elif shmat.shape[2] == 612:
                patch_option = 4 # patch_option = 4 # 612 patches

            # asvf to calculate sunlit and shaded patches
            asvf = np.arccos(np.sqrt(svf))

            anisotropic_feedback = "Sky divided into " + str(int(shmat.shape[2])) + " patches\n \
                                    Anisotropic sky for diffuse shortwave radiation (Perez et al., 1993) and longwave radiation (Martin & Berdahl, 1984)"
            feedback.setProgressText(anisotropic_feedback)
        else:
            feedback.setProgressText("Isotropic sky")
            anisotropic_sky = 0
            diffsh = None
            shmat = None
            vegshmat = None
            vbshvegshmat = None
            asvf = None
            patch_option = 0

        # % Ts parameterisation maps
        if landcover == 1.:
            if np.max(lcgrid) > 21 or np.min(lcgrid) < 1:
                raise QgsProcessingException("The land cover grid includes integer values higher (or lower) than UMEP-formatted" 
                    "land cover grid (should be integer between 1 and 7). If other LC-classes should be included they also need to be included in landcoverclasses_2016a.txt")
            if np.where(lcgrid) == 3 or np.where(lcgrid) == 4:
                raise QgsProcessingException("The land cover grid includes values (decidouos and/or conifer) not appropriate for SOLWEIG-formatted land cover grid (should not include 3 or 4).")

            [TgK, Tstart, alb_grid, emis_grid, TgK_wall, Tstart_wall, TmaxLST, TmaxLST_wall] = Tgmaps_v1(lcgrid, lc_class)
        else:
            TgK = Knight + 0.37
            Tstart = Knight - 3.41
            alb_grid = Knight + albedo_g
            emis_grid = Knight + eground
            TgK_wall = 0.37
            Tstart_wall = -3.41
            TmaxLST = 15.
            TmaxLST_wall = 15.

         # Initialisation of time related variables
        if Ta.__len__() == 1:
            timestepdec = 0
        else:
            timestepdec = dectime[1] - dectime[0]
        timeadd = 0.
        timeaddE = 0.
        timeaddS = 0.
        timeaddW = 0.
        timeaddN = 0.
        firstdaytime = 1.

        WriteMetadataSOLWEIG.writeRunInfo(outputDir, filepath_dsm, gdal_dsm, usevegdem,
                                              filePath_cdsm, trunkfile, filePath_tdsm, lat, lon, utc, landcover,
                                              filePath_lc, metfileexist, inputMet, self.metdata, self.plugin_dir,
                                              absK, absL, albedo_b, albedo_g, ewall, eground, onlyglobal, trunkratio,
                                              transVeg, rows, cols, pos, elvis, cyl, demforbuild, anisotropic_sky)

        feedback.setProgressText("Writing settings for this model run to specified output folder (Filename: RunInfoSOLWEIG_YYYY_DOY_HHMM.txt)")

        # Save svf
        if anisotropic_sky:
            if not poisxy is None:
                patch_characteristics = np.zeros((shmat.shape[2], poisxy.shape[0]))
                for idx in range(poisxy.shape[0]):
                    for idy in range(shmat.shape[2]):
                        # Calculations for patches on sky, shmat = 1 = sky is visible
                        temp_sky = ((shmat[:,:,idy] == 1) & (vegshmat[:,:,idy] == 1))
                        # Calculations for patches that are vegetation, vegshmat = 0 = shade from vegetation
                        temp_vegsh = ((vegshmat[:,:,idy] == 0) | (vbshvegshmat[:,:,idy] == 0))
                        # Calculations for patches that are buildings, shmat = 0 = shade from buildings
                        temp_vbsh = (1 - shmat[:,:,idy]) * vbshvegshmat[:,:,idy]
                        temp_sh = (temp_vbsh == 1)
                        # Sky patch
                        if temp_sky[int(poisxy[idx, 2]), int(poisxy[idx, 1])]:
                            patch_characteristics[idy,idx] = 1.8
                        # Vegetation patch
                        elif (temp_vegsh[int(poisxy[idx, 2]), int(poisxy[idx, 1])]):
                            patch_characteristics[idy,idx] = 2.5
                        # Building patch
                        elif (temp_sh[int(poisxy[idx, 2]), int(poisxy[idx, 1])]):
                            patch_characteristics[idy,idx] = 4.5

        #  If metfile starts at night
        CI = 1.

        # Main function
        feedback.setProgressText("Executing main model")
    
        tmrtplot = np.zeros((rows, cols))
        TgOut1 = np.zeros((rows, cols))

        # Initiate array for I0 values
        if np.unique(DOY).shape[0] > 1:
            unique_days = np.unique(DOY)
            first_unique_day = DOY[DOY == unique_days[0]]
            I0_array = np.zeros((first_unique_day.shape[0]))
        else:
            first_unique_day = DOY.copy()
            I0_array = np.zeros((DOY.shape[0]))

        #numformat = '%d %d %d %d %.5f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f ' \
        #            '%.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f'

        numformat = '%d %d %d %d %.5f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f ' \
                    '%.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f ' \
                        '%.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f'

        for i in np.arange(0, Ta.__len__()):
            feedback.setProgress(int(i * (100. / Ta.__len__()))) # move progressbar forward
            if feedback.isCanceled():
                feedback.setProgressText("Calculation cancelled")
                break
            # Daily water body temperature
            if landcover == 1:
                if ((dectime[i] - np.floor(dectime[i]))) == 0 or (i == 0):
                    Twater = np.mean(Ta[jday[0] == np.floor(dectime[i])])
            # Nocturnal cloudfraction from Offerle et al. 2003
            if (dectime[i] - np.floor(dectime[i])) == 0:
                daylines = np.where(np.floor(dectime) == dectime[i])
                if daylines.__len__() > 1:
                    alt = altitude[0][daylines]
                    alt2 = np.where(alt > 1)
                    rise = alt2[0][0]
                    [_, CI, _, _, _] = clearnessindex_2013b(zen[0, i + rise + 1], jday[0, i + rise + 1],
                                                            Ta[i + rise + 1],
                                                            RH[i + rise + 1] / 100., radG[i + rise + 1], location,
                                                            P[i + rise + 1])  # i+rise+1 to match matlab code. correct?
                    if (CI > 1.) or (CI == np.inf):
                        CI = 1.
                else:
                    CI = 1.

            # radI[i] = radI[i]/np.sin(altitude[0][i] * np.pi/180)

            Tmrt, Kdown, Kup, Ldown, Lup, Tg, ea, esky, I0, CI, shadow, firstdaytime, timestepdec, timeadd, \
                    Tgmap1, Tgmap1E, Tgmap1S, Tgmap1W, Tgmap1N, Keast, Ksouth, Kwest, Knorth, Least, \
                    Lsouth, Lwest, Lnorth, KsideI, TgOut1, TgOut, radIout, radDout, \
                    Lside, Lsky_patch_characteristics, CI_Tg, CI_TgG, KsideD, \
                        dRad, Kside = so.Solweig_2022a_calc(
                        i, dsm, scale, rows, cols, svf, svfN, svfW, svfE, svfS, svfveg,
                        svfNveg, svfEveg, svfSveg, svfWveg, svfaveg, svfEaveg, svfSaveg, svfWaveg, svfNaveg, \
                        vegdsm, vegdsm2, albedo_b, absK, absL, ewall, Fside, Fup, Fcyl, altitude[0][i],
                        azimuth[0][i], zen[0][i], jday[0][i], usevegdem, onlyglobal, buildings, location,
                        psi[0][i], landcover, lcgrid, dectime[i], altmax[0][i], wallaspect,
                        wallheight, cyl, elvis, Ta[i], RH[i], radG[i], radD[i], radI[i], P[i], amaxvalue,
                        bush, Twater, TgK, Tstart, alb_grid, emis_grid, TgK_wall, Tstart_wall, TmaxLST,
                        TmaxLST_wall, first, second, svfalfa, svfbuveg, firstdaytime, timeadd, timestepdec, 
                        Tgmap1, Tgmap1E, Tgmap1S, Tgmap1W, Tgmap1N, CI, TgOut1, diffsh, shmat, vegshmat, vbshvegshmat, 
                        anisotropic_sky, asvf, patch_option)

            # Tmrt, Kdown, Kup, Ldown, Lup, Tg, ea, esky, I0, CI, shadow, firstdaytime, timestepdec, timeadd, \
            #         Tgmap1, Tgmap1E, Tgmap1S, Tgmap1W, Tgmap1N, Keast, Ksouth, Kwest, Knorth, Least, \
            #         Lsouth, Lwest, Lnorth, KsideI, TgOut1, TgOut, radIout, radDout = so.Solweig_2021a_calc(
            #             i, dsm, scale, rows, cols, svf, svfN, svfW, svfE, svfS, svfveg,
            #             svfNveg, svfEveg, svfSveg, svfWveg, svfaveg, svfEaveg, svfSaveg, svfWaveg, svfNaveg,
            #             vegdsm, vegdsm2, albedo_b, absK, absL, ewall, Fside, Fup, Fcyl, altitude[0][i],
            #             azimuth[0][i], zen[0][i], jday[0][i], usevegdem, onlyglobal, buildings, location,
            #             psi[0][i], landcover, lcgrid, dectime[i], altmax[0][i], wallaspect,
            #             wallheight, cyl, elvis, Ta[i], RH[i], radG[i], radD[i], radI[i], P[i], amaxvalue,
            #             bush, Twater, TgK, Tstart, alb_grid, emis_grid, TgK_wall, Tstart_wall, TmaxLST,
            #             TmaxLST_wall, first, second, svfalfa, svfbuveg, firstdaytime, timeadd, timestepdec, 
            #             Tgmap1, Tgmap1E, Tgmap1S, Tgmap1W, Tgmap1N, CI, TgOut1, diffsh, ani)



            # Tmrt, Kdown, Kup, Ldown, Lup, Tg, ea, esky, I0, CI, shadow, firstdaytime, timestepdec, timeadd, \
            # Tgmap1, timeaddE, Tgmap1E, timeaddS, Tgmap1S, timeaddW, Tgmap1W, timeaddN, Tgmap1N, \
            # Keast, Ksouth, Kwest, Knorth, Least, Lsouth, Lwest, Lnorth, KsideI, TgOut1, TgOut, radIout, radDout \
            #     = so.Solweig_2019a_calc(i, dsm, scale, rows, cols, svf, svfN, svfW, svfE, svfS, svfveg,
            #         svfNveg, svfEveg, svfSveg, svfWveg, svfaveg, svfEaveg, svfSaveg, svfWaveg, svfNaveg,
            #         vegdsm, vegdsm2, albedo_b, absK, absL, ewall, Fside, Fup, Fcyl, altitude[0][i],
            #         azimuth[0][i], zen[0][i], jday[0][i], usevegdem, onlyglobal, buildings, location,
            #         psi[0][i], landcover, lcgrid, dectime[i], altmax[0][i], waspect,
            #         wheight, cyl, elvis, Ta[i], RH[i], radG[i], radD[i], radI[i], P[i], amaxvalue,
            #         bush, Twater, TgK, Tstart, alb_grid, emis_grid, TgK_wall, Tstart_wall, TmaxLST,
            #         TmaxLST_wall, first, second, svfalfa, svfbuveg, firstdaytime, timeadd, timeaddE, timeaddS,
            #         timeaddW, timeaddN, timestepdec, Tgmap1, Tgmap1E, Tgmap1S, Tgmap1W, Tgmap1N, CI, TgOut1, diffsh, ani)

            # Save I0 for I0 vs. Kdown output plot to check if UTC is off
            if i < first_unique_day.shape[0]:
                I0_array[i] = I0

            tmrtplot = tmrtplot + Tmrt

            if altitude[0][i] > 0:
                w = 'D'
            else:
                w = 'N'

            # # Write to POIs
            # if not poisxy is None:
            #     for k in range(0, poisxy.shape[0]):
            #         poi_save = np.zeros((1, 35))
            #         poi_save[0, 0] = YYYY[0][i]
            #         poi_save[0, 1] = jday[0][i]
            #         poi_save[0, 2] = hours[i]
            #         poi_save[0, 3] = minu[i]
            #         poi_save[0, 4] = dectime[i]
            #         poi_save[0, 5] = altitude[0][i]
            #         poi_save[0, 6] = azimuth[0][i]
            #         poi_save[0, 7] = radIout
            #         poi_save[0, 8] = radDout
            #         poi_save[0, 9] = radG[i]
            #         poi_save[0, 10] = Kdown[int(poisxy[k, 2]), int(poisxy[k, 1])]
            #         poi_save[0, 11] = Kup[int(poisxy[k, 2]), int(poisxy[k, 1])]
            #         poi_save[0, 12] = Keast[int(poisxy[k, 2]), int(poisxy[k, 1])]
            #         poi_save[0, 13] = Ksouth[int(poisxy[k, 2]), int(poisxy[k, 1])]
            #         poi_save[0, 14] = Kwest[int(poisxy[k, 2]), int(poisxy[k, 1])]
            #         poi_save[0, 15] = Knorth[int(poisxy[k, 2]), int(poisxy[k, 1])]
            #         poi_save[0, 16] = Ldown[int(poisxy[k, 2]), int(poisxy[k, 1])]
            #         poi_save[0, 17] = Lup[int(poisxy[k, 2]), int(poisxy[k, 1])]
            #         poi_save[0, 18] = Least[int(poisxy[k, 2]), int(poisxy[k, 1])]
            #         poi_save[0, 19] = Lsouth[int(poisxy[k, 2]), int(poisxy[k, 1])]
            #         poi_save[0, 20] = Lwest[int(poisxy[k, 2]), int(poisxy[k, 1])]
            #         poi_save[0, 21] = Lnorth[int(poisxy[k, 2]), int(poisxy[k, 1])]
            #         poi_save[0, 22] = Ta[i]
            #         poi_save[0, 23] = TgOut[int(poisxy[k, 2]), int(poisxy[k, 1])]
            #         poi_save[0, 24] = RH[i]
            #         poi_save[0, 25] = esky
            #         poi_save[0, 26] = Tmrt[int(poisxy[k, 2]), int(poisxy[k, 1])]
            #         poi_save[0, 27] = I0
            #         poi_save[0, 28] = CI
            #         poi_save[0, 29] = shadow[int(poisxy[k, 2]), int(poisxy[k, 1])]
            #         poi_save[0, 30] = svf[int(poisxy[k, 2]), int(poisxy[k, 1])]
            #         poi_save[0, 31] = svfbuveg[int(poisxy[k, 2]), int(poisxy[k, 1])]
            #         poi_save[0, 32] = KsideI[int(poisxy[k, 2]), int(poisxy[k, 1])]
            #         # Recalculating wind speed based on powerlaw
            #         WsPET = (1.1 / sensorheight) ** 0.2 * Ws[i]
            #         WsUTCI = (10. / sensorheight) ** 0.2 * Ws[i]
            #         resultPET = p._PET(Ta[i], RH[i], Tmrt[int(poisxy[k, 2]), int(poisxy[k, 1])], WsPET,
            #                             mbody, age, ht, activity, clo, sex)
            #         poi_save[0, 33] = resultPET
            #         resultUTCI = utci.utci_calculator(Ta[i], RH[i], Tmrt[int(poisxy[k, 2]), int(poisxy[k, 1])],
            #                                             WsUTCI)
            #         poi_save[0, 34] = resultUTCI
            #         data_out = outputDir + '/POI_' + str(poiname[k]) + '.txt'
            #         # f_handle = file(data_out, 'a')
            #         f_handle = open(data_out, 'ab')
            #         np.savetxt(f_handle, poi_save, fmt=numformat)
            #         f_handle.close()

            # Write to POIs
            if not poisxy is None:
                for k in range(0, poisxy.shape[0]):
                    poi_save = np.zeros((1, 41))
                    poi_save[0, 0] = YYYY[0][i]
                    poi_save[0, 1] = jday[0][i]
                    poi_save[0, 2] = hours[i]
                    poi_save[0, 3] = minu[i]
                    poi_save[0, 4] = dectime[i]
                    poi_save[0, 5] = altitude[0][i]
                    poi_save[0, 6] = azimuth[0][i]
                    poi_save[0, 7] = radIout
                    poi_save[0, 8] = radDout
                    poi_save[0, 9] = radG[i]
                    poi_save[0, 10] = Kdown[int(poisxy[k, 2]), int(poisxy[k, 1])]
                    poi_save[0, 11] = Kup[int(poisxy[k, 2]), int(poisxy[k, 1])]
                    poi_save[0, 12] = Keast[int(poisxy[k, 2]), int(poisxy[k, 1])]
                    poi_save[0, 13] = Ksouth[int(poisxy[k, 2]), int(poisxy[k, 1])]
                    poi_save[0, 14] = Kwest[int(poisxy[k, 2]), int(poisxy[k, 1])]
                    poi_save[0, 15] = Knorth[int(poisxy[k, 2]), int(poisxy[k, 1])]
                    poi_save[0, 16] = Ldown[int(poisxy[k, 2]), int(poisxy[k, 1])]
                    poi_save[0, 17] = Lup[int(poisxy[k, 2]), int(poisxy[k, 1])]
                    poi_save[0, 18] = Least[int(poisxy[k, 2]), int(poisxy[k, 1])]
                    poi_save[0, 19] = Lsouth[int(poisxy[k, 2]), int(poisxy[k, 1])]
                    poi_save[0, 20] = Lwest[int(poisxy[k, 2]), int(poisxy[k, 1])]
                    poi_save[0, 21] = Lnorth[int(poisxy[k, 2]), int(poisxy[k, 1])]
                    poi_save[0, 22] = Ta[i]
                    poi_save[0, 23] = TgOut[int(poisxy[k, 2]), int(poisxy[k, 1])]
                    poi_save[0, 24] = RH[i]
                    poi_save[0, 25] = esky
                    poi_save[0, 26] = Tmrt[int(poisxy[k, 2]), int(poisxy[k, 1])]
                    poi_save[0, 27] = I0
                    poi_save[0, 28] = CI
                    poi_save[0, 29] = shadow[int(poisxy[k, 2]), int(poisxy[k, 1])]
                    poi_save[0, 30] = svf[int(poisxy[k, 2]), int(poisxy[k, 1])]
                    poi_save[0, 31] = svfbuveg[int(poisxy[k, 2]), int(poisxy[k, 1])]
                    poi_save[0, 32] = KsideI[int(poisxy[k, 2]), int(poisxy[k, 1])]
                    # Recalculating wind speed based on powerlaw
                    WsPET = (1.1 / sensorheight) ** 0.2 * Ws[i]
                    WsUTCI = (10. / sensorheight) ** 0.2 * Ws[i]
                    resultPET = p._PET(Ta[i], RH[i], Tmrt[int(poisxy[k, 2]), int(poisxy[k, 1])], WsPET,
                                        mbody, age, ht, activity, clo, sex)
                    poi_save[0, 33] = resultPET
                    resultUTCI = utci.utci_calculator(Ta[i], RH[i], Tmrt[int(poisxy[k, 2]), int(poisxy[k, 1])],
                                                        WsUTCI)
                    poi_save[0, 34] = resultUTCI
                    poi_save[0, 35] = CI_Tg
                    poi_save[0, 36] = CI_TgG
                    poi_save[0, 37] = KsideD[int(poisxy[k, 2]), int(poisxy[k, 1])]
                    poi_save[0, 38] = Lside[int(poisxy[k, 2]), int(poisxy[k, 1])]
                    poi_save[0, 39] = dRad[int(poisxy[k, 2]), int(poisxy[k, 1])]
                    poi_save[0, 40] = Kside[int(poisxy[k, 2]), int(poisxy[k, 1])]
                    data_out = outputDir + '/POI_' + str(poiname[k]) + '.txt'
                    # f_handle = file(data_out, 'a')
                    f_handle = open(data_out, 'ab')
                    np.savetxt(f_handle, poi_save, fmt=numformat)
                    f_handle.close()

            if hours[i] < 10:
                XH = '0'
            else:
                XH = ''
            if minu[i] < 10:
                XM = '0'
            else:
                XM = ''

            if outputTmrt:
                saveraster(gdal_dsm, outputDir + '/Tmrt_' + str(int(YYYY[0, i])) + '_' + str(int(DOY[i]))
                                + '_' + XH + str(int(hours[i])) + XM + str(int(minu[i])) + w + '.tif', Tmrt)
            if outputKup:
                saveraster(gdal_dsm, outputDir + '/Kup_' + str(int(YYYY[0, i])) + '_' + str(int(DOY[i]))
                                + '_' + XH + str(int(hours[i])) + XM + str(int(minu[i])) + w + '.tif', Kup)
            if outputKdown:
                saveraster(gdal_dsm, outputDir + '/Kdown_' + str(int(YYYY[0, i])) + '_' + str(int(DOY[i]))
                                + '_' + XH + str(int(hours[i])) + XM + str(int(minu[i])) + w + '.tif', Kdown)
            if outputLup:
                saveraster(gdal_dsm, outputDir + '/Lup_' + str(int(YYYY[0, i])) + '_' + str(int(DOY[i]))
                                + '_' + XH + str(int(hours[i])) + XM + str(int(minu[i])) + w + '.tif', Lup)
            if outputLdown:
                saveraster(gdal_dsm, outputDir + '/Ldown_' + str(int(YYYY[0, i])) + '_' + str(int(DOY[i]))
                                + '_' + XH + str(int(hours[i])) + XM + str(int(minu[i])) + w + '.tif', Ldown)
            if outputSh:
                saveraster(gdal_dsm, outputDir + '/Shadow_' + str(int(YYYY[0, i])) + '_' + str(int(DOY[i]))
                                + '_' + XH + str(int(hours[i])) + XM + str(int(minu[i])) + w + '.tif', shadow)
                
            if outputKdiff:
                saveraster(gdal_dsm, outputDir + '/Kdiff_' + str(int(YYYY[0, i])) + '_' + str(int(DOY[i]))
                                + '_' + XH + str(int(hours[i])) + XM + str(int(minu[i])) + w + '.tif', dRad)

            # if outputSstr:
            #     saveraster(gdal_dsm, outputDir + '/Sstr_' + str(int(YYYY[0, i])) + '_' + str(int(DOY[i]))
            #                     + '_' + XH + str(int(hours[i])) + XM + str(int(minu[i])) + w + '.tif', Sstr)

            # Sky view image of patches
            if ((anisotropic_sky == 1) & (i == 0) & (not poisxy is None)):
                    for k in range(poisxy.shape[0]):
                        Lsky_patch_characteristics[:,2] = patch_characteristics[:,k]
                        skyviewimage_out = outputDir + '/POI_' + str(poiname[k]) + '.png'
                        PolarBarPlot(Lsky_patch_characteristics, altitude[0][i], azimuth[0][i], 'Hemisphere partitioning', skyviewimage_out, 0, 5, 0)

        # Save files for Tree Planter
        if outputTreeplanter:
            feedback.setProgressText("Saving files for Tree Planter tool")
            # Save DSM
            copyfile(filepath_dsm, outputDir + '/DSM.tif')

            # Save met file
            #copyfile(inputMet, outputDir + '/metfile.txt')

            # Save CDSM
            if usevegdem == 1:
                copyfile(filePath_cdsm, outputDir + '/CDSM.tif')

            # Saving settings from SOLWEIG for SOLWEIG1D in TreePlanter
            settingsHeader = 'UTC, posture, onlyglobal, landcover, anisotropic, cylinder, albedo_walls, albedo_ground, emissivity_walls, emissivity_ground, absK, absL, elevation, patch_option'
            settingsFmt = '%i', '%i', '%i', '%i', '%i', '%i', '%1.2f', '%1.2f', '%1.2f', '%1.2f', '%1.2f', '%1.2f', '%1.2f', '%i'
            settingsData = np.array([[utc, pos, onlyglobal, landcover, anisotropic_sky, cyl, albedo_b, albedo_g, ewall, eground, absK, absL, alt, patch_option]])
            np.savetxt(outputDir + '/treeplantersettings.txt', settingsData, fmt=settingsFmt, header=settingsHeader, delimiter=' ')

        # Output I0 vs. Kglobal plot
        radG_for_plot = radG[DOY == unique_days[0]]
        hours_for_plot = hours[DOY == unique_days[0]]
        fig, ax = plt.subplots()
        ax.plot(hours_for_plot, I0_array, label='I0')
        ax.plot(hours_for_plot, radG_for_plot, label='Kglobal')
        ax.set_ylabel('Shortwave radiation [$Wm^{-2}$]')
        ax.set_xlabel('Hours')
        ax.set_title('UTC' + str(int(utc)))
        ax.legend()
        fig.savefig(outputDir + '/metCheck.png', dpi=150)

        # Copying met file for SpatialTC
        copyfile(inputMet, outputDir + '/metforcing.txt')
        
        tmrtplot = tmrtplot / Ta.__len__()  # fix average Tmrt instead of sum, 20191022
        saveraster(gdal_dsm, outputDir + '/Tmrt_average.tif', tmrtplot)
        feedback.setProgressText("SOLWEIG: Model calculation finished.")

        return {self.OUTPUT_DIR: outputDir}
    
    def name(self):
        """
        Returns the algorithm name, used for identifying the algorithm. This
        string should be fixed for the algorithm, and must not be localised.
        The name should be unique within each provider. Names should contain
        lowercase alphanumeric characters only and no spaces or other
        formatting characters.
        """
        return 'Outdoor Thermal Comfort: SOLWEIG'

    def displayName(self):
        """
        Returns the translated algorithm name, which should be used for any
        user-visible display of the algorithm name.
        """
        return self.tr('Outdoor Thermal Comfort: SOLWEIG v2022a')

    def group(self):
        """
        Returns the name of the group this algorithm belongs to. This string
        should be localised.
        """
        return self.tr(self.groupId())

    def groupId(self):
        """
        Returns the unique ID of the group this algorithm belongs to. This
        string should be fixed for the algorithm, and must not be localised.
        The group id should be unique within each provider. Group id should
        contain lowercase alphanumeric characters only and no spaces or other
        formatting characters.
        """
        return 'Processor'

    def shortHelpString(self):
        return self.tr('SOLWEIG (v2022a) is a model which can be used to estimate spatial variations of 3D radiation fluxes and '
                       'mean radiant temperature (Tmrt) in complex urban settings. The SOLWEIG model follows the same '
                       'approach commonly adopted to observe Tmrt, with shortwave and longwave radiation fluxes from  '
                       'six directions being individually calculated to derive Tmrt. The model requires a limited number '
                       'of inputs, such as direct, diffuse and global shortwave radiation, air temperature, relative '
                       'humidity, urban geometry and geographical information (latitude, longitude and elevation). '
                       'Additional vegetation and ground cover information can also be used to imporove the estimation of Tmrt.\n'
                       '\n'
                       'Tools to generate sky view factors, wall height and aspect etc. is available in the pre-processing past in UMEP\n'
                       '\n'
                       '------------\n'
                       '\n'
                       'Full manual available via the <b>Help</b>-button.')

    def helpUrl(self):
        url = "https://umep-docs.readthedocs.io/en/latest/processor/Outdoor%20Thermal%20Comfort%20SOLWEIG.html"
        return url

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

    def icon(self):
        cmd_folder = Path(os.path.split(inspect.getfile(inspect.currentframe()))[0]).parent
        icon = QIcon(str(cmd_folder) + "/icons/icon_solweig.png")
        return icon

    def createInstance(self):
        return ProcessingSOLWEIGAlgorithm()
