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

#******************************************************************************
#
# Freewat
# ---------------------------------------------------------
#
#
# Copyright (C) 2014 - 2015 Iacopo Borsi (iacopo.borsi@tea-group.com)
#
# This source 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.
#
# This code is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
# details.
#
# A copy of the GNU General Public License is available on the World Wide Web
# at <http://www.gnu.org/licenses/>. You can also obtain it by writing
# to the Free Software Foundation, 51 Franklin Street, Suite 500 Boston,
# MA 02110-1335 USA.
#
#******************************************************************************

from PyQt4.QtGui import *
from qgis.core import *
from qgis.gui import *
from qgis.utils import *
from PyQt4.QtCore import *
from PyQt4 import QtGui
from PyQt4.QtCore import QVariant
import ftools_utils as ft
from .logger import Logger
import processing
import locale
import math
from osgeo import ogr, gdal
#import osr
import numpy as np
import sys, os
from collections import OrderedDict



# Initialize Qt resources from file resources.py

# Utility: get a vector layer by name
# mutuated by plugin ftools:
def getVectorLayerNames():
    layerMap = QgsMapLayerRegistry.instance().mapLayers()
    layerNames = []
    for name, layer in layerMap.iteritems():
        if layer.type() == QgsMapLayer.VectorLayer:
            layerNames.append(str(layer.name()))

    # Add Layer group name
    tree = QgsProject.instance().layerTreeRoot().children()

    for node in tree:
        if isinstance(node, QgsLayerTreeGroup):
            layerNames.append(node.name())

    return layerNames

def getGroupLayerByName(groupName):
    groupLayer = QgsProject.instance().layerTreeRoot().findGroup(groupName)

    return groupLayer

def getVectorLayerByName(layerName):
    layerMap = QgsMapLayerRegistry.instance().mapLayers()
    for name, layer in layerMap.iteritems():
        if layer.type() == QgsMapLayer.VectorLayer and layer.name() == layerName:
            if layer.isValid():
                return layer
            else:
                return None

def getVectorSingleLine():
    linelayers = []
    for layer in QgsMapLayerRegistry.instance().mapLayers().values():
        if layer.type() == QgsMapLayer.VectorLayer:
            if layer.wkbType() == 2:
                linelayers.append(str(layer.name()))
    return linelayers

def getVectorPolygon():
    polygonLayers = []
    for layer in QgsMapLayerRegistry.instance().mapLayers().values():
        if layer.type() == QgsMapLayer.VectorLayer:
            if layer.wkbType() == 6 or layer.wkbType() == 4:
                polygonLayers.append(str(layer.name()))
    return polygonLayers


def getVectorSingleFeature():
    linelayers = []
    for layers in QgsMapLayerRegistry.instance().mapLayers().values():
        if layers.type() == QgsMapLayer.VectorLayer:
            if layers.wkbType() == 2:
                if layers.featureCount() == 1:
                    linelayers.append(str(layers.name()))
    return linelayers

def getVectorMultiFeature():
    linelayers = []
    for layers in QgsMapLayerRegistry.instance().mapLayers().values():
        if layers.type() == QgsMapLayer.VectorLayer:
            if layers.wkbType() == 2:
                if layers.featureCount() > 1:
                    linelayers.append(str(layers.name()))
    return linelayers


def getFieldNames(layer ):
    fields = layer.pendingFields()
    fieldNames = []
    for field in fields:
        if not field.name() in fieldNames:
            fieldNames.append(unicode(field.name()))
    return sorted(fieldNames, cmp=locale.strcoll)

def getFieldType(layer, fieldName):
    fields = layer.pendingFields()
    for field in fields:
        if field.name() == fieldName:
            return field.typeName()


def getUniqueValuesCount(layer, fieldIndex, useSelection):
    count = 0
    values = []
    if useSelection:
        for f in layer.selectedFeatures():
            if f[fieldIndex] not in values:
                values.append(f[fieldIndex])
                count += 1
    else:
        request = QgsFeatureRequest().setFlags(QgsFeatureRequest.NoGeometry)
        for f in layer.getFeatures(request):
            if f[fieldIndex] not in values:
                values.append(f[fieldIndex])
                count += 1
    return count

