# -*- 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, QVariant
from qgis.core import (QgsProcessing,
                       QgsProcessingAlgorithm,
                       QgsProcessingParameterPoint,
                       QgsProcessingParameterString,
                       QgsProcessingParameterBoolean,
                       QgsProcessingParameterNumber,
                       QgsProcessingParameterFolderDestination,
                       QgsProcessingParameterFeatureSink,
                       QgsProcessingParameterRasterLayer,
                       QgsProcessingParameterEnum,
                       QgsProcessingParameterFeatureSource,
                       QgsProcessingParameterVectorDestination,
                       QgsProcessingParameterExtent,
                       QgsProcessingException,
                       QgsVectorLayer,
                       QgsFeature,
                       QgsGeometry,
                       QgsPointXY,
                       QgsVectorFileWriter)

from qgis.PyQt.QtGui import QIcon
from osgeo import gdal, osr
from osgeo.gdalconst import *
import os
import numpy as np
import inspect
from pathlib import Path
import sys
from ..util import misc
from ..util import landCoverFractions_v1 as land


class ProcessingLandCoverFractionPointAlgorithm(QgsProcessingAlgorithm):
    """
    This algorithm is a processing version of Image Morphometric Calculator Point
    """

    USE_POINTLAYER = 'USE_POINTLAYER'
    INPUT_POINT = 'INPUT_POINT'
    INPUT_POINTLAYER = 'INPUT_POINTLAYER'
    INPUT_DISTANCE = 'INPUT_DISTANCE'
    INPUT_INTERVAL = 'INPUT_INTERVAL'
    INPUT_LCGRID = 'INPUT_LCGRID'
    FILE_PREFIX = 'FILE_PREFIX'
    OUTPUT_DIR = 'OUTPUT_DIR'
    OUTPUT_POLYGON = 'OUTPUT_POLYGON'
    SAVE_POINT = 'SAVE_POINT'
    OUTPUT_POINT = 'OUTPUT_POINT'
    
    def initAlgorithm(self, config):
        
        self.addParameter(QgsProcessingParameterPoint(self.INPUT_POINT,
            self.tr('Point of interest'), optional=True))
        # self.addParameter(QgsProcessingParameterBoolean(self.USE_POINTLAYER,
        #     self.tr("Obtain point of interest from point in vector layer"), defaultValue=False))
        self.addParameter(QgsProcessingParameterFeatureSource(self.INPUT_POINTLAYER,
            self.tr('Point vector layer'), [QgsProcessing.TypeVectorPoint], optional=True))
        self.addParameter(QgsProcessingParameterNumber(self.INPUT_DISTANCE, 
            self.tr('Search distance (meter)'),
            QgsProcessingParameterNumber.Integer,
            QVariant(200), False, minValue=0))
        self.addParameter(QgsProcessingParameterNumber(self.INPUT_INTERVAL, 
            self.tr('Wind direction search interval (degree)'), 
            QgsProcessingParameterNumber.Double,
            QVariant(5), False, minValue=0.1, maxValue=360.))
        self.addParameter(QgsProcessingParameterRasterLayer(self.INPUT_LCGRID,
            self.tr('UMEP formatted land cover grid (see help for more info)'), None, False))
        self.addParameter(QgsProcessingParameterString(self.FILE_PREFIX, 
            self.tr('File prefix')))
        self.addParameter(QgsProcessingParameterBoolean(self.SAVE_POINT,
            self.tr("Save point of interest as new vector layer"), defaultValue=False))
        self.addParameter(QgsProcessingParameterVectorDestination(self.OUTPUT_POINT,
            self.tr("Output point layer (point of interest)")))
        self.addParameter(QgsProcessingParameterVectorDestination(self.OUTPUT_POLYGON,
            self.tr("Output polygon layer (area of interest)")))
        self.addParameter(QgsProcessingParameterFolderDestination(self.OUTPUT_DIR, 
            self.tr('Output folder')))

        self.plugin_dir = os.path.dirname(__file__)
        if not (os.path.isdir(self.plugin_dir + '/data')):
            os.mkdir(self.plugin_dir + '/data')

    def processAlgorithm(self, parameters, context, feedback):
        # InputParameters
        # usePointlayer = self. parameterAsBool(parameters, self.USE_POINTLAYER, context)
        inputPoint = None
        inputPointLayer = self.parameterAsVectorLayer(parameters, self.INPUT_POINTLAYER, context)
        inputDistance = self.parameterAsDouble(parameters, self.INPUT_DISTANCE, context)
        inputInterval = self.parameterAsDouble(parameters, self.INPUT_INTERVAL, context)
        lcgrid = self.parameterAsRasterLayer(parameters, self.INPUT_LCGRID, context)
        filePrefix = self.parameterAsString(parameters, self.FILE_PREFIX, context)
        outputDir = self.parameterAsString(parameters, self.OUTPUT_DIR, context)
        outputPolygon = self.parameterAsOutputLayer(parameters, self.OUTPUT_POLYGON, context)
        savePoint = self.parameterAsBool(parameters, self.SAVE_POINT, context)
        outputPoint = None
        
        if parameters['OUTPUT_DIR'] == 'TEMPORARY_OUTPUT':
            if not (os.path.isdir(outputDir)):
                os.mkdir(outputDir)

        # Get POI
        if inputPointLayer is None:
            feedback.setProgressText("Point location obtained manually")
            inputPoint = self.parameterAsPoint(parameters, self.INPUT_POINT, context)
            x = float(inputPoint[0])
            y = float(inputPoint[1])
        else:
            # inputPointLayer = self.parameterAsVectorLayer(parameters, self.INPUT_POINTLAYER, context)
            feedback.setProgressText("Point location obtained from point in vector layer")
            if not inputPointLayer.selectedFeatures():
                if inputPointLayer.featureCount() != 1:
                    raise QgsProcessingException("Point vector layer includes more than one object. Select one feature point and try again.")
                else:
                    for poi in inputPointLayer.getFeatures():
                        x = poi.geometry().asPoint()[0]
                        y = poi.geometry().asPoint()[1]
            else:
                if inputPointLayer.selectedFeatureCount() != 1:
                    raise QgsProcessingException("More than one point of interest is selected. Select only one feature point and try again.")
                for poi in inputPointLayer.selectedFeatures():
                    x = poi.geometry().asPoint()[0]
                    y = poi.geometry().asPoint()[1]

        feedback.setProgressText("x = " + str(x))
        feedback.setProgressText("y = " + str(y))
                    
        r = inputDistance
        
        dsmlayer = self.parameterAsRasterLayer(parameters, self.INPUT_LCGRID, context)
        if dsmlayer is None:
            raise QgsProcessingException("No valid building land cover raster layer is selected")

        provider = dsmlayer.dataProvider()
        filePath_dsm_build = str(provider.dataSourceUri())
        bigraster = gdal.Open(filePath_dsm_build)
        bbox = (x - r, y + r, x + r, y - r)
        gdal.Translate(self.plugin_dir + '/data/clipdsm.tif', bigraster, projWin=bbox)
        bigraster = None
        dataset = gdal.Open(self.plugin_dir + '/data/clipdsm.tif')
        dsm = dataset.ReadAsArray().astype(float)
        sizex = dsm.shape[0]
        sizey = dsm.shape[1]

        geotransform = dataset.GetGeoTransform()
        scale = 1 / geotransform[1]
        degree = float(inputInterval)
        nd = dataset.GetRasterBand(1).GetNoDataValue()
        nodata_test = (dsm == nd)
        
        if nodata_test.any():
            raise QgsProcessingException("NoData values present within the radius (" + str(int(r)) + 
                                " m) from point of interest. Extend your raster(s) or more point.")
        else:
            # immorphresult = morph.imagemorphparam_v2(dsm, dem, scale, 1, degree, feedback, 1)
            landcoverresult = land.landcover_v1(dsm, 1, degree, feedback, 1)

        


        # save to file
        pre = filePrefix
        header = 'Wd Paved Buildings EvergreenTrees DecidiousTrees Grass Baresoil Water'
        numformat = '%3d %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f'
        arr = np.concatenate((landcoverresult["deg"], landcoverresult["lc_frac"]), axis=1)
        np.savetxt(outputDir + '/' + pre + '_' + 'LCFPoint_anisotropic.txt', arr,
                   fmt=numformat, delimiter=' ', header=header, comments='')

        header = 'Paved Buildings EvergreenTrees DecidiousTrees Grass Baresoil Water'
        numformat = '%5.3f %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f'
        arr2 = np.array(landcoverresult["lc_frac_all"])
        np.savetxt(outputDir + '/' + pre + '_' + 'LCFPoint_isotropic.txt', arr2,
                    fmt=numformat, delimiter=' ', header=header, comments='')

        dataset = None
        dataset2 = None
        dataset3 = None
        
        crs = str(dsmlayer.crs().authid())
        self.create_poly_layer(outputPolygon, x, y, r, crs)
        
        if savePoint:
            outputPoint = self.parameterAsOutputLayer(parameters, self.OUTPUT_POINT, context)
            self.create_point_layer(outputPoint, x, y, crs)
            results = {self.OUTPUT_DIR: outputDir, self.OUTPUT_POLYGON: outputPolygon, self.OUTPUT_POINT: outputPoint}
        else:
            results = {self.OUTPUT_DIR: outputDir, self.OUTPUT_POLYGON: outputPolygon}

        return results

    def create_point_layer(self, outputPoint, x, y, crs):
        uri = "Point?field=id:integer&field=x:double&field=y:double&index=yes&crs=" + crs
        self.poiLayer = QgsVectorLayer(uri, "Point of Interest", "memory")
        self.provider = self.poiLayer.dataProvider()
        # create the feature
        fc = int(self.provider.featureCount())
        feature = QgsFeature()
        feature.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(x, y)))
        feature.setAttributes([fc, x, y])
        self.poiLayer.startEditing()
        self.poiLayer.addFeature(feature)  # ,True
        self.poiLayer.commitChanges()

        fields = self.poiLayer.fields()

        writer = QgsVectorFileWriter(outputPoint, "CP1250", fields, self.provider.wkbType(),
                                     self.provider.crs(), "ESRI shapefile")
        self.poiLayer.selectAll()
        selection = self.poiLayer.selectedFeatures()

        for feature in selection:
            writer.addFeature(feature)
        del writer
    
    def create_poly_layer(self, outputPolygon, x, y, r, crs):
        uri = "Polygon?field=id:integer&index=yes&crs=" + crs
        self.polyLayer = QgsVectorLayer(uri, "Study area", "memory")
        self.provider = self.polyLayer.dataProvider()

        # create buffer feature
        fc = int(self.provider.featureCount())
        featurepoly = QgsFeature()

        # Assign feature the buffered geometry
        radius = r
        featurepoly.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(x, y)).buffer(radius, 1000)) # fix issue #400
        featurepoly.setAttributes([fc])
        self.polyLayer.startEditing()
        self.polyLayer.addFeature(featurepoly)
        self.polyLayer.commitChanges()
        prov = self.provider
        fields = self.polyLayer.fields()

        writer = QgsVectorFileWriter(outputPolygon, "CP1250", fields, prov.wkbType(),
                                     prov.crs(), "ESRI shapefile")
        self.polyLayer.selectAll()
        selection = self.polyLayer.selectedFeatures()

        for feature in selection:
            writer.addFeature(feature)
        del writer
    
    def name(self):
        return 'Urban Land Cover: Land Cover Fraction (Point)'

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

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

    def groupId(self):
        return 'Pre-Processor'

    def shortHelpString(self):
        return self.tr('The Land Cover Fraction (Point) plugin calculates land cover fractions required for UMEP (see Land Cover '
        'Reclassifier) from a point location based on a land cover raster grid. A land cover grid suitable for the processor in '
        'UMEP can be derived using the Land Cover Classifier. The fraction will vary depending on what angle (wind direction) '
        'you are interested in. Thus, this plugin is able to derive the land cover fractions for different directions.\n'
        '-------------\n'
        'Full manual available via the <b>Help</b>-button.')

    def helpUrl(self):
        url = "https://umep-docs.readthedocs.io/en/latest/pre-processor/Urban%20Land%20Cover%20Land%20Cover%20Fraction%20(Point).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/LandCoverFractionPointIcon.png")
        return icon

    def createInstance(self):
        return ProcessingLandCoverFractionPointAlgorithm()