from __future__ import absolute_import
from builtins import str
from PyQt5.QtCore import *
from PyQt5.QtGui import *
from PyQt5 import QtCore, QtGui
from PyQt5.QtWidgets import *
from qgis.core import *
from qgis.gui import *
from .Ui_RUSLE_slope import Ui_RUSLE_slope
from . import GdalTools_utils as Utils
from osgeo import gdal
import os
import osgeo
try:
    from PyQt5.QtCore import QString
except ImportError:
    QString = str
from processing.core.Processing import Processing
from . import ftools_utils
from . import GdalTools_utils as Utils
import numpy as np

class RUSLE_slope(QDialog, Ui_RUSLE_slope):
    
    def __init__(self, iface):
        QDialog.__init__(self)
        self.iface = iface
        self.setupUi(self)
        
        self.selectButton_contour.clicked.connect(self.fillInputFileEdit)
        self.inputLayerCombo_contour.currentIndexChanged.connect(self.fillInputAttrib)
        self.selectButton_dem.clicked.connect(self.fillInputRasterEdit)
        self.selectButton_out.clicked.connect( self.fillOutputFileEdit) 
        self.buttonBox.accepted.connect( self.convert)
        self.selectButton_mask.clicked.connect(self.fillInputFileEdit2)
        
        # connect help
        #self.buttonBox.button(QtGui.QDialogButtonBox.Help).clicked.connect( self.help)
        
    # HELP
    def help(self):
        QMessageBox.about(self, "RUSLE", """<p><b>LS factor</b></p> 
        <p><b>Definition</b>: LS factor accounts for the effect of topography on erosion (Renard et al., 1997). The length factor (L) is a ratio 
        which compares the soil loss with that from a field of specified length of 22.6 m (Wischmeier et al., 1978). The slope factor (S) is a 
        ratio which compares the soil loss with that from a field of specified slope of 9%. The factors L and S are usually evaluated together 
        in the erosion prediction. These LS values could be obtained directly in the field or through contour maps or Digital Elevation Models 
        (DEM) (Renard et al., 1997). There are different equations to calculate LS factor. Wischmeier and Smith (1978) proposed that slope length 
        and slope steepness can be used in a single index, which expresses the ratio of soil loss.</p>
        <p><b>Method</b></p> 
        <p>Input file = Contour lines (shapefile) and Mask file OR DEM (raster) and Mask file. </p>
        <p>Parameters = In case of contours define Attribute with elevation values and cell size. </p>
        <p>Output file = LS factor raster file </p>""")          

    # INPUT VECTOR FILE
    def fillInputFileEdit(self):
        lastUsedFilter = Utils.FileFilter.lastUsedVectorFilter()
        inputFile, encoding = Utils.FileDialog.getOpenFileName(self, self.tr( "Select the input file" ), Utils.FileFilter.allVectorsFilter(), lastUsedFilter, True)
        self.lastEncoding = encoding
        self.loadFields(inputFile)
        self.inputLayerCombo_contour.addItem(inputFile)
        check = QFile(inputFile)
    
    def fillInputAttrib(self, inputFile):
        self.layer = QgsVectorLayer(self.inputLayerCombo_contour.currentText(), 'inputFile'  , "ogr")
        self.lineAttrib.clear()    
        changedField = ftools_utils.getFieldList(self.layer)
        for f in changedField:
            if f.type() == QVariant.Int or f.type() == QVariant.String or f.type()==QVariant.Double:
                self.lineAttrib.addItem(str(f.name()))        
    
    def loadFields(self, vectorFile = QString()):
        if vectorFile == None:
            return
                
        try:
            (fields, names) = Utils.getVectorFields(vectorFile)
        except Exception as e:
            QErrorMessage(self).showMessage( str(e) )
            
        self.inputLayerCombo_contour.clearEditText()
        ncodec = QTextCodec.codecForName(self.lastEncoding)
        
    def fillInputFileEdit2(self):
        lastUsedFilter = Utils.FileFilter.lastUsedVectorFilter()
        inputFile, encoding = Utils.FileDialog.getOpenFileName(self, self.tr( "Select the input file" ), Utils.FileFilter.allVectorsFilter(), lastUsedFilter, True)
        self.lastEncoding = encoding
        self.loadFields(inputFile)
        self.inputLayerCombo_mask.addItem(inputFile)
        check = QFile(inputFile) 
          
    # ------------------------------ // ------------------------------------ // ----------------------------
    
    # INPUT RASTER FILE
    def fillInputRasterEdit(self):
        lastUsedFilter = Utils.FileFilter.lastUsedRasterFilter()
        inputFile, encoding = Utils.FileDialog.getOpenFileName(self, self.tr( "Select the input DEM" ), Utils.FileFilter.allRastersFilter(), lastUsedFilter, True)
        self.lastEncoding = encoding
        #self.loadFields(inputFile)
        self.inputLayerCombo_dem.addItem(inputFile)
        #check = QFile(inputFile)    
            
    # OUTPUT RASTER FILE
    def fillOutputFileEdit(self):
        lastUsedFilter = Utils.FileFilter.lastUsedRasterFilter()
        outputFile = Utils.FileDialog.getSaveFileName(self, self.tr( "Select the raster file to save the results to" ), Utils.FileFilter.allRastersFilter(), lastUsedFilter )
        if outputFile==None:
            return
        Utils.FileFilter.setLastUsedRasterFilter(lastUsedFilter)
                
        #self.outputFormat = Utils.fillRasterOutputFormat( lastUsedFilter, outputFile )
        self.inputLayerCombo_out.setText(outputFile)   
                
    # -------------------------------- // --------------------------------- // -------------------------------
    
    # CONVERT VECTOR TO GRID
    def convert(self):
        self.buttonBox.button(QDialogButtonBox.Ok).setDefault(False)
        # if contour exist
        if self.inputLayerCombo_contour.currentText()!="":
            inputLayer = self.inputLayerCombo_contour.currentText()
            inputMask = self.inputLayerCombo_mask.currentText()
            # layer information
            layer = QgsVectorLayer(inputLayer, inputLayer , "ogr")
            vectorlayer_vector =  layer.dataProvider()
            # extent
            extent_rect = vectorlayer_vector.extent()
            xmin = extent_rect.xMinimum()
            xmax = extent_rect.xMaximum()
            ymin = extent_rect.yMinimum()
            ymax = extent_rect.yMaximum()
            extent = str(xmin) + "," + str(xmax) + "," + str(ymin) + "," + str(ymax)
            # elevation attribute
            Elevation = self.lineAttrib.currentText()
            # cellsize
            cellSize = int(self.linePix.value())
            outPath = self.inputLayerCombo_out.text()
          
            Processing.initialize()
            # grid directory (qgis2)
            filedir = QFileInfo(QgsApplication.qgisUserDatabaseFilePath()).path() + "/raster_rusle.tif"
           
            
            #Processing.runAlgorithm("grass7:v.to.rast.attribute", None, inputLayer, 0,Elevation, extent, cellSize, -1.0, 0.0001, filedir)
            Processing.runAlgorithm("grass7:v.to.rast", {
                'input': inputLayer,
                'type': [0, 1, 3], 'where': '', 'use': 0, 'attribute_column': Elevation, 'rgb_column': None,
                'label_column': None, 'value': 1, 'memory': 300, 'output': filedir,
                'GRASS_REGION_PARAMETER': None, 'GRASS_REGION_CELLSIZE_PARAMETER': cellSize,
                'GRASS_RASTER_FORMAT_OPT': '',
                'GRASS_RASTER_FORMAT_META': '', 'GRASS_SNAP_TOLERANCE_PARAMETER': -1,
                'GRASS_MIN_AREA_PARAMETER': 0.0001})
            out = filedir
            userDir = QFileInfo(QgsApplication.qgisUserDatabaseFilePath()).path() + "/grid_topo_rusle.tif"
            Processing.runAlgorithm("grass:r.surf.contour", None, out, extent, cellSize, userDir)
            Processing.runAlgorithm("grass7:r.surf.contour", {'input': out,
                                                     'output': userDir, 'GRASS_REGION_PARAMETER': extent + ' [EPSG:3763]',
                                                     'GRASS_REGION_CELLSIZE_PARAMETER': cellSize,
                                                     'GRASS_RASTER_FORMAT_OPT': '', 'GRASS_RASTER_FORMAT_META': ''})
            outGrid = userDir
            
            # clip grid(interpolation) with polygon (mask)
            maskgrid = QFileInfo(QgsApplication.qgisUserDatabaseFilePath()).path() + "/maskgrid_rusle.tif"
            #Processing.runAlgorithm("saga:clipgridwithpolygon", None, outGrid, inputMask, maskgrid)
            Processing.runAlgorithm("gdal:cliprasterbymasklayer", {'INPUT': outGrid,
                                                                   'MASK': inputMask,
                                                                   'SOURCE_CRS': None, 'TARGET_CRS': None,
                                                                   'NODATA': None,
                                                                   'ALPHA_BAND': False, 'CROP_TO_CUTLINE': True,
                                                                   'KEEP_RESOLUTION': False,
                                                                   'SET_RESOLUTION': False,
                                                                   'X_RESOLUTION': None, 'Y_RESOLUTION': None,
                                                                   'MULTITHREADING': False, 'OPTIONS': '',
                                                                   'DATA_TYPE': 0,
                                                                   'EXTRA': '', 'OUTPUT': maskgrid})
            
            list = []
            
            # slope steepness (percentage)
            grid_mask = maskgrid
            userSlope = QFileInfo(QgsApplication.qgisUserDatabaseFilePath()).path() + "/slope_steepness.sdat"
            #Processing.runAlgorithm("grass:r.slope.aspect",None, grid_mask, 1, 1, 1.0, 0.0, extent, cellSize, userSlope, None, None, None, None, None, None, None, None)
            Processing.runAlgorithm("saga:slopeaspectcurvature",
                                    {'ELEVATION': grid_mask, 'METHOD': 6, 'UNIT_SLOPE': 2,
                                     'UNIT_ASPECT': 0, 'SLOPE': userSlope, 'ASPECT': 'TEMPORARY_OUTPUT',
                                     'C_GENE': 'TEMPORARY_OUTPUT', 'C_PLAN': 'TEMPORARY_OUTPUT',
                                     'C_PROF': 'TEMPORARY_OUTPUT',
                                     'C_TANG': 'TEMPORARY_OUTPUT', 'C_LONG': 'TEMPORARY_OUTPUT',
                                     'C_CROS': 'TEMPORARY_OUTPUT',
                                     'C_MINI': 'TEMPORARY_OUTPUT', 'C_MAXI': 'TEMPORARY_OUTPUT',
                                     'C_TOTA': 'TEMPORARY_OUTPUT',
                                     'C_ROTO': 'TEMPORARY_OUTPUT'})
            slope_steepness = userSlope
            
            #list.append(str(slope_steepness))
            
            # slope lenght (meters)
            slope_lenght = QFileInfo(QgsApplication.qgisUserDatabaseFilePath()).path() + "/slope_lenght.sdat"
            #Processing.runAlgorithm("saga:slopelength", None, grid_mask, slope_lenght)
            Processing.runAlgorithm("saga:slopelength", {
                'DEM': grid_mask,
                'LENGTH': slope_lenght})
            slope_lenght_tif = slope_lenght
           
            list.append(str(slope_lenght_tif))
            a = ";".join(list)
            
            
            # LS factor
            #Processing.runAlgorithm("saga:rastercalculator", None, slope_steepness, ";".join(list), "(0.065+0.0456*a + 0.006541*(a^2)) * ((b/22.1)^0.5)", outPath) 
            #Processing.runAlgorithm("grass:r.mapcalculator",None,slope_steepness,slope_lenght_tif,None,None,None,None,"(0.065+0.0456*A + 0.006541*(A^2)) * ((B/22.1)^0.5)",extent,cellSize,outPath)
            Processing.runAlgorithm("grass7:r.mapcalc.simple", {
                'a': slope_steepness,
                'b': slope_lenght_tif,
                'c': None, 'd': None, 'e': None, 'f': None,
                'expression': '(0.065+0.0456*A + 0.006541*(A^2)) * ((B/22.1)^0.5)', 'output': outPath,
                'GRASS_REGION_PARAMETER': extent + ' [EPSG:3763]',
                'GRASS_REGION_CELLSIZE_PARAMETER': cellSize, 'GRASS_RASTER_FORMAT_OPT': '',
                'GRASS_RASTER_FORMAT_META': ''})

            # add result into canvas
            file_info = QFileInfo(outPath)
            if file_info.exists():
                layer_name = file_info.baseName()
            else:
                return False
            rlayer_new = QgsRasterLayer(outPath, layer_name)
            # add result into canvas
            QgsProject.instance().addMapLayer(rlayer_new)
            self.canvas = QgsMapCanvas()
            # self.canvas.scene().addItem(vdlayer)
            QgsProject.instance().addMapLayer(rlayer_new)
 
        # if dem exist
        if self.inputLayerCombo_contour.currentText()=="":
            gdal.AllRegister() 
            inputLayer_dem = self.inputLayerCombo_dem.currentText()
            inputMask = self.inputLayerCombo_mask.currentText()
           
            gdalRaster = gdal.Open(str(inputLayer_dem))
            
            x = gdalRaster.RasterXSize
            y = gdalRaster.RasterYSize
            geo = gdalRaster.GetGeoTransform()  
            minx = geo[0]
            maxy = geo[3]
            maxx = minx + geo[1]*x
            miny = maxy + geo[5]*y
            extent_raster = str(minx) + "," + str(maxx) + "," + str(miny) + "," + str(maxy)  
            pixelSize = geo[1]   
            inputMask = self.inputLayerCombo_mask.currentText()
            outPath = self.inputLayerCombo_out.text()
            
            list_dem = []
            
            Processing.initialize()
        
            # slope steepness from DEM
            userSlope_dem = QFileInfo(QgsApplication.qgisUserDatabaseFilePath()).path() + "/slope_steepness_dem_rusle.sdat"
            #Processing.runAlgorithm("grass:r.slope.aspect",None, inputLayer_dem, 1, 1, 1.0, 0.0, extent_raster, pixelSize, userSlope_dem, None, None, None, None, None, None, None, None)       
            #Processing.runAlgorithm("saga:slopeaspectcurvature", None, inputLayer_dem, 6,2,0, userSlope_dem, None, None, None, None, None, None, None, None,None,None,None)
            Processing.runAlgorithm("saga:slopeaspectcurvature",
                           {'ELEVATION': inputLayer_dem, 'METHOD': 6, 'UNIT_SLOPE': 2,
                            'UNIT_ASPECT': 0, 'SLOPE': userSlope_dem, 'ASPECT': 'TEMPORARY_OUTPUT',
                            'C_GENE': 'TEMPORARY_OUTPUT', 'C_PLAN': 'TEMPORARY_OUTPUT', 'C_PROF': 'TEMPORARY_OUTPUT',
                            'C_TANG': 'TEMPORARY_OUTPUT', 'C_LONG': 'TEMPORARY_OUTPUT', 'C_CROS': 'TEMPORARY_OUTPUT',
                            'C_MINI': 'TEMPORARY_OUTPUT', 'C_MAXI': 'TEMPORARY_OUTPUT', 'C_TOTA': 'TEMPORARY_OUTPUT',
                            'C_ROTO': 'TEMPORARY_OUTPUT'})
            outSlope_dem = userSlope_dem
    
                       
            # slope lenght (meters) from DEM
            slope_lenght_dem = QFileInfo(QgsApplication.qgisUserDatabaseFilePath()).path() + "/slope_lenght_dem.sdat"
            #Processing.runAlgorithm("saga:slopelength", None, inputLayer_dem, slope_lenght_dem)
            Processing.runAlgorithm("saga:slopelength", {
                'DEM': inputLayer_dem,
                'LENGTH': slope_lenght_dem})
            slope_lenght_tif_dem = slope_lenght_dem
           
            
            # LS factor
            ls_factor = QFileInfo(QgsApplication.qgisUserDatabaseFilePath()).path() + "/ls_factor_factor.tif"
            #Processing.runAlgorithm("saga:rastercalculator", None, userSlope_dem, slope_lenght_dem, "(0.065+0.0456*a + 0.006541*(a^2)) * ((b/22.1)^0.5)", False,8, ls_factor) 
            #Processing.runAlgorithm("grass:r.mapcalculator",None,outSlope_dem, slope_lenght_tif_dem,None,None,None,None,"(0.065+0.0456*A + 0.006541*(A^2)) * ((B/22.1)^0.5)",extent_raster,pixelSize,ls_factor)
            Processing.runAlgorithm("grass7:r.mapcalc.simple", {
                'a': outSlope_dem,
                'b': slope_lenght_tif_dem,
                'c': None, 'd': None, 'e': None, 'f': None,
                'expression': '(0.065+0.0456*A + 0.006541*(A^2)) * ((B/22.1)^0.5)', 'output': ls_factor,
                'GRASS_REGION_PARAMETER': extent_raster + ' [EPSG:3763]',
                'GRASS_REGION_CELLSIZE_PARAMETER': pixelSize, 'GRASS_RASTER_FORMAT_OPT': '', 'GRASS_RASTER_FORMAT_META': ''})
            ls_factor_factor = ls_factor
            
            # clip grid(interpolation) with polygon (mask)
            #Processing.runAlgorithm("saga:clipgridwithpolygon", None, ls_factor_factor, inputMask, outPath)
            Processing.runAlgorithm("gdal:cliprasterbymasklayer", {'INPUT': ls_factor_factor,
                                                          'MASK': inputMask,
                                                          'SOURCE_CRS': None, 'TARGET_CRS': None, 'NODATA': None,
                                                          'ALPHA_BAND': False, 'CROP_TO_CUTLINE': True,
                                                          'KEEP_RESOLUTION': False, 'SET_RESOLUTION': False,
                                                          'X_RESOLUTION': None, 'Y_RESOLUTION': None,
                                                          'MULTITHREADING': False, 'OPTIONS': '', 'DATA_TYPE': 0,
                                                          'EXTRA': '', 'OUTPUT': outPath})
                        
            # LS factor
            #Processing.runAlgorithm("saga:rastercalculator", None, ";".join(list_dem), "(0.065+0.0456*a + 0.006541*(a^2)) * ((b/22.1)^0.5)", outPath)         

            # add result into canvas
            file_info = QFileInfo(outPath)
            if file_info.exists():
                layer_name = file_info.baseName()
            else:
                return False
            rlayer_new = QgsRasterLayer(outPath, layer_name)
            # add result into canvas
            QgsProject.instance().addMapLayer(rlayer_new)
            self.canvas = QgsMapCanvas()
            # self.canvas.scene().addItem(vdlayer)
            QgsProject.instance().addMapLayer(rlayer_new)