# ---------------------------
def SelectByRegionAndSave_PointOnPolygon(inPoly, inPts):
        # INPUT:
        # inPoly: Multipolygon Layer where you want perform a selection
        # inPts:  Multipoint that intersects the polygon layer
        # OUTPUT:
        # newlayer: the vector layer (multipolygon) where the selection is saved (temporary file)
        inputLayer = ft.getVectorLayerByName(inPoly)
        selectLayer = ft.getVectorLayerByName(inPts)
        inputProvider = inputLayer.dataProvider()
        selectProvider = selectLayer.dataProvider()
        feat = QgsFeature()
        infeat = QgsFeature()
        geom = QgsGeometry()
        selectedSet = []
        index = ft.createIndex(inputProvider)

        selectFit = selectProvider.getFeatures()
        while selectFit.nextFeature(feat):
                geom = QgsGeometry(feat.geometry())
                intersects = index.intersects(geom.boundingBox())
                for id in intersects:
                    inputProvider.getFeatures( QgsFeatureRequest().setFilterFid( int(id) ) ).nextFeature( infeat )
                    tmpGeom = QgsGeometry( infeat.geometry() )
                    if geom.intersects(tmpGeom):
                        selectedSet.append(infeat.id())

        outputSelection = inputLayer.setSelectedFeatures(selectedSet)


        vectorLayer = processing.getObject(inputLayer.name())
        provider = vectorLayer.dataProvider()
        fields = provider.fields()


        newName = 'selectedGrid'

        newlayer = QgsVectorLayer("MultiPolygon", "temporary_selection", "memory")

        pr = newlayer.dataProvider()

        pr.addAttributes(fields.toList())
        newlayer.startEditing()


        features = processing.features(vectorLayer)
        for feat in features:
            newlayer.addFeature(feat)

        return newlayer

def dirDialog(parent):
    settings = QSettings()
    dirName = settings.value("")
    encode = settings.value("")
    fileDialog = QgsEncodingFileDialog(parent, "Browse a directory", dirName, encode)
    fileDialog.setFileMode(QFileDialog.DirectoryOnly)
    fileDialog.setAcceptMode(QFileDialog.AcceptSave)
    fileDialog.setConfirmOverwrite(False)
    if not fileDialog.exec_() == QDialog.Accepted:
            return None, None
    folders = fileDialog.selectedFiles()
    settings.setValue("", QFileInfo(unicode(folders[0])).absolutePath())
    return (unicode(folders[0]), unicode(fileDialog.encoding()))
    print 'Ecco la DIR --------------->>> ', fileDialog.encoding()

def fileDialog(parent):

    myfile = QFileDialog.getOpenFileName(parent, "Choose a file", "", "")
    fileInfo = QFileInfo(myfile)

    return unicode(fileInfo.canonicalFilePath()) # fileName()

def getRasterNames():
    layerMap = QgsMapLayerRegistry.instance().mapLayers()
    layerNames = []
    for name, layer in layerMap.iteritems():
        if layer.type() == QgsMapLayer.RasterLayer:
            layerNames.append(str(layer.name()))
    return layerNames


def getRasterByName(layerName):
    layerMap = QgsMapLayerRegistry.instance().mapLayers()
    for name, layer in layerMap.iteritems():
        if layer.type() == QgsMapLayer.RasterLayer and layer.name() == layerName:
            if layer.isValid():
                return layer
            else:
                return None



