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

class RUSLE_runoff(QDialog, Ui_RUSLE_runoff):
    
    def __init__(self, iface):
        QDialog.__init__(self)
        self.iface = iface
        self.setupUi(self)
        
        self.selectButton_meteo.clicked.connect(self.fillInputFileEdit)
        self.selectMask.clicked.connect(self.fillInputMask)        
        self.buttonBox.accepted.connect( self.convert) 
        self.selectButton_out.clicked.connect( self.fillOutputFileEdit)
        self.inputLayerCombo_meteo.currentIndexChanged(QString).connect(self.fillInputAttrib)
        
        # connect help
        #self.buttonBox.button(QtGui.QDialogButtonBox.Help).clicked.connect( self.help)
        
    #   HELP
    def help(self):
        QMessageBox.about(self, "RUSLE", """<p><b>R factor</b></p> 
        <p><b>Definition:</b> Wischmeier and Smith (1978) defines rainfall runoff erosivity factor as the factor that indi-cates the soil loss potential of a given storm event and represents the effect of rainfall intensi-ty on soil erosion. When factors other than rainfall are held constant, soil losses from fields could be computed by the erosion index (EI30). The EI30 is the summation for each rain-storm of the kinetic energy (expressed in MJ/ha), multiplied by the greatest amount of rain in any 30 minute period expressed in mm/h (Wischmeier et al., 1971).  However, the EI30 value is usually not available, which limits the spatio-temporal application of the factor R. In order to overcome this problem, this factor has been estimated from several algorithms based on the monthly rainfall. 
        In this application, the R factor only incorporates one method, rainfall derived EI30, accord-ing to Loureiro and Coutinho (2001). The R factor is composed by two scripts (main code and interface window). To obtain the R factor, two methods are incorporated. The R factor can be obtained through field observa-tions or can be calculated from meteorological data. In this case, the input file is the meteoro-logical data shapefile (point type).
        In order to obtain the R factor, the points are interpolated through an interpolation method chosen by the user, such as: Inverse Distance Weighting (Bonham-Carter, 1994), Kriging (Olaya, 2004), Cubic Spline Approximation (SAGA) (Haber et al., 2001) and Spatial approx-imation using spline with tension (GRASS) (GRASS GIS, 2013).
        In order to define the study area extents, a shapefile input mask is defined.</p>
        <p><b>Method</b></p> 
        <p>Input file = Meteorological data (points file with minimal number of points of 3); Mask file. </p>
        <p>Parameters = Interpolation method; Field to define the R column OR fields to define Rain10 and Days10 columns; Cell size. </p>
        <p>Output file = R factor raster file </p>
        <p>Wischmeier, W.H., Smith, D.D., 1978. Predicting Rainfall Erosion Losses: A Guide to Con-servation Planning with Universal Soil Loss Equation (USLE). Agriculture Handbook, De-partment of Agriculture, Washington, DC, 703.</p>
        <p>Loureiro, N.S., Coutinho, M.A., 2001. A new procedure to estimate the RUSLE EI30 index, based on monthly rainfall data and applied to the Algarve region, Portugal. Journal of Hy-drology. 250, 12-18.</p>
        <p>Bonham-Carter, G.F., 1994. Geographic Information Systems for Geoscientists, Modelling with GIS. Computer Method in Geosciences. 13, 152-153.</p>
        <p>Haber, J., Zeilfelder, F., Davydov, O., Seidel, H.P., 2001. Smooth Approximation and Ren-dering of Large Scattered Data Sets, In: Th Ertl, K Joy and A Varshney (Eds) Proceedings of IEEE Visualization 2001, IEEE Computer Society. 571, 341-347.</p>
        <p>Olaya, V., 2004. A gentle introduction to SAGA GIS, Edition 1.1, Rev. December 9.</p>
        <p>GRASS GIS, 2013. The worlds leading Free GIS software. Available at: http://grass.osgeo.org/ (accessed in: February, 2015).</p>
        <p>Wischmeier, W. H., Johnson C.B., Cross., B.V., 1971. A soil erodibility nomograph for farmland and construction sites. Journal of Soil and Water Conservation. 26, 189 193.</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_meteo.addItem(inputFile)
        check = QFile(inputFile) 
        
    def fillInputAttrib(self, inputFile):
        self.layer = QgsVectorLayer(unicode(self.inputLayerCombo_meteo.currentText()).encode('utf8'), inputFile  , "ogr")
        self.linedays.clear()
        self.linerain.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.linerain.addItem(unicode(f.name()))
                self.linedays.addItem(unicode(f.name()))
                self.lineR.addItem(unicode(f.name()))
                
    def loadFields(self, vectorFile = QString()):
        if vectorFile == None:
            return
        try:
            (fields, names) = Utils.getVectorFields(vectorFile)
        except Exception, e:
            QErrorMessage(self).showMessage( str(e) )          
        self.linerain.clearEditText()
        self.linedays.clearEditText()
        ncodec = QTextCodec.codecForName(self.lastEncoding)  
      
    # INPUT MASK SHAPEFILE
    def fillInputMask(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.inputMaskCombo.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) 
          
     # -------------------------------- // --------------------------------- // -------------------------------
     
# POINTS INTERPOLATION
    def convert(self):    
        
        gdal.AllRegister()
    # ------------------------ FIRST METHOD -------------------------------------------------
        self.buttonBox.button(QDialogButtonBox.Ok).setDefault(False)
        inputLayer = self.inputLayerCombo_meteo.currentText()
        inputMask = self.inputMaskCombo.currentText()
        # layer information
        layer = QgsVectorLayer(unicode(inputLayer).encode('utf8'), inputLayer , "ogr")  
        vectorlayer_vector =  layer.dataProvider()
        layer_mask = QgsVectorLayer(unicode(inputMask).encode('utf8'), inputMask , "ogr")  
        vectorlayer_mask =  layer_mask.dataProvider()        
        # mask extent
        extent_rect = vectorlayer_mask.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)
        # cellsize
        cellSize = int(self.linePix.value())
        outPath = self.inputLayerCombo_out.text() 
        # size of atributte table == number of points
        count = layer.featureCount()
        
        
        if self.comboBoxErosivity.currentText()=='R values':
            R_values = self.lineR.currentText()
            if self.comboBoxMethod.currentText()=='Kriging':
                Processing.initialize()
                # grid directory (qgis2)
                kriging_interpolation_R = QFileInfo(QgsApplication.qgisUserDbFilePath()).path() + "/kriging_interpolation_R.sdat"        
                Processing.runAlgorithm("saga:ordinarykrigingglobal", None, inputLayer, R_values, True, 0, 1, False, 100, False, 0.0, 10, 1000, 1.0, 0.1, 1.0, 0.5, cellSize, True, extent, kriging_interpolation_R, None)
                #kriging_int_R = kriging_interpolation_R + "." + "tif"
                
                mask_R = QFileInfo(QgsApplication.qgisUserDbFilePath()).path() + "/mask_R.sdat"
                # clip grid(interpolation) with polygon (mask)
                Processing.runAlgorithm("saga:clipgridwithpolygon", None, kriging_interpolation_R, inputMask, outPath)
            if self.comboBoxMethod.currentText()=="Inverse Distance Weighting":
                # grid directory (qgis2)
                idw_interpolation_R = QFileInfo(QgsApplication.qgisUserDbFilePath()).path() + "/idw_interpolation_R"        
                Processing.runAlgorithm("grass:v.surf.idw", None, inputLayer, count, 2.0, R_values, False, extent, cellSize, -1.0, 0.0001, idw_interpolation_R)
                idw_int_R = idw_interpolation_R + "." + "tif"
                
                int_mask_R = QFileInfo(QgsApplication.qgisUserDbFilePath()).path() + "/int_mask_R.sdat"
                # clip grid(interpolation) with polygon (mask)
                Processing.runAlgorithm("saga:clipgridwithpolygon", None, idw_int_R, inputMask, outPath)
            if self.comboBoxMethod.currentText()=="Cubic spline approximation (SAGA)":
                # grid directory (qgis2)
                cubicSpline_interpolation_R = QFileInfo(QgsApplication.qgisUserDbFilePath()).path() + "/cubicSpline_interpolation_R.sdat"        
                Processing.runAlgorithm("saga:cubicsplineapproximation", None, inputLayer, R_values, 0, 3, count, 5, 140.0, extent, cellSize, cubicSpline_interpolation_R)
                #cubicSpline_int_R = cubicSpline_interpolation_R + "." + "tif"
                
                # clip grid(interpolation) with polygon (mask)
                Processing.runAlgorithm("saga:clipgridwithpolygon", None, cubicSpline_int_R, inputMask, outPath)
            if self.comboBoxMethod.currentText()=="Spatial approximation using spline with tension (GRASS)":
                # grid directory (qgis2)
                rst_interpolation_R = QFileInfo(QgsApplication.qgisUserDbFilePath()).path() + "/rst_interpolation_R"        
                Processing.runAlgorithm("grass:v.surf.rst", None, inputLayer, "", None, R_values, 40, 40, 300, 0.001, 2.5, 1, 0, 0, False, False, extent, cellSize, -1, 0.0001, rst_interpolation_R, None, None, None, None, None)
                rst_int_R = rst_interpolation_R + "." + "tif"
                
                # clip grid(interpolation) with polygon (mask)
                Processing.runAlgorithm("saga:clipgridwithpolygon", None, rst_int_R, inputMask, outPath)        
 
        # if method == loureiro e coutinho (2001)                    
        if self.comboBoxErosivity.currentText()=='EI30 = 7.05rain(10)-88.92days(10) Loureiro e Coutinho 2001':
            rain_10 = self.linerain.currentText()
            #QMessageBox.about(self,'RUSLE', str(rain_10))
            days_10 = self.linedays.currentText()  
            #QMessageBox.about(self,'RUSLE', str(days_10))

            if self.comboBoxMethod.currentText()=="Kriging":
                # rain kriging
                Processing.initialize()
                # grid directory (qgis2)
                kriging_interpolation_rain = QFileInfo(QgsApplication.qgisUserDbFilePath()).path() + "/kriging_interpolation_rain.sdat"        
                Processing.runAlgorithm("saga:ordinarykrigingglobal", None, inputLayer, rain_10, True, 0, 1, False, 100, False, 0.0, 10, 1000, 1.0, 0.1, 1.0, 0.5, cellSize, True, extent, kriging_interpolation_rain, None)
                #kriging_int_rain = kriging_interpolation_rain + "." + "tif"
                mask_rain = QFileInfo(QgsApplication.qgisUserDbFilePath()).path() + "/mask_rain.sdat"
                # clip grid(interpolation) with polygon (mask)
                Processing.runAlgorithm("saga:clipgridwithpolygon", None, kriging_interpolation_rain, inputMask, mask_rain)
                #int_mask_rain = mask_rain + "." + "tif"    
                # days kriging
                kriging_interpolation_days = QFileInfo(QgsApplication.qgisUserDbFilePath()).path() + "/kriging_interpolation_days.sdat"        
                Processing.runAlgorithm("saga:ordinarykrigingglobal", None, inputLayer, days_10, True, 0, 1, False, 100, False, 0.0, 10, 1000, 1.0, 0.1, 1.0, 0.5, cellSize, True, extent, kriging_interpolation_days, None)
                #kriging_int_days = kriging_interpolation_days + "." + "tif"
                mask_days = QFileInfo(QgsApplication.qgisUserDbFilePath()).path() + "/mask_days.sdat"
                # clip grid(interpolation) with polygon (mask)
                Processing.runAlgorithm("saga:clipgridwithpolygon", None, kriging_int_days, inputMask, mask_days)
                #int_mask_days = mask_days + "." + "tif"  
                
                # Loureiro e Coutinho formula EI30
                #EI30_kriging = QFileInfo(QgsApplication.qgisUserDbFilePath()).path() + "/EI30_kriging"
                Processing.runAlgorithm("grass:r.mapcalculator", None, int_mask_rain, mask_days, None, None, None, None,  "7.05*A - 88.92*B", extent, cellSize, outPath)
                #EI30_tif = EI30 + "." + "tif"             
                
            if self.comboBoxMethod.currentText()=="Inverse Distance Weighting":
                Processing.initialize()
                # rain idw
                # grid directory (qgis2)
                idw_interpolation_rain = QFileInfo(QgsApplication.qgisUserDbFilePath()).path() + "/idw_interpolation_rain"    
                Processing.runAlgorithm("grass:v.surf.idw", None, inputLayer, count, 2.0, rain_10, False, extent, cellSize, -1.0, 0.0001, idw_interpolation_rain)
                idw_int_rain = idw_interpolation_rain + "." + "tif"
                int_mask_rain = QFileInfo(QgsApplication.qgisUserDbFilePath()).path() + "/int_mask_rain.sdat"
                # clip grid(interpolation) with polygon (mask)
                Processing.runAlgorithm("saga:clipgridwithpolygon", None, idw_int_rain, inputMask, int_mask_rain)
                #mask_idw_rain = int_mask_rain + "." + "tif"
                # days idw
                idw_interpolation_days = QFileInfo(QgsApplication.qgisUserDbFilePath()).path() + "/idw_interpolation_days"        
                Processing.runAlgorithm("grass:v.surf.idw", None, inputLayer, count, 2.0, days_10, False, extent, cellSize, -1.0, 0.0001, idw_interpolation_days)
                idw_int_days = idw_interpolation_days + "." + "tif"
                int_mask_days = QFileInfo(QgsApplication.qgisUserDbFilePath()).path() + "/int_mask_days.sdat"
                # clip grid(interpolation) with polygon (mask)
                Processing.runAlgorithm("saga:clipgridwithpolygon", None, idw_int_days, inputMask, int_mask_days)  
                #mask_idw_days = int_mask_days + "." + "tif"
                
                # Loureiro e Coutinho formula EI30
                EI30 = QFileInfo(QgsApplication.qgisUserDbFilePath()).path() + "/EI30"
                Processing.runAlgorithm("saga:rastercalculator", None, int_mask_rain, int_mask_days, "7.05*a - 88.92*b", False, 7, outPath)
                EI30_tif = EI30 + "." + "tif"    
                
                
                
 
            if self.comboBoxMethod.currentText()=="Cubic spline approximation (SAGA)":
                # cubic spline rain
                # grid directory (qgis2)
                cubicSpline_interpolation_rain = QFileInfo(QgsApplication.qgisUserDbFilePath()).path() + "/cubicSpline_interpolation_rain.sdat"        
                Processing.runAlgorithm("saga:cubicsplineapproximation", None, inputLayer, rain_10, 3, 20, 5, 140.0, extent, cellSize, cubicSpline_interpolation_rain)
                #cubicSpline_int_rain = cubicSpline_interpolation_rain + "." + "tif"          
                cubic_mask_rain = QFileInfo(QgsApplication.qgisUserDbFilePath()).path() + "/cubic_mask_rain.sdat"
                # clip grid(interpolation) with polygon (mask)
                Processing.runAlgorithm("saga:clipgridwithpolygon", None, cubicSpline_interpolation_rain, inputMask, cubic_mask_rain)
                #mask_cubic_rain = cubic_mask_rain + "." + "tif"
                # days cubic spline
                cubicSpline_interpolation_days = QFileInfo(QgsApplication.qgisUserDbFilePath()).path() + "/cubicSpline_interpolation_days.sdat"        
                Processing.runAlgorithm("saga:cubicsplineapproximation", None, inputLayer, days_10, 3, 20, 5, 140.0, extent, cellSize, cubicSpline_interpolation_days)
                #cubicSpline_int_days = cubicSpline_interpolation_days + "." + "tif"                
                cubic_mask_days = QFileInfo(QgsApplication.qgisUserDbFilePath()).path() + "/cubic_mask_days.sdat"
                # clip grid(interpolation) with polygon (mask)
                Processing.runAlgorithm("saga:clipgridwithpolygon", None, cubicSpline_interpolation_days, inputMask, cubic_mask_days)  
                #mask_cubic_days = cubic_mask_days + "." + "tif"
                
                # Loureiro e Coutinho formula EI30
                #EI30 = QFileInfo(QgsApplication.qgisUserDbFilePath()).path() + "/EI30"
                Processing.runAlgorithm("grass:r.mapcalculator", None, cubic_mask_rain, cubic_mask_days, None, None, None, None,  "7.05*A - 88.92*B", extent, cellSize, outPath)
                #EI30_tif = EI30 + "." + "tif"                             
              
            if self.comboBoxMethod.currentText()=="Spatial approximation using spline with tension (GRASS)":
                # grid directory (qgis2)
                rst_interpolation_rain = QFileInfo(QgsApplication.qgisUserDbFilePath()).path() + "/rst_interpolation_rain"        
                Processing.runAlgorithm("grass:v.surf.rst", None, inputLayer, "", None, rain_10, 40, 40, 300, 0.001, 2.5, 1, 0, 0, False, False, extent, cellSize, -1, 0.0001, rst_interpolation_rain, None, None, None, None, None)
                rst_int_rain = rst_interpolation_rain + "." + "tif"
                
                # clip grid(interpolation) with polygon (mask)
                spatial_mask_rain = QFileInfo(QgsApplication.qgisUserDbFilePath()).path() + "/spatial_mask_rain.sdat"
                Processing.runAlgorithm("saga:clipgridwithpolygon", None, rst_int_rain, inputMask, spatial_mask_rain)
                #mask_spatial_rain = spatial_mask_rain + "." + "tif"
                
                # days spatial approximation
                rst_interpolation_days = QFileInfo(QgsApplication.qgisUserDbFilePath()).path() + "/rst_interpolation_days"        
                Processing.runAlgorithm("grass:v.surf.rst", None, inputLayer, "", None, days_10, 40, 40, 300, 0.001, 2.5, 1, 0, 0, False, False, extent, cellSize, -1, 0.0001, rst_interpolation_days, None, None, None, None, None)
                rst_int_days = rst_interpolation_days + "." + "tif"   
                
                # clip grid(interpolation) with polygon (mask)
                spatial_mask_days = QFileInfo(QgsApplication.qgisUserDbFilePath()).path() + "/spatial_mask_days.sdat"
                Processing.runAlgorithm("saga:clipgridwithpolygon", None, rst_int_days, inputMask, spatial_mask_days)
                #mask_spatial_days = spatial_mask_days + "." + "tif"
                
                # Loureiro e Coutinho formula EI30
                #EI30 = QFileInfo(QgsApplication.qgisUserDbFilePath()).path() + "/EI30"
                Processing.runAlgorithm("grass:r.mapcalculator", None, spatial_mask_rain, spatial_mask_days, None, None, None, None,  "7.05*A - 88.92*B", extent, cellSize, outPath)
                #EI30_tif = EI30 + "." + "tif"                
 

        # add result into canvas
        
        
        import subprocess
                            
        test = subprocess.call(['gdalinfo.exe', '-stats', str(outPath)], shell=True)   
        QMessageBox.about(self,'RUSLE', str(test))
        
        file_info = QFileInfo(outPath)
        QMessageBox.about(self,'RUSLE', str(file_info))
        if file_info.exists():
            QMessageBox.about(self,'RUSLE', str(file_info))
            layer_name = file_info.baseName()
        else:
            return False        
        dataset = gdal.Open(str(outPath))
        QMessageBox.about(self,'RUSLE', str(outPath))
        QMessageBox.about(self,'RUSLE', str(dataset))
        band = dataset.GetRasterBand(1)
     
        min = band.GetMinimum()
        max = band.GetMaximum() 
        QMessageBox.about(self,'RUSLE', str(min))
        QMessageBox.about(self,'RUSLE', str(max))
        fileInfo = QFileInfo(outPath)
        baseName = fileInfo.baseName()
        rlayer = QgsRasterLayer(outPath, baseName)                
    
        #renderer = rlayer.renderer()
        #provider = rlayer.dataProvider()
        #extent = rlayer.extent()
        
        #ver = provider.hasStatistics(1, QgsRasterBandStats.All)
        #QMessageBox.about(self,'RUSLE', str(ver))
        
        #stats = provider.bandStatistics(1, QgsRasterBandStats.All,extent, 0)
        #QMessageBox.about(self,'RUSLE', str(stats))
        
        #if ver is not False:
            #print "minimumValue = ", stats.minimumValue
        
            #print "maximumValue = ", stats.maximumValue
        
        if (min < 0):
            min = 0  
        
        else: 
            min= min
        
        #max = stats.maximumValue
        range = float(max) - float(min)
        add = range//float(2)
        interval = float(min) + add
        
        colDic = {'red':'#ff0000', 'yellow':'#ffff00','blue':'#0000ff'}
        
        valueList =[float(min), interval, float(max)]
        
        lst = [ QgsColorRampShader.ColorRampItem(valueList[0], QColor(colDic['red'])), 
                QgsColorRampShader.ColorRampItem(valueList[1], QColor(colDic['yellow'])), 
                QgsColorRampShader.ColorRampItem(valueList[2], QColor(colDic['blue']))]
        
        myRasterShader = QgsRasterShader()
        myColorRamp = QgsColorRampShader()
        
        myColorRamp.setColorRampItemList(lst)
        myColorRamp.setColorRampType(QgsColorRampShader.INTERPOLATED)
        myRasterShader.setRasterShaderFunction(myColorRamp)
        
        myPseudoRenderer = QgsSingleBandPseudoColorRenderer(rlayer.dataProvider(), rlayer.type(), myRasterShader)
        
        rlayer.setRenderer(myPseudoRenderer)
        
        rlayer.triggerRepaint()                
        
        QgsMapLayerRegistry.instance().addMapLayer(rlayer)        



        
        
        #file_info = QFileInfo(outPath)
        #if file_info.exists():
            #layer_name = file_info.baseName()
        #else:
            #return False
        #rlayer_new = QgsRasterLayer(outPath, layer_name)
        #if rlayer_new.isValid():
            #QgsMapLayerRegistry.instance().addMapLayer(rlayer_new)
            ##layer = QgsMapCanvasLayer(rlayer_new)
            ##layerList = [layer]
            ##extent = self.iface.canvas.setExtent(rlayer_new.extent())
            ##self.iface.canvas.setLayerSet(layerList)
            ##self.iface.canvas.setVisible(True)         
            #return True
        #else:
            #return False    
        #QMessageBox.information(self, self.tr( "Finished" ), self.tr( "Aquifer media completed." ) )            