def copyRasterToField(inputGrid, inputRaster, outField_name="Z", precise=False):
    # Open data
    gridLayer = QgsMapLayerRegistry.instance().mapLayersByName(inputGrid)[0]
    geometryType = gridLayer.geometryType()

    if geometryType != 2:
        raise Exception("Input vector must be polygon")

    fieldIdx = gridLayer.fieldNameIndex(outField_name)

    if fieldIdx == -1:
        raise Exception("%s field does not exist" % outField_name)

    outField = gridLayer.dataProvider().fields()[fieldIdx]

    if outField.type() != QVariant.Double:
        raise Exception("%s field type must be Double" % outField_name)

    rasterLayer = QgsMapLayerRegistry.instance().mapLayersByName(inputRaster)[0]
    rasterBBox = rasterLayer.extent()
    xOrigin = rasterBBox.xMinimum()
    yOrigin = rasterBBox.yMaximum()
    cellSizeX = rasterLayer.rasterUnitsPerPixelX()
    cellSizeY = rasterLayer.rasterUnitsPerPixelY()
    nCellsXGDAL = rasterLayer.width()
    nCellsYGDAL = rasterLayer.height()

    gridLayer.startEditing()

    for gridCell in gridLayer.getFeatures():
        gridCellBBOX = gridCell.geometry().boundingBox()
        featureBBox = gridCellBBOX.intersect(rasterBBox)
        if featureBBox.isEmpty():
            continue

        offsetX,offsetY,nCellsX,nCellsY = cellInfoForBBox(featureBBox, rasterBBox, cellSizeX, cellSizeY)

        if (( offsetX + nCellsX ) > nCellsXGDAL ):
          nCellsX = nCellsXGDAL - offsetX
        if (( offsetY + nCellsY ) > nCellsYGDAL ):
          nCellsY = nCellsYGDAL - offsetY

        valueSum, count = statisticsFromMiddlePointTest(rasterLayer, gridCell.geometry(), offsetX, offsetY,
                                                        nCellsX, nCellsY, cellSizeX, cellSizeY, precise)
        mean = valueSum/count

        gridLayer.changeAttributeValue(gridCell.id(), fieldIdx, mean)

        #print valueSum,valueSum/count,count

    gridLayer.commitChanges()


def cellInfoForBBox(featureBBox, rasterBBox, cellSizeX, cellSizeY):
    offsetX = (featureBBox.xMinimum() - rasterBBox.xMinimum()) / cellSizeX
    offsetY = (rasterBBox.yMaximum() - featureBBox.yMaximum()) / cellSizeY

    maxColumn = ((featureBBox.xMaximum() - rasterBBox.xMinimum()) / cellSizeX) + 1
    maxRow = ((rasterBBox.yMaximum() - featureBBox.yMinimum()) / cellSizeY) + 1

    nCellsX = maxColumn - offsetX
    nCellsY = maxRow - offsetY

    return int(offsetX), int(offsetY), int(nCellsX), int(nCellsY)


def statisticsFromMiddlePointTest(rasterLayer, poly, pixelOffsetX, pixelOffsetY,
                                  nCellsX, nCellsY, cellSizeX, cellSizeY, precise):
    rasterBBox = rasterLayer.extent()
    polyRectAligned_xmin = rasterBBox.xMinimum() + (pixelOffsetX * cellSizeX)
    polyRectAligned_ymin = rasterBBox.yMaximum() - (pixelOffsetY * cellSizeY) - (nCellsY * cellSizeY)
    polyRectAligned_xmax = rasterBBox.xMinimum() + (pixelOffsetX * cellSizeX) + (nCellsX * cellSizeX)
    polyRectAligned_ymax = rasterBBox.yMaximum() - (pixelOffsetY * cellSizeY)
    polyRectAligned = QgsRectangle(polyRectAligned_xmin, polyRectAligned_ymin,
                                   polyRectAligned_xmax, polyRectAligned_ymax)

    cellsBlock = rasterLayer.dataProvider().block(1, polyRectAligned, nCellsX, nCellsY)
    valuesSum = 0.0
    count = 0

    cellCenterY = rasterBBox.yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2
    for y in range(nCellsY):
        cellCenterX = rasterBBox.xMinimum() + pixelOffsetX * cellSizeX + cellSizeX / 2
        for x in range(nCellsX):
            keep = True
            if precise:
                keep = False
                cellCenter = QgsPoint(cellCenterX, cellCenterY)
                if poly.contains(cellCenter):
                    keep = True
            if keep:
                val = cellsBlock.value(y, x)
                if not math.isnan(val):
                    valuesSum += val
                count += 1
            cellCenterX += cellSizeX
        cellCenterY -= cellSizeY

    return valuesSum, count

def statisticsFromPreciseIntersection():
    pass

# --
def createAttributes(vectorLayer, fields):
    provider = vectorLayer.dataProvider()
    fieldsToCreate = []
    for field in fields:
        if provider.fieldNameIndex(field) != -1:
            continue
        else:
            print "New field %s" % field
            new_field = QgsField(field, QVariant.Double)
            fieldsToCreate.append(new_field)
    provider.addAttributes(fieldsToCreate)
    # hack to update fields
    vectorLayer.startEditing()
    vectorLayer.commitChanges()
# ---
def getModelsInfoLists(layerNameList):
    # From the list of layers loaded in QGis TOC, search for modeltable_ layers
    # and retrieve model names and working path
    # OUTPUT: a tuple of 2 lists

    isok = 0
    modelNameList = []
    pathList = []
    for mName in layerNameList:
        if mName[0:10] == 'modeltable':
            isok = 1
            modelNameTable = getVectorLayerByName(mName)
            for f in modelNameTable.getFeatures():
                nameTemp = f['name']
                pathTemp = f['working_dir']
                modelNameList.append(nameTemp)
                pathList.append(pathTemp)
    ## Message Error if no model is found
    if isok == 0:
        Logger.warning(None, 'No model found!!',
                                    'There is no model table in TOC '
                                    'You have to create a MODEL before '
                                    'running this tool ! ' )
    if isok >= 1:
        return (modelNameList, pathList)
# --
def getModelInfoByName(modelName):
    # from the model name, get relative Working Path and number of stress periods
    layerNameList = getVectorLayerNames()
    (modelNameList, pathList) = getModelsInfoLists(layerNameList)
    i = 0;
    for mName in modelNameList:
        if mName ==  modelName:
           pathfile = pathList[i]
        i = i+1

    timelayer = getVectorLayerByName('timetable_' + modelName)
    ftit = timelayer.getFeatures()
    nsp = 0
    for i in ftit:
        nsp = nsp + 1

    return (pathfile, nsp)
# --
def getModelUnitsByName(modelName):
    # from the model name, get relative Working Path and number of stress periods
    layerNameList = getVectorLayerNames()
    (modelNameList, pathList) = getModelsInfoLists(layerNameList)
    for mName in modelNameList:
        if mName ==  modelName:
           modelNameTable = getVectorLayerByName('modeltable_' + mName )
           for f in modelNameTable.getFeatures():
                    lenuni   = f['length_unit']
                    itmuni = f["time_unit"]

    return (lenuni, itmuni)
# --
def getModelNlayByName(modelName):
    # from the model name, get relative Working Path and number of stress periods
    layerNameList = getVectorLayerNames()
    (modelNameList, pathList) = getModelsInfoLists(layerNameList)

    lpflayer = getVectorLayerByName('lpf_' + modelName)
    ftit = lpflayer.getFeatures()
    nlay = 0
    for i in ftit:
        nlay  += 1

    return nlay
# --
# --
def getTransportModelsByName(modelName):
    # from the model name, get relative Transport Models
    layerNameList = getVectorLayerNames()
    transport = getVectorLayerByName('transport_' + modelName)
    model_list = []

    for ft in transport.getFeatures():
        nameTemp = ft['name']
        model_list.append(nameTemp)
    return model_list

# --
# --
def getModelNspecByName(flowmodelName, transportName):
    # from the model name, get relative number of species, mobile species and Unit of Mass

    trplayer = getVectorLayerByName('species_' + transportName)
    ftit = trplayer.getFeatures()
    # Get n. of species and immobile species, if any
    nspec = 0
    no_mob = 0
    for ft in ftit:
        nspec  += 1
        yes_no_mob = ft['is_mobile']
        if (yes_no_mob == 'no') or (yes_no_mob == 'no'):
            no_mob +=1

    mobile = nspec - no_mob

    trptable = getVectorLayerByName('transport_'+ flowmodelName )

    for f in trptable.getFeatures():
        if f['name'] == transportName:
            mass = f['unit']

    return (nspec, mobile, mass)
# --

def dataFromCSV(csvfile, newname  = 'csv layer', uploadQgis = True):

    # INPUT:
    # csvfile = name of csv file (complete path)
    # newname = name of the new vector layer; default is 'csv layer'
    # uploadQgis = option for uploading the new vector layer in Qgis TOC; default is True

    # OUTPUT:
    # csvl: a vector layer correspondig to CSV file
    # csv_dict: a dictionary, whose KEY is feature iterator (0,1,2, etc.)
    #           and argument are the features value of that record

    # CSV table loader
    uri = QUrl.fromLocalFile(csvfile)
    uri.addQueryItem("type","csv")
    uri.addQueryItem("geomType","none")
    # uri.addQueryItem("delimiter",",")
    # uri.addQueryItem("skipLines","0")
    # uri.addQueryItem("xField","field_1")
    # uri.addQueryItem("yField","field_2")
    csvl = QgsVectorLayer(uri.toString(), newname, "delimitedtext")

    if uploadQgis and csvl.isValid():
		QgsMapLayerRegistry.instance().addMapLayer(csvl)

    # transform csv layer in dictionary
    dp = csvl.dataProvider()
    csv_dict = {}

    i = 0
    for f in csvlayer.getFeatures():
        ft_lst = []
        for fld in dp.fields():
            ft_lst.append(f[fld.name()])
        del fld
        csv_dict[i+1] = ft_lst[1:]
        i += 1

    return (csvl , csv_dict)
# -----
def getUniqueValuesArray(nparray):
    '''
    nparray  = a numpy array
    arrayout = the nmarray with unique values
    '''
    bb, idx  = np.unique(nparray, return_index=True)
    arrayout = bb[np.argsort(idx)]
    return arrayout


def load_table(widget):
    '''
    Converts the item of a QTableWidget table in a QgsVectorLayer loadable in QGIS with proper fields names and value types
    Input: widget has to be a QTableWidget Item!

    WARNING! Method available only from QGIS >= 2.14
    '''
    # create temporary memory layer (this version will work only with QGIS >= 2.14)
    csvl = QgsVectorLayer('Point?crs=EPSG:4326', 'foo' , "memory")
    # data provider
    pr = csvl.dataProvider()

    nfull = 0
    for i in range(0, widget.rowCount()):
        itm = widget.item(i,0)
        if itm is None:
            pass
        else:
            if itm.text() != '':
                nfull += 1


    # create empty dictionary
    dd = OrderedDict()
    # fill the dictionary with keys = column hader names of the Qwidget
    for c in range(widget.columnCount()):
        cn = widget.horizontalHeaderItem(c)
        dd[cn.text()] = [0 for i in range(nfull)]

    # populate the dictionary values withe the QTableWiget attributes
    for i in range(0, nfull):
        itm = widget.item(i,0)
        if itm is None :
            pass
        else:
            for v, k in enumerate(dd.keys()):
                temp = widget.item(i,v)
                dd[k][i] = temp.text()


    # create list of test_dict values that contains the filed type of the table
    def conversion(dd):
        '''
        test for QTableWigetItem types and convert them into QVariant (Int, Double, String)
        Input:
        dd = dictionary; key = field name, values = field values
        '''
        test_dict = {}
        test_dict = dd.fromkeys(dd) # test_dict has the same keys of dd but values are all empty

        for i, k in enumerate(dd.keys()):
            try:
                # test if the values of the original dictionary are INTEGERS
                map(int, dd.values()[i])
                test_dict[k] = QVariant.Int
            except:
                try:
                    # test if the values of the original dictionary are FLOATS
                    map(float, dd.values()[i])
                    test_dict[k] = QVariant.Double
                except:
                    # test if the values of the original dictionary are STRINGS
                    map(str, dd.values()[i])
                    test_dict[k] = QVariant.String
        return test_dict

    test_dict = conversion(dd)
    var = test_dict.values()

    # start editing of the QgsVectorLayer
    csvl.startEditing()

    # populate the Attribute table of the QgsVectorLayer with filed names = dictionary keys and field types = list of type
    for v, k in enumerate(dd.keys()):
        pr.addAttributes([QgsField(k, var[v])])

    # update attribute table
    csvl.updateFields()

    # populate the attribute table with the dictionary values
    ll = []
    for v in range(len(dd.values()[0])):
        ft = QgsFeature()
        # ft.setGeometry(QgsGeometry.fromPoint(QgsPoint(0,0)))
        attr_list = []
        for ik, k in enumerate(dd.keys()):
            value = dd[k][v]
            # Convert from QVariant type to python type
            if var[ik] == QVariant.Int :
                value = int(value)
            elif var[ik] == QVariant.Double :
                value = float(value)
            else:
                value = str(value)
            # Update list of attributes for this feature
            attr_list.append(value)

        ft.setAttributes(attr_list)
        # fake geometry for each point, just as a workaround (see Warning at the beginning of the method). This line is not necessary wih QGIS >= 2.14
        ft.setGeometry(QgsGeometry.fromPoint(QgsPoint(1,1)))
        ll.append(ft)

    # commit the attributes to the provider
    pr.addFeatures(ll)

    # commit the final changes
    csvl.commitChanges()

    # uncomment to load the table in the legend
    # QgsMapLayerRegistry.instance().addMapLayer(csvl)

    return csvl


def read_data(fname, rows, first_block, last_block):
    '''
    creates a dictionary with keys = TS and values = np arrays for each TS of the last rows of each block
    fname = (str) path of the ET_ARRAY.OUT
    rows = (int) number of the model rows
    first_block = (int) first TS, corresponds to diff_ts
    last_block = (int) last TS, corresponds to diff_th
    '''
    with open(fname) as f:
        data = f.readlines()
    blocks = {}
    for block in range(first_block, last_block + 1):
        start = (rows * 2 + 4) * block + (rows * 2 + 4) - rows
        blocks[block] = np.genfromtxt((line[:-1] for line in data[start:start+rows]), autostrip=True)
    return blocks

def deselectAll():
    '''
    just takes all the vector layers loaded in the TOC and deselect all features.
    same thing as pushing the deselect button in the interface
    '''
    for id_, layer in QgsMapLayerRegistry.instance().mapLayers().iteritems():
        if layer.type() == layer.VectorLayer:
            layer.removeSelection()

def ModelPath():
    '''
    returns the complete path of the model
    '''
    # get all the layers in the TOC as dictionary
    layerMap = QgsMapLayerRegistry.instance().mapLayers()
    for k,v in layerMap.items():
        # look for table modeltable
        if 'modeltable' in k:
            vl = v
    mp = os.path.dirname(vl.dataProvider().dataSourceUri())
    mp = mp + '/'
    mp = mp[8:]
    return mp


def load_obh(filename):
    '''
    filename = path of the .obh file
    outputs:
        unique_sp = list of uniques SP
        total_sp = list of all the SP in the file loaded
    '''

    f = open(filename, 'r')

    # skip the first 3 useless lines
    line = f.readlines()[3:]

    total_sp = []

    for i in line:
        t = i.strip().split()
        try:
            total_sp.append(int(t[4]))
        except IndexError:
            pass

    unique_sp = []

    for itm in total_sp:
        if itm not in unique_sp:
            unique_sp.append(itm)

    unique_sp.sort()

    return unique_sp, total_sp

    f.close()


def pop_message(text, value):
    '''
    simple methos to call to pop up a message window

    text: String. Text message that will be displayed in the window
    value: String. Type of window that will be opened, available values are:
        1. information for and information box
        2. warning for a warning box
        3. critical for a sever error
    '''
    if value == 'information':
        Logger.information(None, 'Information', text)
    elif value == 'warning':
        Logger.warning(None, 'Warning', text)

    elif value == 'critical':
        Logger.critical(None, 'Critical', text)



class ComboStyledItemDelegate(QStyledItemDelegate):
    """this class is necessary for OS X and linux to dosplay comboboxes of
    checkboxes correctly"""

    def __init__(self, parent=None):
        super(ComboStyledItemDelegate, self).__init__(parent)

    def paint(self, painter_, option_, index_):
        option_.showDecorationSelected = False
        QStyledItemDelegate.paint(self, painter_, option_, index_);

