# -*- 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.QtCore import *
from PyQt4.QtGui import *
from qgis.core import *
from qgis.gui import *
from qgis.utils import *
from pyspatialite import dbapi2 as sqlite3
import processing
from freewat_utils import getVectorLayerByName, createAttributes, getVectorLayerNames, getModelsInfoLists, deselectAll
from freewat_formulas import cond, lerp, timeoffsetCalc, timeoffsetCalcIsodate
from sqlite_utils import uploadQgisVectorLayer, checkIfTableExists
import sys, os
import math
import datetime

# ---------------------------
# Create a new MODEL: its DB SQLite, along with its modeltable and its timetable
def createModel(modelname, workingDir, isChild, lengthString, timeString, timeParameters, initDateString, initTimeString,  crs):

    QApplication.processEvents()
    # Set file name

    dbname = workingDir + '/' + modelname + ".sqlite"
    # creating/connecting SQL database object
    con = sqlite3.connect(dbname)
    # con = sqlite3.connect(":memory:") if you want write it in RAM
    con.enable_load_extension(True)
    cur = con.cursor()
    # Initialise spatial db
    # Insert a control for existing table or not
    cur.execute("PRAGMA foreign_keys = ON")
    # cur.execute("SELECT load_extension('libspatialite.so.5');")
    # Initializing Spatial MetaData: GEOMETRY_COLUMNS and SPATIAL_REF_SYS
    cur.execute("SELECT InitSpatialMetaData(1);")


    # Create new model table
    # Create new SQL table

    modeltableName = 'modeltable_' + modelname

    SQLstring = 'CREATE TABLE "%s" ("name" varchar(20), "length_unit" varchar(20), "time_unit" varchar(20), "is_child" integer, "working_dir" varchar(99), "initial_date" varchar(99), "initial_time" varchar(99), "crs" varchar(99))' %modeltableName

    cur.execute(SQLstring)

    #comment line for no more geometry column
    #cur.execute("SELECT AddGeometryColumn('%s', 'the_geom', 4326, 'POINT', 'XY');"%modeltableName)

    # Insert values in Model Table
    parameters = [modelname, lengthString, timeString, isChild, workingDir, initDateString, initTimeString, crs]
    sqlstr = 'INSERT INTO %s'%modeltableName
    cur.execute(sqlstr + ' (name,length_unit,time_unit,is_child,working_dir,initial_date, initial_time, crs) VALUES (?, ?, ?, ?, ?,?, ?, ?);', (parameters[0], parameters[1], parameters[2], parameters[3], parameters[4], parameters[5], parameters[6], parameters[7]))

    # Create and initialize the model time table
    # Create new SQL table
    timeName = 'timetable_%s'%modelname

    sql2 = 'CREATE TABLE "%s" '%timeName
    sql22 = sql2 + ' ("ID" INTEGER NOT NULL PRIMARY KEY AUTOINCREMENT, "sp" integer,"length" float, "ts" integer, "multiplier" float, "state" varchar(4) );'

    cur.execute(sql22)

    #comment line for no more geometry column
    #sql3 = " SELECT AddGeometryColumn( '%s', 'the_geom', 4326, 'POINT', 'XY'); "%timeName
    #cur.execute(sql3)

    # Initialize values in Model Table

    sql4 = 'INSERT INTO %s'%timeName
    sql44 = sql4 + '   ( sp, length, ts, multiplier, state) VALUES ( ?, ?, ?, ?, ?);'

    cur.execute(sql44, (1, timeParameters[0], timeParameters[1], timeParameters[2], timeParameters[3] ))
    # Close cursor
    cur.close()
    # Save the changes
    con.commit()
    # Close connection
    con.close()

    # Add the model table into QGis map
    uri = QgsDataSourceURI()
    uri.setDatabase(dbname)
    schema = ''
    table = modeltableName
    #geom_column = 'the_geom'
    geom_column = None
    uri.setDataSource(schema, table, geom_column)
    display_name = modeltableName
    tableLayer = QgsVectorLayer(uri.uri(), display_name, 'spatialite')

    QgsMapLayerRegistry.instance().addMapLayer(tableLayer)


     # Add the model time table into QGis map
    uri2 = QgsDataSourceURI()
    uri2.setDatabase(dbname)
    uri2.setDataSource(schema, timeName, geom_column)
    display_name2 = timeName
    timeLayer = QgsVectorLayer(uri2.uri(), display_name2, 'spatialite')

    QgsMapLayerRegistry.instance().addMapLayer(timeLayer)

# Create a new MDO layer
def createMDO(gridLayer, dbName, newName, fieldNameList, fieldTypeList, fieldDefaultList):

    provider = gridLayer.dataProvider() #provider = vectorLayer.dataProvider()

    #retreive information of the modeltable and the crs field
    layerNameList = getVectorLayerNames()
    (modelNameList, pathList) =  getModelsInfoLists(layerNameList)

    #match the exact model trough the list of all models
    db = [os.path.basename(dbName).replace('.sqlite', '')]
    for i in db:
        if i in modelNameList:
            modelName = 'modeltable_' + str(i)


    model = getVectorLayerByName(modelName)
    for ft in model.getFeatures():
        crs = ft['crs']

    #vl = QgsVectorLayer("Polygon?crs=epsg:4326", "temporary_points", "memory")
    vl = QgsVectorLayer("Polygon?crs=" + crs, "temporary_points", "memory")
    pr = vl.dataProvider()
    fld = gridLayer.dataProvider().fields()

    idfield = 0
    for f in fld:
    	pr.addAttributes([f])

    for feature in processing.features(gridLayer):
        att = feature.attributes()
        pr.addFeatures([feature])
    vl.commitChanges()
    vl.startEditing()

    # Add new fields
    for i in range(0, len(fieldNameList)):
        newfield = QgsField(fieldNameList[i], fieldTypeList[i])
        res = pr.addAttributes( [newfield] )
    vl.updateFields()

    # Add default values for new fields
    fields = vl.dataProvider().fields()
    idval = 0
    for field in fieldNameList:
        idx = fields.indexFromName(field)#retrieve field index from name
        for f in vl.getFeatures():
            attrs = f.attributes()
            for attr in attrs:
                feat = f.id()
                vl.changeAttributeValue(feat, idx, fieldDefaultList[idval])
                QApplication.processEvents()
        idval += 1

    vl.commitChanges()


    # Upload the vlayer into DB SQlite
    uploadQgisVectorLayer(dbName, vl, newName)

    # Retrieve the Spatialite layer and add it to mapp
    uri = QgsDataSourceURI()
    uri.setDatabase(dbName)
    schema = ''
    table = newName
    geom_column = "Geometry"
    uri.setDataSource(schema, table, geom_column)
    display_name = table

    wlayer = QgsVectorLayer(uri.uri(), display_name, 'spatialite')

    QgsMapLayerRegistry.instance().addMapLayer(wlayer)

    deselectAll()


# Create a new MDO layer applied to the whole grid (faster!):
def createMDOLarge(gridLayer, dbName, newName, fieldNameList, fieldTypeList, fieldDefaultList):

    provider = gridLayer.dataProvider() #provider = vectorLayer.dataProvider()

    #retreive information of the modeltable and the crs field
    layerNameList = getVectorLayerNames()
    (modelNameList, pathList) =  getModelsInfoLists(layerNameList)

    #match the exact model trough the list of all models
    db = [os.path.basename(dbName).replace('.sqlite', '')]
    for i in db:
        if i in modelNameList:
            modelName = 'modeltable_' + str(i)

    model = getVectorLayerByName(modelName)
    for ft in model.getFeatures():
        crs = ft['crs']

    #vl = QgsVectorLayer("Polygon?crs=epsg:4326", "temporary_points", "memory")
    vl = QgsVectorLayer("Polygon?crs=" + crs, "temporary_points", "memory")
    pr = vl.dataProvider()
    fld = gridLayer.dataProvider().fields()

    # Add the original grid fields
    idfield = 0
    for f in fld:
    	pr.addAttributes([f])

    # Add new fields
    for i in range(0, len(fieldNameList)):
        newfield = QgsField(fieldNameList[i], fieldTypeList[i])
        res = pr.addAttributes( [newfield] )
    vl.updateFields()

    # Now, put the layer into DB:
    # Upload the vlayer into DB SQlite
    uploadQgisVectorLayer(dbName, vl, newName)

    ## --- Inlcude new features, using default values
    # build the sql string : first the GRID field names
    sqlinsertField = """INSERT INTO %s ("""%newName
    sqlinsertField += "Geometry, "
    for f in fld:
        sqlinsertField += "%s , "%f.name()

    # then the specific fields (here we have to manage with the comma and ) at the end !
    for j, nf in enumerate(fieldNameList):
        if j < ( len(fieldDefaultList) - 1):
            sqlinsertField += "%s , "%nf
        else:
            sqlinsertField += "%s ) "%nf

    # create the connection SQL database object
    con = sqlite3.connect(dbName)
    cur = con.cursor()

    # complete the sql string with values
    # get values from grid table
    rs = cur.execute('SELECT * FROM %s'%gridLayer.name())
    grrows = cur.fetchall()
    geom_rs = cur.execute('SELECT AsText(geometry) FROM %s'%gridLayer.name())
    geomrows = cur.fetchall()

    for j,g in enumerate(grrows):

        sqlinsert = sqlinsertField
        crsOnlyCode = str(crs[5:])
        sqlinsert += " VALUES ( MultipolygonFromText( ' %s ', %s) , "%(geomrows[j][0], crsOnlyCode)
        attrs = []
        attrs.extend([grrows[j][0],0, grrows[j][3], grrows[j][4]])
        attrs.extend(fieldDefaultList)
        ftins = tuple(attrs)

        nf = len(vl.fields()) -1
        for i,ft in enumerate(attrs):
            # Here, again, we have to specify string in case of the last feature
            if i == nf :
                sqlinsert += "?);"
            # ... for all the others :
            else:
                sqlinsert += "?,"

        cur.execute(sqlinsert, ftins)

    # Close SpatiaLiteDB
    cur.close()
    # Save the changes
    con.commit()
    # Close connection
    con.close()

    # Retrieve the Spatialite layer and add it to mapp
    uri = QgsDataSourceURI()
    uri.setDatabase(dbName)
    schema = ''
    table = newName
    geom_column = "Geometry"
    uri.setDataSource(schema, table, geom_column)
    display_name = table

    # Get the DB layer as a Qgs layer object
    wlayer = QgsVectorLayer(uri.uri(), display_name, 'spatialite')

    QgsMapLayerRegistry.instance().addMapLayer(wlayer)

    deselectAll()

##


def createModelLayer(gridLayer, pathFile, modelname, layerName, parameters):

    name_fields = ['BORDER' ,'ACTIVE','TOP', 'BOTTOM', 'THICKNESS', 'STRT', 'KX', 'KY', 'KZ', 'SS', 'SY', 'NT', 'NE', 'WETDRY']
    type_fields = [QVariant.Int , QVariant.Int, QVariant.Double,  QVariant.Double , QVariant.Double, QVariant.Double, QVariant.Double, QVariant.Double, QVariant.Double, QVariant.Double, QVariant.Double, QVariant.Double, QVariant.Double, QVariant.Double ]

    dbName = pathFile + '/' + modelname + '.sqlite'

    # Convert info into MODFLOW flags:
    layType    = parameters[1]
    layAverage  = parameters[2]
    chani   = parameters[3]
    layWet  = parameters[4]
    top = parameters[5]
    bottom = parameters[6]

    default_fields = [0, 1, top, bottom, (top - bottom ), 1.0, 0.001, 0.001, 0.0001, 0.001, 0.1, 1, 1, -0.01]

    createMDOLarge(gridLayer,dbName,layerName, name_fields, type_fields, default_fields)



    # Create or update LPF table of the Model
    nameTable = "lpf_"+ modelname



    # Define values of new feature for LPF
    # -- Create list ft
    ft = [layerName,layType, layAverage,chani,layWet]
    # --
    #Convert in a tuple
    ftinsert = tuple(ft)

    # Check if LPF table exists:
    lpflayer = getVectorLayerByName(nameTable)


    if lpflayer == None :
        # Connect to DB
        # creating/connecting SQL database object
        con = sqlite3.connect(dbName)
        # con = sqlite3.connect(":memory:") if you want write it in RAM
        con.enable_load_extension(True)
        cur = con.cursor()
        # Create new LPF table in DB
##        SQLcreate = 'CREATE TABLE %s ("name" varchar(20), "type" integer,' \
##                '"layavg" integer, "chani" integer, "laywet" integer );'% nameTable
        SQLcreate = 'CREATE TABLE %s ("name" varchar(20), "type" varchar(20),' \
                '"layavg" varchar(20), "chani" integer, "laywet" varchar(20) );'% nameTable

        cur.execute(SQLcreate)
        #comment line for no more geometry column
        #sql2 = "SELECT AddGeometryColumn( '%s', 'Geometry', 4326, 'POINT', 'XY');"% nameTable
        #cur.execute(sql2)

        # Insert values in Table
        sql3 = 'INSERT INTO %s '%nameTable +  '(name,type,layavg,chani,laywet)'

        cur.execute(sql3 + 'VALUES (?, ?, ?, ?, ?);', ftinsert)
        # Close SpatiaLiteDB
        cur.close()
        # Save the changes
        con.commit()
        # Close connection
        con.close()


        # Retrieve the Spatialite layer and add it to mapp
        uri = QgsDataSourceURI()
        uri.setDatabase(dbName)
        schema = ''
        table = nameTable
        #geom_column = "Geometry"
        geom_column = None
        uri.setDataSource(schema, table, geom_column)
        display_name = nameTable

        wlayer = QgsVectorLayer(uri.uri(), display_name, 'spatialite')


        # set the Widget of the column as ValusMap, that is as combobox with selectable values
        wlayer.setEditorWidgetV2(1, 'ValueMap')
        wlayer.setEditorWidgetV2Config(1, {'confined':'confined', 'convertible':'convertible'})

        wlayer.setEditorWidgetV2(2, 'ValueMap')
        wlayer.setEditorWidgetV2Config(2, {'harmonic':'harmonic', 'logarithmic':'logarithmic', 'arithmetic-mean':'arithmetic-mean'})

        wlayer.setEditorWidgetV2(4, 'ValueMap')
        wlayer.setEditorWidgetV2Config(4, {'No':'No', 'Yes':'Yes'})



        QgsMapLayerRegistry.instance().addMapLayer(wlayer)

    else:
        # Update existing LPF table
        # Add features related to the new model layer
        dp = lpflayer.dataProvider()
        fields = lpflayer.pendingFields()
        ftnew = QgsFeature(fields)
        for i in range(0,len(ft)):
            ftnew[i] = ft[i]

        dp.addFeatures( [ ftnew ] )

# --- Create a CHD Layer
def createChdLayer(gridLayer, name, pathFile, modelName, nsp):

    QApplication.processEvents()

    layerName = name + "_chd"
    dbName = pathFile + '/' + modelName + '.sqlite'

    name_fields = []
    type_fields = []

    name_fields.append('from_lay')
    type_fields.append(QVariant.Int)
    name_fields.append('to_lay')
    type_fields.append(QVariant.Int)

    # Initialize Default values for fields:
    default_fields = [1 , 1]

    for i in range(1,nsp+1):
        name_fields.append(str(i) + '_shead')
        type_fields.append(QVariant.Double)
        default_fields.append(10)

        name_fields.append(str(i) + '_ehead')
        type_fields.append(QVariant.Double)
        default_fields.append(10)

    # Create MDO for CHD:
    createMDO(gridLayer, dbName, layerName, name_fields, type_fields, default_fields)

# ---------------------------
# --- Create a Well Layer
def createWellLayer(gridLayer, name, pathFile, modelName, nsp):

    QApplication.processEvents()

    layerName = name + "_well"
    dbName = pathFile + '/' + modelName + '.sqlite'

    name_fields = []
    type_fields = []

    name_fields.append('from_lay')
    type_fields.append(QVariant.Int)
    name_fields.append('to_lay')
    type_fields.append(QVariant.Int)
    name_fields.append('active')
    type_fields.append(QVariant.Int)
    name_fields.append('group')
    type_fields.append(QVariant.Int)
    name_fields.append('use')
    type_fields.append(QVariant.String)

    default_fields = [1, 1, 1, 1, '']

    for i in range(1,nsp+1):
        name_fields.append('sp_'+ str(i) )
        type_fields.append(QVariant.Double)
        default_fields.append(-100)


    # Create MDO for WEL:
    createMDO(gridLayer, dbName, layerName, name_fields, type_fields, default_fields)

# ---------------------------
# --- Create a MNW Layer
def createMnwLayer(gridLayer, name, pathFile, modelName, nsp):

    layerName = name + "_mnw"
    dbName = pathFile + '/' + modelName + '.sqlite'

    name_fields = []
    type_fields = []

    name_fields.append('WELLID')
    type_fields.append(QVariant.String)
    # name_fields.append('to_lay')
    # type_fields.append(QVariant.Int)
    # name_fields.append('active')
    # type_fields.append(QVariant.Int)

    default_fields =  ['well_1']  #,1,1]

    for i in range(1,nsp+1):
        name_fields.append('Qw_'+ str(i) )
        type_fields.append(QVariant.Double)
        default_fields.append(-100)


    # Create MDO for WEL:
    createMDO(gridLayer, dbName, layerName, name_fields, type_fields, default_fields)


# ---------------------------
# --- Create a RCH Layer
def createRchLayer(gridLayer, name, pathFile, modelName, nsp):

    QApplication.processEvents()

    layerName = name + "_rch"
    dbName = pathFile + '/' + modelName + '.sqlite'
    #
    name_fields = []
    type_fields = []
    default_fields = []

    for i in range(1,nsp+1):
        name_fields.append( 'sp_' + str(i) + '_rech')
        type_fields.append(QVariant.Double)
        default_fields.append(0.01)

        name_fields.append('sp_' + str(i) + '_irch')
        type_fields.append(QVariant.Int)
        default_fields.append(1)

    createMDOLarge(gridLayer, dbName, layerName, name_fields, type_fields, default_fields)
# ---------------------------

def createChdLayerFromLine(newName, dbName, gridLayer, riverLayer, csvlayer, nsp, from_lay, to_lay):

    if riverLayer.wkbType() != 2:
        raise Exception("Input river layer must by single part linestring")

    if riverLayer.featureCount() > 1:
        raise Exception("Input river layer must contain only one feature")

    # trasformo csv layer in dictionary
    dp = csvlayer.dataProvider()
    csv_dict = {}

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

    # convert the values of the dictionary in usable float numbers if they contain ',' and are string of unicode
    for k in csv_dict:
        for i in range(len(csv_dict[k])):
            if type(csv_dict[k][i]) == str or type(csv_dict[k][i]) == unicode and ',' in csv_dict[k][i]:
                csv_dict[k][i] = float(csv_dict[k][i].replace(',', '.'))
                QApplication.processEvents()


    # retrieve information of the modeltable and the crs field
    layerNameList = getVectorLayerNames()
    (modelNameList, pathList) =  getModelsInfoLists(layerNameList)

    #match the exact model trough the list of all models
    db = [os.path.basename(dbName).replace('.sqlite', '')]
    for i in db:
        if i in modelNameList:
            modelName = 'modeltable_' + str(i)

    model = getVectorLayerByName(modelName)
    for ft in model.getFeatures():
        crs = ft['crs']

    provider = gridLayer.dataProvider()
    grid = QgsVectorLayer("Polygon?crs=" + crs, "temporary_layer", "memory")
    pr = grid.dataProvider()

    # fields of original grid layer
    fld = gridLayer.dataProvider().fields()
    for f in fld:
    	pr.addAttributes([f])

    # get the selection of grid cells that intersect the river layer
    for i in gridLayer.getFeatures():
        for l in riverLayer.getFeatures():
            if i.geometry().intersects(l.geometry()):
                gridLayer.select(i.id())

    # create a request with the selection
    request = QgsFeatureRequest().setFilterFids(gridLayer.selectedFeaturesIds())

    # features from the original grid layer with the request
    for feature in gridLayer.getFeatures(request):
        pr.addFeatures([feature])

    grid.startEditing()
    grid.commitChanges()

    # Add the new fields 'layer' and 'segment' (Integer )
    newfield = QgsField('from_lay', QVariant.Int)
    pr.addAttributes([newfield])
    newfield = QgsField('to_lay', QVariant.Int)
    pr.addAttributes([newfield])

    # Define List of new fields name
    fieldsList = []
    for i in range(1, nsp + 1):
        fieldsList.append(str(i) + '_shead')
        fieldsList.append( str(i) + '_ehead')


    # Call the method to write new fields to be added
    createAttributes(grid, fieldsList)

    #  Start editing to write record
    grid.startEditing()

    changedAttributesDict = {}
    ft_sel_id = []

    # pass the request to getFeatures (speeds up the process)
    for feature in grid.getFeatures():
    # for feature in grid.getFeatures(request):
        QApplication.processEvents()
        gridPoly = feature.geometry()
        inCellRiverLength = 0.0
        inCoeff = 0.0
        outCoeff = 0.0
        # I take the only one river
        riverFeature = [f for f in riverLayer.getFeatures(QgsFeatureRequest(0))][0]
        riverLine = riverFeature.geometry()
        riverLength = riverLine.length()
        riverPoints = riverLine.asPolyline()
        #
        if riverLine.intersects(gridPoly):
            ft_sel_id.append(feature.id())
            inCellRiver = riverLine.intersection(gridPoly)
            if inCellRiver.wkbType() == QGis.WKBLineString:
                linepoints = inCellRiver.asPolyline()
                # first point
                pointIn = linepoints[0]
                # last point
                pointOut = linepoints[len(linepoints)-1]
            elif inCellRiver.wkbType() == QGis.WKBMultiLineString:
                multilinepoints = inCellRiver.asMultiPolyline()
                # first point of first list
                pointIn = multilinepoints[0][0]
                lastGroup = multilinepoints[len(multilinepoints)-1]
                # last point of last list
                pointOut = lastGroup[len(lastGroup)-1]

            _, _, beforeVertexIndexpointIn, _, _ = riverLine.closestVertex(QgsPoint(pointIn))
            _, _, beforeVertexIndexpointOut, _, _ = riverLine.closestVertex(QgsPoint(pointOut))
            beforeVertexpointIn = riverPoints[beforeVertexIndexpointIn]
            beforeVertexpointOut = riverPoints[beforeVertexIndexpointOut]
            if beforeVertexIndexpointIn == -1:
                beforeVertexpointIn = pointIn
            if beforeVertexIndexpointOut == -1:
                beforeVertexpointOut = pointOut

            if beforeVertexIndexpointIn == -1:
                beforeVertexIndexpointIn = 0
            if beforeVertexIndexpointOut == -1:
                beforeVertexIndexpointOut = 0

            stretchIn = riverPoints[0: beforeVertexIndexpointIn+1]
            stretchIn.append(pointIn)
            stretchInGeom = QgsGeometry.fromPolyline(stretchIn)
            currentInLength = stretchInGeom.length()
            stretchOut = riverPoints[0: beforeVertexIndexpointOut+1]
            stretchOut.append(pointOut)

            stretchOutGeom = QgsGeometry.fromPolyline(stretchOut)
            currentOutLength = stretchOutGeom.length()
            inCoeff = currentInLength / riverLength
            outCoeff = currentOutLength / riverLength
            changedAttributes = {}

            # Insert value to from_lay field
            fieldIdx = grid.dataProvider().fieldNameIndex('from_lay')
            changedAttributes[fieldIdx] = from_lay

            # Insert value to to_lay field
            fieldIdx = grid.dataProvider().fieldNameIndex('to_lay')
            changedAttributes[fieldIdx] = to_lay

            for sp, field in csv_dict.iteritems():
                QApplication.processEvents()
                # shead_in, shead_out, ehead_in, ehead_out
                # river head
                valueIn = field[0]
                valueOut = field[1]
                valueIn_int = lerp(valueIn, valueOut, inCoeff)
                valueOut_int = lerp(valueIn, valueOut, outCoeff)
                avg = (valueIn_int + valueOut_int)/2

                fieldIdx = grid.dataProvider().fieldNameIndex('{}_shead'.format(sp))
                changedAttributes[fieldIdx] = avg

                # river bottom
                valueIn = field[2]
                valueOut = field[3]
                valueIn_int = lerp(valueIn, valueOut, inCoeff)
                valueOut_int = lerp(valueIn, valueOut, outCoeff)
                avg = (valueIn_int + valueOut_int)/2

                fieldIdx = grid.dataProvider().fieldNameIndex('{}_ehead'.format(sp))
                changedAttributes[fieldIdx] = avg

            changedAttributesDict[feature.id()] = changedAttributes

    grid.dataProvider().changeAttributeValues(changedAttributesDict)
    grid.commitChanges()

    # Select only intersected cells (whose IDs are in ft_sel_id)
    grid.select(ft_sel_id)


    # Upload the vlayer into DB SQlite
    newName = newName + '_chd'
    uploadQgisVectorLayer(dbName, grid, newName, srid= grid.crs().postgisSrid(), selected = True)

    # Retrieve the Spatialite layer and add it to mapp
    uri = QgsDataSourceURI()
    uri.setDatabase(dbName)
    schema = ''
    table = newName
    geom_column = "Geometry"
    uri.setDataSource(schema, table, geom_column)
    display_name = table

    wlayer = QgsVectorLayer(uri.uri(), display_name, 'spatialite')

    QgsMapLayerRegistry.instance().addMapLayer(wlayer)

    deselectAll()


def createRivLayer(newName, dbName, gridLayer, riverLayer, csvlayer, width, layer, xyz, nsp ):
    # --- Input:
    # newName: name of the new MDO for river
    # dbName: name of model DB (complete path)
    # gridLayer: grid layer
    # csvlayer: csv table with parameters (for each stress period)
    #           [sp, rh_in, rh_out, bt_in, bt_out, krb_in, krb_out, thk_in, thk_out]
    # river: line layer
    # width: river width, given as scalar (constant)
    # layer: layer where RIV is applied  (integer)
    # xyz: river segment id (integer)
    # nsp: stress period number

    # Check if river layer has the correct format
    # To do: change check of Wkb into Geometry type



    if riverLayer.wkbType() != 2:
        raise Exception("Input river layer must by single part linestring")

    if riverLayer.featureCount() > 1:
        raise Exception("Input river layer must contain only one feature")

    ##
    # trasformo csv layer in dictionary
    dp = csvlayer.dataProvider()
    csv_dict = {}

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

    # convert the values of the dictionary in usable float numbers if they contain ',' and are string of unicode
    for k in csv_dict:
        for i in range(len(csv_dict[k])):
            if type(csv_dict[k][i]) == str or type(csv_dict[k][i]) == unicode and ',' in csv_dict[k][i]:
                csv_dict[k][i] = float(csv_dict[k][i].replace(',', '.'))
                QApplication.processEvents()


    #retreive information of the modeltable and the crs field

    layerNameList = getVectorLayerNames()
    (modelNameList, pathList) =  getModelsInfoLists(layerNameList)

    #match the exact model trough the list of all models
    db = [os.path.basename(dbName).replace('.sqlite', '')]
    for i in db:
        if i in modelNameList:
            modelName = 'modeltable_' + str(i)

    model = getVectorLayerByName(modelName)
    for ft in model.getFeatures():
        crs = ft['crs']


    #
    provider = gridLayer.dataProvider()
    grid = QgsVectorLayer("Polygon?crs=" + crs, "temporary_layer", "memory")
    pr = grid.dataProvider()

    # fields of original grid layer
    fld = gridLayer.dataProvider().fields()
    for f in fld:
    	pr.addAttributes([f])


    # get the selection of grid cells that intersect the river layer
    for i in gridLayer.getFeatures():
        for l in riverLayer.getFeatures():
            if i.geometry().intersects(l.geometry()):
                gridLayer.select(i.id())


    # create a request with the selection
    request = QgsFeatureRequest().setFilterFids(gridLayer.selectedFeaturesIds())


    # features from the original grid layer
    for feature in gridLayer.getFeatures(request):
        pr.addFeatures([feature])

    grid.startEditing()
    grid.commitChanges()

    # Add the new fields 'layer' and 'segment' (Integer )
    newfield = QgsField('layer', QVariant.Int)
    pr.addAttributes( [newfield] )
    newfield = QgsField('segment', QVariant.Int)
    pr.addAttributes( [newfield] )


    # Define List of new fields name
    fieldsList = ['length']
    for i in range(1,nsp+1):
        fieldsList.append('stage_'+ str(i) )
        fieldsList.append('rbot_'+ str(i) )
        fieldsList.append('cond_'+ str(i) )


    # Call the method to write new fields to be added
    createAttributes(grid, fieldsList)

    #  Start editing to write record
    grid.startEditing()


    changedAttributesDict = {}
    #
    ft_sel_id = []

    # pass the request to getFeatures (speeds up the process)
    for feature in grid.getFeatures():
    # for feature in grid.getFeatures(request):
        QApplication.processEvents()
        gridPoly = feature.geometry()
        inCellRiverLength = 0.0
        inCoeff = 0.0
        outCoeff = 0.0
        # I take the only one river
        riverFeature = [f for f in riverLayer.getFeatures(QgsFeatureRequest(0))][0]
        riverLine = riverFeature.geometry()
        riverLength = riverLine.length()
        riverPoints = riverLine.asPolyline()
        #
        if riverLine.intersects(gridPoly):
            #
            ft_sel_id.append(feature.id())
            #
            inCellRiver = riverLine.intersection(gridPoly)
            # print "Grid ID %s - River %s - Lenght in cell: %s" % (feature.id(),riverFeature.id(),inCellRiver.length())
            if inCellRiver.wkbType() == QGis.WKBLineString:
                linepoints = inCellRiver.asPolyline()
                # first point
                pointIn = linepoints[0]
                # last point
                pointOut = linepoints[len(linepoints)-1]
                #print "Ingresso: %s %s" % (pointIn.x(),pointIn.y())
                #print "Uscita: %s %s" % (pointOut.x(),pointOut.y())
            elif inCellRiver.wkbType() == QGis.WKBMultiLineString:
                multilinepoints = inCellRiver.asMultiPolyline()
                # first point of first list
                pointIn = multilinepoints[0][0]
                lastGroup = multilinepoints[len(multilinepoints)-1]
                # last point of last list
                pointOut = lastGroup[len(lastGroup)-1]
                #print "Ingresso: %s %s" % (pointIn.x(),pointIn.y())
                #print "Uscita: %s %s" % (pointOut.x(),pointOut.y())

            _, _, beforeVertexIndexpointIn, _, _ = riverLine.closestVertex(QgsPoint(pointIn))
            _, _, beforeVertexIndexpointOut, _, _ = riverLine.closestVertex(QgsPoint(pointOut))
            beforeVertexpointIn = riverPoints[beforeVertexIndexpointIn]
            beforeVertexpointOut = riverPoints[beforeVertexIndexpointOut]
            #print "Indice vertice precedente Ingresso: %s" % (beforeVertexIndexpointIn)
            #print "Indice vertice precedente Uscita: %s" % (beforeVertexIndexpointOut)
            if beforeVertexIndexpointIn == -1:
                beforeVertexpointIn = pointIn
            if beforeVertexIndexpointOut == -1:
                beforeVertexpointOut = pointOut
            #print "Vertice precedente Ingresso: %s %s" % (beforeVertexpointIn.x(),beforeVertexpointIn.y())
            #print "Vertice precedente Uscita: %s %s" % (beforeVertexpointOut.x(),beforeVertexpointOut.y())

            if beforeVertexIndexpointIn == -1:
                beforeVertexIndexpointIn = 0
            if beforeVertexIndexpointOut == -1:
                beforeVertexIndexpointOut = 0

            stretchIn = riverPoints[0: beforeVertexIndexpointIn+1]
            stretchIn.append(pointIn)
            stretchInGeom = QgsGeometry.fromPolyline(stretchIn)
            currentInLength = stretchInGeom.length()
            stretchOut = riverPoints[0: beforeVertexIndexpointOut+1]
            stretchOut.append(pointOut)

##            try:
            stretchOutGeom = QgsGeometry.fromPolyline(stretchOut)
            currentOutLength = stretchOutGeom.length()
            inCoeff = currentInLength / riverLength
            outCoeff = currentOutLength / riverLength
##                print "Parziale ingresso: %s" % (currentInLength)
##                print "Parziale uscita: %s" % (currentOutLength)
##                print "Coefficiente interpolazioni ingresso: %s" % (float(currentInLength)/riverLine.length())
##                print "Coefficiente interpolazioni ingresso: %s" % (float(currentOutLength)/riverLine.length())
##                print "Grid ID: %s" % (feature.id())
            changedAttributes = {}
            # Insert value to layer field (layer)
            fieldIdx = grid.dataProvider().fieldNameIndex('layer')
            changedAttributes[fieldIdx] = layer

            # Insert value to segment field (segment)
            fieldIdx = grid.dataProvider().fieldNameIndex('segment')
            changedAttributes[fieldIdx] = xyz

            # Insert value to length field
            fieldIdx = grid.dataProvider().fieldNameIndex('length')
            changedAttributes[fieldIdx] = inCellRiver.length()
            Lreach = inCellRiver.length()

            for sp, field in csv_dict.iteritems():
                QApplication.processEvents()
                # rh_in, rh_out, bt_in, bt_out, krb_in, krb_out, thk_in, thk_out
                # river head
                valueIn = field[0]
                valueOut = field[1]
                valueIn_int = lerp(valueIn, valueOut, inCoeff)
                valueOut_int = lerp(valueIn, valueOut, outCoeff)
                avg = (valueIn_int + valueOut_int)/2

                fieldIdx = grid.dataProvider().fieldNameIndex('stage_%i' % sp)
                changedAttributes[fieldIdx] = avg

                # river bottom
                valueIn = field[2]
                valueOut = field[3]
                valueIn_int = lerp(valueIn, valueOut, inCoeff)
                valueOut_int = lerp(valueIn, valueOut, outCoeff)
                avg = (valueIn_int + valueOut_int)/2

                fieldIdx = grid.dataProvider().fieldNameIndex('rbot_%i' % sp)
                changedAttributes[fieldIdx] = avg

                # river thickness
                valueIn = field[6]
                valueOut = field[7]
                valueIn_int = lerp(valueIn, valueOut, inCoeff)
                valueOut_int = lerp(valueIn, valueOut, outCoeff)
                th_avg = (valueIn_int + valueOut_int)/2

                # river conductance
                valueIn = field[4]
                valueOut = field[5]
                valueIn_int = lerp(valueIn, valueOut, inCoeff)
                valueOut_int = lerp(valueIn, valueOut, outCoeff)
                avg = (valueIn_int + valueOut_int)/2

                cnd = cond(avg,Lreach,width,th_avg)

                fieldIdx = grid.dataProvider().fieldNameIndex('cond_%i' % sp)
                changedAttributes[fieldIdx] = cnd


            changedAttributesDict[feature.id()] = changedAttributes
##            except:
##                pass
                #print "Errore FID: %s" % feature.id()

    grid.dataProvider().changeAttributeValues(changedAttributesDict)
    grid.commitChanges()

    # Select only intersected cells (whose IDs are in ft_sel_id)
    grid.select(ft_sel_id)


    # Upload the vlayer into DB SQlite
    newName = newName + '_riv'
    uploadQgisVectorLayer(dbName, grid, newName, srid= grid.crs().postgisSrid(), selected = True )

    # Retrieve the Spatialite layer and add it to mapp
    uri = QgsDataSourceURI()
    uri.setDatabase(dbName)
    schema = ''
    table = newName
    geom_column = "Geometry"
    uri.setDataSource(schema, table, geom_column)
    display_name = table

    wlayer = QgsVectorLayer(uri.uri(), display_name, 'spatialite')

    QgsMapLayerRegistry.instance().addMapLayer(wlayer)

    deselectAll()


# ---------------------------

def createMultiRivLayer(newName, dbName, gridLayer, riverLayer, nsp):
    # --- Input:
    # newName: name of the new MDO for river
    # dbName: name of model DB (complete path)
    # gridLayer: grid layer
    # riverLayer: line layer
    # nsp: stress period number

    # retrieve information of the modeltable and the crs field
    layerNameList = getVectorLayerNames()
    (modelNameList, pathList) =  getModelsInfoLists(layerNameList)

    #match the exact model trough the list of all models
    db = [os.path.basename(dbName).replace('.sqlite', '')]
    for i in db:
        if i in modelNameList:
            modelName = 'modeltable_' + str(i)

    model = getVectorLayerByName(modelName)
    for ft in model.getFeatures():
        crs = ft['crs']

    #
    provider = gridLayer.dataProvider()
    grid = QgsVectorLayer("Polygon?crs=" + crs, "temporary_layer", "memory")
    pr = grid.dataProvider()


    # fields of original grid layer
    fld = gridLayer.dataProvider().fields()
    for f in fld:
    	pr.addAttributes([f])

    # get the selection of grid cells that intersect the river layer
    for i in gridLayer.getFeatures():
        for l in riverLayer.getFeatures():
            if i.geometry().intersects(l.geometry()):
                gridLayer.select(i.id())


    # create a request with the selection
    request = QgsFeatureRequest().setFilterFids(gridLayer.selectedFeaturesIds())


    # features from the original grid layer
    for feature in gridLayer.getFeatures(request):
        pr.addFeatures([feature])


    grid.startEditing()
    grid.commitChanges()



    # Add the new field layer (Integer )
    newfield = QgsField('layer', QVariant.Int)
    pr.addAttributes( [newfield] )


    # Add the new field segment' (Integer )
    newfield = QgsField('segment', QVariant.Int)
    pr.addAttributes( [newfield] )



    # Define List of new fields name
    fieldsList = ['length']
    for i in range(1, nsp+1):
        fieldsList.append('stage_'+ str(i) )
        fieldsList.append('rbot_'+ str(i) )
        fieldsList.append('cond_'+ str(i) )


    # Call the method to write new fields to be added
    createAttributes(grid, fieldsList)

    #  Start editing to write record
    grid.startEditing()

    changedAttributesDict = {}


    # takes multi features of the same river
    riverFeature = [f for f in riverLayer.getFeatures(QgsFeatureRequest(0))][0]
    riverFeature = [f for f in riverLayer.getFeatures()]
    riverLine = [f.geometry() for f in riverFeature]
    riverLength = [l.length() for l in riverLine]
    riverPoints = [p.asPolyline() for p in riverLine]
    #

    ft_sel_id = []
    ft_dict = {}

    for ids, segment in enumerate(riverLine):

        for feature in grid.getFeatures():
            gridPoly = feature.geometry()
            inCellRiverLength = 0.0
            inCoeff = 0.0
            outCoeff = 0.0


            if segment.intersects(gridPoly):
                ft_sel_id.append(feature.id())


                #
                inCellRiver = segment.intersection(gridPoly)
                # print "Grid ID %s - River %s - Lenght in cell: %s" % (feature.id(),riverFeature.id(),inCellRiver.length())
                if inCellRiver.wkbType() == QGis.WKBLineString:
                    linepoints = inCellRiver.asPolyline()
                    # first point
                    pointIn = linepoints[0]
                    # last point
                    pointOut = linepoints[len(linepoints)-1]
                    #print "Ingresso: %s %s" % (pointIn.x(),pointIn.y())
                    #print "Uscita: %s %s" % (pointOut.x(),pointOut.y())
                elif inCellRiver.wkbType() == QGis.WKBMultiLineString:
                    multilinepoints = inCellRiver.asMultiPolyline()
                    # first point of first list
                    pointIn = multilinepoints[0][0]
                    lastGroup = multilinepoints[len(multilinepoints)-1]
                    # last point of last list
                    pointOut = lastGroup[len(lastGroup)-1]
                    #print "Ingresso: %s %s" % (pointIn.x(),pointIn.y())
                    #print "Uscita: %s %s" % (pointOut.x(),pointOut.y())

                _, _, beforeVertexIndexpointIn, _, _ = segment.closestVertex(QgsPoint(pointIn))
                _, _, beforeVertexIndexpointOut, _, _ = segment.closestVertex(QgsPoint(pointOut))

                # for point in riverPoints:
                point = segment.asPolyline()
                beforeVertexpointIn = point[beforeVertexIndexpointIn]
                beforeVertexpointOut = point[beforeVertexIndexpointOut]
                if beforeVertexIndexpointIn == -1:
                    beforeVertexIndexpointIn = 0
                if beforeVertexIndexpointOut == -1:
                    beforeVertexIndexpointOut = 0

                stretchIn = point[0: beforeVertexIndexpointIn+1]
                stretchIn.append(pointIn)
                stretchInGeom = QgsGeometry.fromPolyline(stretchIn)
                currentInLength = stretchInGeom.length()
                stretchOut = point[0: beforeVertexIndexpointOut+1]
                stretchOut.append(pointOut)
                stretchOutGeom = QgsGeometry.fromPolyline(stretchOut)
                currentOutLength = stretchOutGeom.length()
                lung = segment.length()
                inCoeff = currentInLength / lung
                outCoeff = currentOutLength / lung
                changedAttributes = {}


                # Insert value to layer field (layer)
                fieldIdx = grid.dataProvider().fieldNameIndex('layer')
                changedAttributes[fieldIdx] = 1

                # Insert value to segment field (segment)
                fieldIdx = grid.dataProvider().fieldNameIndex('xyz')
                changedAttributes[fieldIdx] = ids + 1

                # Insert value to length field
                fieldIdx = grid.dataProvider().fieldNameIndex('length')
                changedAttributes[fieldIdx] = inCellRiver.length()

                for i in range(1,nsp+1):
                    # Insert value to segment field (segment)
                    fieldIdx = grid.dataProvider().fieldNameIndex('stage_' +  str(i))
                    changedAttributes[fieldIdx] = 0.1

                    fieldIdx = grid.dataProvider().fieldNameIndex('rbot_' +  str(i))
                    changedAttributes[fieldIdx] = 0.1

                    fieldIdx = grid.dataProvider().fieldNameIndex('cond_' +  str(i))
                    changedAttributes[fieldIdx] = 0.001


                if feature.id() not in ft_dict.keys():
                    ft_dict[feature.id()] = [[feature, changedAttributes]]
                else:
                    ft_dict[feature.id()].append([feature, changedAttributes])

                changedAttributesDict[feature.id()] = changedAttributes



    # Now, add (possible) multiple cells.
    for j, f in enumerate(ft_dict.keys()):
        if len(ft_dict[f]) > 1:
            # this is a multi-cell feature - take all the cells except the last one (already considered!)
            for i, ff in enumerate(ft_dict[f][:-1]):
                ## -- Create futures as copies of this feature
                # Values of dictionary are list of 2-sized list: [feature, changedAttributesDictionary]
                feat = QgsFeature(grid.fields())
                feat.setGeometry(ff[0].geometry())
                # Add attributes row and col
                attributesList = [0, 0, ff[0]['ROW'],ff[0]['COL']]
                feat.setAttributes(ff[0].attributes())

                # Add attributes stored in the dictionary
                for fieldIdx in ff[1].keys():
                    feat.setAttribute(fieldIdx, ff[1][fieldIdx])
                # Add feature
                (res, outFeats) = grid.dataProvider().addFeatures( [ feat ] )


    grid.dataProvider().changeAttributeValues(changedAttributesDict)
    grid.commitChanges()
#
    # Select only intersected cells (whose IDs are in ft_sel_id)
    # grid.select(ft_sel_id)
#
    # Upload the vlayer into DB SQlite
    newName = newName + '_riv'
    uploadQgisVectorLayer(dbName, grid, newName, srid= grid.crs().postgisSrid())

    # Retrieve the Spatialite layer and add it to mapp
    uri = QgsDataSourceURI()
    uri.setDatabase(dbName)
    schema = ''
    table = newName
    geom_column = "Geometry"
    uri.setDataSource(schema, table, geom_column)
    display_name = table

    wlayer = QgsVectorLayer(uri.uri(), display_name, 'spatialite')

    QgsMapLayerRegistry.instance().addMapLayer(wlayer)

    deselectAll()

# --- Create a Drain Layer
def createDrnLayer(newName, dbName, gridLayer, riverLayer, csvlayer, width, layer, xyz, nsp ):
    # --- Input:
    # newName: name of the new MDO for river
    # dbName: name of model DB (complete path)
    # gridLayer: grid layer
    # csvlayer: csv table with parameters (for each stress period)
    #           [sp, elev_in, elev_out, krb_in, krb_out, thk_in, thk_out]
    # riverLayer: line layer
    # width: river width, given as scalar (constant)
    # layer: layer where RIV is applied  (integer)
    # xyz: river segment id (integer)
    # nsp: stress period number

    # Check if river layer has the correct format
    # To do: change check of Wkb into Geometry type

    #retreive information of the modeltable and the crs field

    ##
    # trasformo csv layer in dictionary
    dp = csvlayer.dataProvider()
    csv_dict = {}

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

    # convert the values of the dictionary in usable float numbers if they contain ',' and are string of unicode
    for k in csv_dict:
        for i in range(len(csv_dict[k])):
            if type(csv_dict[k][i]) == str or type(csv_dict[k][i]) == unicode and ',' in csv_dict[k][i]:
                csv_dict[k][i] = float(csv_dict[k][i].replace(',', '.'))
                QApplication.processEvents()

    layerNameList = getVectorLayerNames()
    (modelNameList, pathList) =  getModelsInfoLists(layerNameList)

    #match the exact model trough the list of all models
    db = [os.path.basename(dbName).replace('.sqlite', '')]
    for i in db:
        if i in modelNameList:
            modelName = 'modeltable_' + str(i)

    model = getVectorLayerByName(modelName)
    for ft in model.getFeatures():
        crs = ft['crs']

    # Create a copy of Grid Layer as temporary layer
    #
    provider = gridLayer.dataProvider()
    grid = QgsVectorLayer("Polygon?crs=" + crs, "temporary_layer", "memory")
    pr = grid.dataProvider()

    # fields of original grid layer
    fld = gridLayer.dataProvider().fields()
    for f in fld:
    	pr.addAttributes([f])

    # get the selection of grid cells that intersect the river layer
    for i in gridLayer.getFeatures():
        for l in riverLayer.getFeatures():
            if i.geometry().intersects(l.geometry()):
                gridLayer.select(i.id())


    # create a request with the selection
    request = QgsFeatureRequest().setFilterFids(gridLayer.selectedFeaturesIds())


    # features from the original grid layer
    for feature in gridLayer.getFeatures(request):
        pr.addFeatures([feature])

    grid.startEditing()
    grid.commitChanges()

    # Add the new fields layer and segment (Integer)
    newfield = QgsField('layer', QVariant.Int)
    pr.addAttributes( [newfield] )
    newfield = QgsField('segment', QVariant.Int)
    pr.addAttributes( [newfield] )


    # Define List of new fields name
    fieldsList = ['length']
    for i in range(1,nsp+1):
        fieldsList.append('elev_'+ str(i) )
        fieldsList.append('cond_'+ str(i) )




    # Call the method to write new fields to be added
    createAttributes(grid, fieldsList)

    #  Start editing to write record
    grid.startEditing()

    changedAttributesDict = {}
    #
    ft_sel_id = []
    for feature in grid.getFeatures():
        gridPoly = feature.geometry()
        inCellRiverLength = 0.0
        inCoeff = 0.0
        outCoeff = 0.0
        # I take the only one river
        riverFeature = [f for f in riverLayer.getFeatures(QgsFeatureRequest(0))][0]
        riverLine = riverFeature.geometry()
        riverLength = riverLine.length()
        riverPoints = riverLine.asPolyline()
        #
        if riverLine.intersects(gridPoly):
            #
            ft_sel_id.append(feature.id())
            #
            inCellRiver = riverLine.intersection(gridPoly)
            # print "Grid ID %s - River %s - Lenght in cell: %s" % (feature.id(),riverFeature.id(),inCellRiver.length())
            if inCellRiver.wkbType() == QGis.WKBLineString:
                linepoints = inCellRiver.asPolyline()
                # first point
                pointIn = linepoints[0]
                # last point
                pointOut = linepoints[len(linepoints)-1]
                #print "Ingresso: %s %s" % (pointIn.x(),pointIn.y())
                #print "Uscita: %s %s" % (pointOut.x(),pointOut.y())
            elif inCellRiver.wkbType() == QGis.WKBMultiLineString:
                multilinepoints = inCellRiver.asMultiPolyline()
                # first point of first list
                pointIn = multilinepoints[0][0]
                lastGroup = multilinepoints[len(multilinepoints)-1]
                # last point of last list
                pointOut = lastGroup[len(lastGroup)-1]
                #print "Ingresso: %s %s" % (pointIn.x(),pointIn.y())
                #print "Uscita: %s %s" % (pointOut.x(),pointOut.y())

            _, _, beforeVertexIndexpointIn, _, _ = riverLine.closestVertex(QgsPoint(pointIn))
            _, _, beforeVertexIndexpointOut, _, _ = riverLine.closestVertex(QgsPoint(pointOut))
            beforeVertexpointIn = riverPoints[beforeVertexIndexpointIn]
            beforeVertexpointOut = riverPoints[beforeVertexIndexpointOut]
            #print "Indice vertice precedente Ingresso: %s" % (beforeVertexIndexpointIn)
            #print "Indice vertice precedente Uscita: %s" % (beforeVertexIndexpointOut)
            if beforeVertexIndexpointIn == -1:
                beforeVertexpointIn = pointIn
            if beforeVertexIndexpointOut == -1:
                beforeVertexpointOut = pointOut
            #print "Vertice precedente Ingresso: %s %s" % (beforeVertexpointIn.x(),beforeVertexpointIn.y())
            #print "Vertice precedente Uscita: %s %s" % (beforeVertexpointOut.x(),beforeVertexpointOut.y())

            if beforeVertexIndexpointIn == -1:
                beforeVertexIndexpointIn = 0
            if beforeVertexIndexpointOut == -1:
                beforeVertexIndexpointOut = 0

            stretchIn = riverPoints[0: beforeVertexIndexpointIn+1]
            stretchIn.append(pointIn)
            stretchInGeom = QgsGeometry.fromPolyline(stretchIn)
            currentInLength = stretchInGeom.length()
            stretchOut = riverPoints[0: beforeVertexIndexpointOut+1]
            stretchOut.append(pointOut)

##            try:
            stretchOutGeom = QgsGeometry.fromPolyline(stretchOut)
            currentOutLength = stretchOutGeom.length()
            inCoeff = currentInLength / riverLength
            outCoeff = currentOutLength / riverLength
##                print "Parziale ingresso: %s" % (currentInLength)
##                print "Parziale uscita: %s" % (currentOutLength)
##                print "Coefficiente interpolazioni ingresso: %s" % (float(currentInLength)/riverLine.length())
##                print "Coefficiente interpolazioni ingresso: %s" % (float(currentOutLength)/riverLine.length())
##                print "Grid ID: %s" % (feature.id())
            changedAttributes = {}
            # Insert value to layer field (layer)
            fieldIdx = grid.dataProvider().fieldNameIndex('layer')
            changedAttributes[fieldIdx] = layer

            # Insert value to segment field (segment)
            fieldIdx = grid.dataProvider().fieldNameIndex('segment')
            changedAttributes[fieldIdx] = xyz

            # Insert value to length field
            fieldIdx = grid.dataProvider().fieldNameIndex('length')
            changedAttributes[fieldIdx] = inCellRiver.length()
            Lreach = inCellRiver.length()

            for sp, field in csv_dict.iteritems():
                # elev_in, elev_out, krb_in, krb_out, thk_in, thk_out

                # elevation
                valueIn = field[0]
                valueOut = field[1]
                valueIn_int = lerp(valueIn, valueOut, inCoeff)
                valueOut_int = lerp(valueIn, valueOut, outCoeff)
                avg = (valueIn_int + valueOut_int)/2

                fieldIdx = grid.dataProvider().fieldNameIndex('elev_%i' % sp)
                changedAttributes[fieldIdx] = avg

                # thickness
                valueIn = field[4]
                valueOut = field[5]
                valueIn_int = lerp(valueIn, valueOut, inCoeff)
                valueOut_int = lerp(valueIn, valueOut, outCoeff)
                th_avg = (valueIn_int + valueOut_int)/2

                # conductance
                valueIn = field[2]
                valueOut = field[3]
                valueIn_int = lerp(valueIn, valueOut, inCoeff)
                valueOut_int = lerp(valueIn, valueOut, outCoeff)
                avg = (valueIn_int + valueOut_int)/2

                cnd = cond(avg,Lreach,width,th_avg)

                fieldIdx = grid.dataProvider().fieldNameIndex('cond_%i' % sp)
                changedAttributes[fieldIdx] = cnd

                QApplication.processEvents()



            changedAttributesDict[feature.id()] = changedAttributes
##            except:
##                pass
                #print "Errore FID: %s" % feature.id()

    grid.dataProvider().changeAttributeValues(changedAttributesDict)
    grid.commitChanges()

    # Select only intersected cells (whose IDs are in ft_sel_id)
    grid.select(ft_sel_id)


    # Upload the vlayer into DB SQlite
    newName = newName + '_drn'
    uploadQgisVectorLayer(dbName, grid, newName, srid= grid.crs().postgisSrid(), selected = True )

    # Retrieve the Spatialite layer and add it to mapp
    uri = QgsDataSourceURI()
    uri.setDatabase(dbName)
    schema = ''
    table = newName
    geom_column = "Geometry"
    uri.setDataSource(schema, table, geom_column)
    display_name = table

    wlayer = QgsVectorLayer(uri.uri(), display_name, 'spatialite')

    QgsMapLayerRegistry.instance().addMapLayer(wlayer)

    deselectAll()
# ---------------------------


def createMultiDrnLayer(newName, dbName, gridLayer, riverLayer, nsp):
    # --- Input:
    # newName: name of the new MDO for river
    # dbName: name of model DB (complete path)
    # gridLayer: grid layer
    # riverLayer: line layer
    # nsp: stress period number

    # retrieve information of the modeltable and the crs field

    layerNameList = getVectorLayerNames()
    (modelNameList, pathList) =  getModelsInfoLists(layerNameList)

    #match the exact model trough the list of all models
    db = [os.path.basename(dbName).replace('.sqlite', '')]
    for i in db:
        if i in modelNameList:
            modelName = 'modeltable_' + str(i)

    model = getVectorLayerByName(modelName)
    for ft in model.getFeatures():
        crs = ft['crs']

    #
    provider = gridLayer.dataProvider()
    grid = QgsVectorLayer("Polygon?crs=" + crs, "temporary_layer", "memory")
    pr = grid.dataProvider()


    # fields of original grid layer
    fld = gridLayer.dataProvider().fields()
    for f in fld:
    	pr.addAttributes([f])

    # get the selection of grid cells that intersect the river layer
    for i in gridLayer.getFeatures():
        for l in riverLayer.getFeatures():
            if i.geometry().intersects(l.geometry()):
                gridLayer.select(i.id())


    # create a request with the selection
    request = QgsFeatureRequest().setFilterFids(gridLayer.selectedFeaturesIds())


    # features from the original grid layer
    for feature in gridLayer.getFeatures(request):
        pr.addFeatures([feature])


    grid.startEditing()
    grid.commitChanges()


    # Add the new field layer (Integer )
    newfield = QgsField('layer', QVariant.Int)
    pr.addAttributes( [newfield] )

    # Add the new field segment' (Integer )
    newfield = QgsField('segment', QVariant.Int)
    pr.addAttributes( [newfield] )


    # Define List of new fields name
    fieldsList = ['length']
    for i in range(1, nsp + 1):
        fieldsList.append('elev_'+ str(i) )
        fieldsList.append('cond_'+ str(i) )


    # Call the method to write new fields to be added
    createAttributes(grid, fieldsList)

    #  Start editing to write record
    grid.startEditing()

    changedAttributesDict = {}


    # takes multi features of the same river
    riverFeature = [f for f in riverLayer.getFeatures(QgsFeatureRequest(0))][0]
    riverFeature = [f for f in riverLayer.getFeatures()]
    riverLine = [f.geometry() for f in riverFeature]
    riverLength = [l.length() for l in riverLine]
    riverPoints = [p.asPolyline() for p in riverLine]
    #

    ft_sel_id = []
    ft_dict = {}

    for ids, segment in enumerate(riverLine):

        for feature in grid.getFeatures():
            gridPoly = feature.geometry()
            inCellRiverLength = 0.0
            inCoeff = 0.0
            outCoeff = 0.0


            if segment.intersects(gridPoly):
                ft_sel_id.append(feature.id())


                #
                inCellRiver = segment.intersection(gridPoly)
                # print "Grid ID %s - River %s - Lenght in cell: %s" % (feature.id(),riverFeature.id(),inCellRiver.length())
                if inCellRiver.wkbType() == QGis.WKBLineString:
                    linepoints = inCellRiver.asPolyline()
                    # first point
                    pointIn = linepoints[0]
                    # last point
                    pointOut = linepoints[len(linepoints)-1]
                    #print "Ingresso: %s %s" % (pointIn.x(),pointIn.y())
                    #print "Uscita: %s %s" % (pointOut.x(),pointOut.y())
                elif inCellRiver.wkbType() == QGis.WKBMultiLineString:
                    multilinepoints = inCellRiver.asMultiPolyline()
                    # first point of first list
                    pointIn = multilinepoints[0][0]
                    lastGroup = multilinepoints[len(multilinepoints)-1]
                    # last point of last list
                    pointOut = lastGroup[len(lastGroup)-1]
                    #print "Ingresso: %s %s" % (pointIn.x(),pointIn.y())
                    #print "Uscita: %s %s" % (pointOut.x(),pointOut.y())

                _, _, beforeVertexIndexpointIn, _, _ = segment.closestVertex(QgsPoint(pointIn))
                _, _, beforeVertexIndexpointOut, _, _ = segment.closestVertex(QgsPoint(pointOut))

                # for point in riverPoints:
                point = segment.asPolyline()
                beforeVertexpointIn = point[beforeVertexIndexpointIn]
                beforeVertexpointOut = point[beforeVertexIndexpointOut]
                if beforeVertexIndexpointIn == -1:
                    beforeVertexIndexpointIn = 0
                if beforeVertexIndexpointOut == -1:
                    beforeVertexIndexpointOut = 0

                stretchIn = point[0: beforeVertexIndexpointIn+1]
                stretchIn.append(pointIn)
                stretchInGeom = QgsGeometry.fromPolyline(stretchIn)
                currentInLength = stretchInGeom.length()
                stretchOut = point[0: beforeVertexIndexpointOut+1]
                stretchOut.append(pointOut)
                stretchOutGeom = QgsGeometry.fromPolyline(stretchOut)
                currentOutLength = stretchOutGeom.length()
                lung = segment.length()
                inCoeff = currentInLength / lung
                outCoeff = currentOutLength / lung
                changedAttributes = {}


                # Insert value to layer field (layer)
                fieldIdx = grid.dataProvider().fieldNameIndex('layer')
                changedAttributes[fieldIdx] = 1

                # Insert value to segment field (segment)
                fieldIdx = grid.dataProvider().fieldNameIndex('xyz')
                changedAttributes[fieldIdx] = ids + 1

                # Insert value to length field
                fieldIdx = grid.dataProvider().fieldNameIndex('length')
                changedAttributes[fieldIdx] = inCellRiver.length()

                for i in range(1,nsp+1):
                    # Insert value to segment field (segment)
                    fieldIdx = grid.dataProvider().fieldNameIndex('elev_' +  str(i))
                    changedAttributes[fieldIdx] = 0.1

                    fieldIdx = grid.dataProvider().fieldNameIndex('cond_' +  str(i))
                    changedAttributes[fieldIdx] = 0.001


                if feature.id() not in ft_dict.keys():
                    ft_dict[feature.id()] = [[feature, changedAttributes]]
                else:
                    ft_dict[feature.id()].append([feature, changedAttributes])

                changedAttributesDict[feature.id()] = changedAttributes



    # Now, add (possible) multiple cells.
    for j, f in enumerate(ft_dict.keys()):
        if len(ft_dict[f]) > 1:
            # this is a multi-cell feature - take all the cells except the last one (already considered!)
            for i, ff in enumerate(ft_dict[f][:-1]):
                ## -- Create futures as copies of this feature
                # Values of dictionary are list of 2-sized list: [feature, changedAttributesDictionary]
                feat = QgsFeature(grid.fields())
                feat.setGeometry(ff[0].geometry())
                # Add attributes row and col
                attributesList = [0, 0, ff[0]['ROW'],ff[0]['COL']]
                feat.setAttributes(ff[0].attributes())

                # Add attributes stored in the dictionary
                for fieldIdx in ff[1].keys():
                    feat.setAttribute(fieldIdx, ff[1][fieldIdx])
                # Add feature
                (res, outFeats) = grid.dataProvider().addFeatures( [ feat ] )


    grid.dataProvider().changeAttributeValues(changedAttributesDict)
    grid.commitChanges()

    # Upload the vlayer into DB SQlite
    newName = newName + '_drn'
    uploadQgisVectorLayer(dbName, grid, newName, srid= grid.crs().postgisSrid())

    # Retrieve the Spatialite layer and add it to mapp
    uri = QgsDataSourceURI()
    uri.setDatabase(dbName)
    schema = ''
    table = newName
    geom_column = "Geometry"
    uri.setDataSource(schema, table, geom_column)
    display_name = table

    wlayer = QgsVectorLayer(uri.uri(), display_name, 'spatialite')

    QgsMapLayerRegistry.instance().addMapLayer(wlayer)

    deselectAll()




# --- Create a GHB Layer
def createGhbLayer(newName, dbName, gridLayer, lineLayer, csvlayer, bndDistance, from_lay, to_lay, xyz, nsp):
    # --- Input:
    # newName: name of the new MDO for river
    # dbName: name of model DB (complete path)
    # gridLayer: grid layer
    # csvlayer: csv table with parameters (for each stress period)
    #           [sp, rh_in, rh_out, bt_in, bt_out, krb_in, krb_out, thk_in, thk_out]
    # lineLayer: line layer
    # bndDistance: boundary distance, given as scalar (constant)
    # from_lay: layer (from) where GHB is applied  (integer)
    # to_lay:  layer (to) where GHB is applied  (integer)
    # xyz: river segment id (integer)
    # nsp: stress period number

    # Check if river layer has the correct format
    # To do: change check of Wkb into Geometry type

    ##
    # trasformo csv layer in dictionary
    dp = csvlayer.dataProvider()
    csv_dict = {}

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

    #retreive information of the modeltable and the crs field

    layerNameList = getVectorLayerNames()
    (modelNameList, pathList) =  getModelsInfoLists(layerNameList)

    #match the exact model trough the list of all models
    db = [os.path.basename(dbName).replace('.sqlite', '')]
    for i in db:
        if i in modelNameList:
            modelName = 'modeltable_' + str(i)

    model = getVectorLayerByName(modelName)
    for ft in model.getFeatures():
        crs = ft['crs']

    # Create a copy of Grid Layer as temporary layer
    #
    provider = gridLayer.dataProvider()
    grid = QgsVectorLayer("Polygon?crs=" + crs, "temporary_layer", "memory")
    pr = grid.dataProvider()

    # fields of original grid layer
    fld = gridLayer.dataProvider().fields()
    for f in fld:
    	pr.addAttributes([f])
    # features from the original grid layer
    for feature in processing.features(gridLayer):
        att = feature.attributes()
        pr.addFeatures([feature])
        QApplication.processEvents()

        # Retrieve feature's area and edge
        # ATTENTION !!!!!!
        # ASSUMPTION: grid is made by squared cells !!!!
        area = feature.geometry().area()
        edge = math.sqrt(area)

    grid.startEditing()
    grid.commitChanges()

    # Add the new fields 'layer' and 'segment' (Integer)
    newfield = QgsField('from_lay', QVariant.Int)
    pr.addAttributes( [newfield] )
    newfield = QgsField('to_lay', QVariant.Int)
    pr.addAttributes( [newfield] )
    newfield = QgsField('segment', QVariant.Int)
    pr.addAttributes( [newfield] )


    # Define List of new fields name
    fieldsList = []
    for i in range(1,nsp+1):
        fieldsList.append('bhead_'+ str(i) )
        fieldsList.append('cond_'+ str(i) )


    # Call the method to write new fields to be added
    createAttributes(grid, fieldsList)

    #  Start editing to write record
    grid.startEditing()

    changedAttributesDict = {}
    #
    ft_sel_id = []
    for feature in grid.getFeatures():
        QApplication.processEvents()
        gridPoly = feature.geometry()
        inCellRiverLength = 0.0
        inCoeff = 0.0
        outCoeff = 0.0
        # I take the only one river
        riverFeature = [f for f in lineLayer.getFeatures(QgsFeatureRequest(0))][0]
        riverLine = riverFeature.geometry()
        riverLength = riverLine.length()
        riverPoints = riverLine.asPolyline()
        #
        if riverLine.intersects(gridPoly):
            #
            ft_sel_id.append(feature.id())
            #
            inCellRiver = riverLine.intersection(gridPoly)
            # print "Grid ID %s - River %s - Lenght in cell: %s" % (feature.id(),riverFeature.id(),inCellRiver.length())
            if inCellRiver.wkbType() == QGis.WKBLineString:
                linepoints = inCellRiver.asPolyline()
                # first point
                pointIn = linepoints[0]
                # last point
                pointOut = linepoints[len(linepoints)-1]
                #print "Ingresso: %s %s" % (pointIn.x(),pointIn.y())
                #print "Uscita: %s %s" % (pointOut.x(),pointOut.y())
            elif inCellRiver.wkbType() == QGis.WKBMultiLineString:
                multilinepoints = inCellRiver.asMultiPolyline()
                # first point of first list
                pointIn = multilinepoints[0][0]
                lastGroup = multilinepoints[len(multilinepoints)-1]
                # last point of last list
                pointOut = lastGroup[len(lastGroup)-1]
                #print "Ingresso: %s %s" % (pointIn.x(),pointIn.y())
                #print "Uscita: %s %s" % (pointOut.x(),pointOut.y())

            _, _, beforeVertexIndexpointIn, _, _ = riverLine.closestVertex(QgsPoint(pointIn))
            _, _, beforeVertexIndexpointOut, _, _ = riverLine.closestVertex(QgsPoint(pointOut))
            beforeVertexpointIn = riverPoints[beforeVertexIndexpointIn]
            beforeVertexpointOut = riverPoints[beforeVertexIndexpointOut]
            #print "Indice vertice precedente Ingresso: %s" % (beforeVertexIndexpointIn)
            #print "Indice vertice precedente Uscita: %s" % (beforeVertexIndexpointOut)
            if beforeVertexIndexpointIn == -1:
                beforeVertexpointIn = pointIn
            if beforeVertexIndexpointOut == -1:
                beforeVertexpointOut = pointOut
            #print "Vertice precedente Ingresso: %s %s" % (beforeVertexpointIn.x(),beforeVertexpointIn.y())
            #print "Vertice precedente Uscita: %s %s" % (beforeVertexpointOut.x(),beforeVertexpointOut.y())

            if beforeVertexIndexpointIn == -1:
                beforeVertexIndexpointIn = 0
            if beforeVertexIndexpointOut == -1:
                beforeVertexIndexpointOut = 0

            stretchIn = riverPoints[0: beforeVertexIndexpointIn+1]
            stretchIn.append(pointIn)
            stretchInGeom = QgsGeometry.fromPolyline(stretchIn)
            currentInLength = stretchInGeom.length()
            stretchOut = riverPoints[0: beforeVertexIndexpointOut+1]
            stretchOut.append(pointOut)

##            try:
            stretchOutGeom = QgsGeometry.fromPolyline(stretchOut)
            currentOutLength = stretchOutGeom.length()
            inCoeff = currentInLength / riverLength
            outCoeff = currentOutLength / riverLength
##                print "Parziale ingresso: %s" % (currentInLength)
##                print "Parziale uscita: %s" % (currentOutLength)
##                print "Coefficiente interpolazioni ingresso: %s" % (float(currentInLength)/riverLine.length())
##                print "Coefficiente interpolazioni ingresso: %s" % (float(currentOutLength)/riverLine.length())
##                print "Grid ID: %s" % (feature.id())
            changedAttributes = {}
            # Insert value to from_lay field
            fieldIdx = grid.dataProvider().fieldNameIndex('from_lay')
            changedAttributes[fieldIdx] = from_lay

            # Insert value to to_lay field
            fieldIdx = grid.dataProvider().fieldNameIndex('to_lay')
            changedAttributes[fieldIdx] = to_lay

            # Insert value to segment field (segment)
            fieldIdx = grid.dataProvider().fieldNameIndex('segment')
            changedAttributes[fieldIdx] = xyz

            # Insert value to length field
            fieldIdx = grid.dataProvider().fieldNameIndex('length')
            changedAttributes[fieldIdx] = inCellRiver.length()
            Lreach = inCellRiver.length()

            for sp, field in csv_dict.iteritems():
                QApplication.processEvents()
                # bhe_in, bhe_out, hc_in, hc_out, thick_in, thick_out

                # boundary head
                bhe_in = field[0]
                bhe_out = field[1]
                bhe_in_int = lerp(bhe_in, bhe_out, inCoeff)
                bhe_out_int = lerp(bhe_in, bhe_out, outCoeff)
                bhe_avg = (bhe_in_int + bhe_out_int) / 2

                fieldIdx = grid.dataProvider().fieldNameIndex('bhead_%i' % sp)
                changedAttributes[fieldIdx] = bhe_avg

                # thickness for the conductance
                thick_in = field[4]
                thick_out = field[5]
                thick_in_int = lerp(thick_in, thick_out, inCoeff)
                thick_out_int = lerp(thick_in, thick_out, outCoeff)
                thick_avg = (thick_in_int + thick_out_int) / 2

                # conductance
                cond_in = field[2]
                cond_out = field[3]
                cond_in_int = lerp(cond_in, cond_out, inCoeff)
                cond_out_int = lerp(cond_in, cond_out, outCoeff)
                cond_avg = (cond_in_int + cond_out_int) / 2

                cnd = cond(cond_avg, thick_avg, edge, bndDistance)

                fieldIdx = grid.dataProvider().fieldNameIndex('cond_%i' % sp)
                changedAttributes[fieldIdx] = cnd

            changedAttributesDict[feature.id()] = changedAttributes

    grid.dataProvider().changeAttributeValues(changedAttributesDict)
    grid.commitChanges()

    # Select only intersected cells (whose IDs are in ft_sel_id)
    grid.select(ft_sel_id)


    # Upload the vlayer into DB SQlite
    newName = newName + '_ghb'
    uploadQgisVectorLayer(dbName, grid, newName, srid= grid.crs().postgisSrid(), selected = True )

    # Retrieve the Spatialite layer and add it to mapp
    uri = QgsDataSourceURI()
    uri.setDatabase(dbName)
    schema = ''
    table = newName
    geom_column = "Geometry"
    uri.setDataSource(schema, table, geom_column)
    display_name = table

    wlayer = QgsVectorLayer(uri.uri(), display_name, 'spatialite')

    QgsMapLayerRegistry.instance().addMapLayer(wlayer)

    deselectAll()
# ---------------------------



def createMultiGhbLayer(newName, dbName, gridLayer, riverLayer, nsp):
    # --- Input:
    # newName: name of the new MDO for river
    # dbName: name of model DB (complete path)
    # gridLayer: grid layer
    # riverLayer: line layer
    # nsp: stress period number

    # retrieve information of the modeltable and the crs field

    layerNameList = getVectorLayerNames()
    (modelNameList, pathList) =  getModelsInfoLists(layerNameList)

    #match the exact model trough the list of all models
    db = [os.path.basename(dbName).replace('.sqlite', '')]
    for i in db:
        if i in modelNameList:
            modelName = 'modeltable_' + str(i)

    model = getVectorLayerByName(modelName)
    for ft in model.getFeatures():
        crs = ft['crs']

    #
    provider = gridLayer.dataProvider()
    grid = QgsVectorLayer("Polygon?crs=" + crs, "temporary_layer", "memory")
    pr = grid.dataProvider()


    # fields of original grid layer
    fld = gridLayer.dataProvider().fields()
    for f in fld:
    	pr.addAttributes([f])

    # get the selection of grid cells that intersect the river layer
    for i in gridLayer.getFeatures():
        for l in riverLayer.getFeatures():
            if i.geometry().intersects(l.geometry()):
                gridLayer.select(i.id())


    # create a request with the selection
    request = QgsFeatureRequest().setFilterFids(gridLayer.selectedFeaturesIds())


    # features from the original grid layer
    for feature in gridLayer.getFeatures(request):
        pr.addFeatures([feature])


    grid.startEditing()
    grid.commitChanges()

    # Add the new fields 'layer' and 'xyz' (Integer)
    newfield = QgsField('from_lay', QVariant.Int)
    pr.addAttributes( [newfield] )
    newfield = QgsField('to_lay', QVariant.Int)
    pr.addAttributes( [newfield] )
    newfield = QgsField('segment', QVariant.Int)
    pr.addAttributes( [newfield] )


    # Define List of new fields name
    fieldsList = []
    for i in range(1,nsp+1):
        fieldsList.append('bhead_'+ str(i))
        fieldsList.append('cond_'+ str(i))

    # Call the method to write new fields to be added
    createAttributes(grid, fieldsList)


    #  Start editing to write record
    grid.startEditing()

    changedAttributesDict = {}


    # takes multi features of the same river
    riverFeature = [f for f in riverLayer.getFeatures(QgsFeatureRequest(0))][0]
    riverFeature = [f for f in riverLayer.getFeatures()]
    riverLine = [f.geometry() for f in riverFeature]
    riverLength = [l.length() for l in riverLine]
    riverPoints = [p.asPolyline() for p in riverLine]
    #

    ft_sel_id = []
    ft_dict = {}

    for ids, segment in enumerate(riverLine):

        for feature in grid.getFeatures():
            gridPoly = feature.geometry()
            inCellRiverLength = 0.0
            inCoeff = 0.0
            outCoeff = 0.0


            if segment.intersects(gridPoly):
                ft_sel_id.append(feature.id())


                #
                inCellRiver = segment.intersection(gridPoly)
                # print "Grid ID %s - River %s - Lenght in cell: %s" % (feature.id(),riverFeature.id(),inCellRiver.length())
                if inCellRiver.wkbType() == QGis.WKBLineString:
                    linepoints = inCellRiver.asPolyline()
                    # first point
                    pointIn = linepoints[0]
                    # last point
                    pointOut = linepoints[len(linepoints)-1]
                    #print "Ingresso: %s %s" % (pointIn.x(),pointIn.y())
                    #print "Uscita: %s %s" % (pointOut.x(),pointOut.y())
                elif inCellRiver.wkbType() == QGis.WKBMultiLineString:
                    multilinepoints = inCellRiver.asMultiPolyline()
                    # first point of first list
                    pointIn = multilinepoints[0][0]
                    lastGroup = multilinepoints[len(multilinepoints)-1]
                    # last point of last list
                    pointOut = lastGroup[len(lastGroup)-1]
                    #print "Ingresso: %s %s" % (pointIn.x(),pointIn.y())
                    #print "Uscita: %s %s" % (pointOut.x(),pointOut.y())

                _, _, beforeVertexIndexpointIn, _, _ = segment.closestVertex(QgsPoint(pointIn))
                _, _, beforeVertexIndexpointOut, _, _ = segment.closestVertex(QgsPoint(pointOut))

                # for point in riverPoints:
                point = segment.asPolyline()
                beforeVertexpointIn = point[beforeVertexIndexpointIn]
                beforeVertexpointOut = point[beforeVertexIndexpointOut]
                if beforeVertexIndexpointIn == -1:
                    beforeVertexIndexpointIn = 0
                if beforeVertexIndexpointOut == -1:
                    beforeVertexIndexpointOut = 0

                stretchIn = point[0: beforeVertexIndexpointIn+1]
                stretchIn.append(pointIn)
                stretchInGeom = QgsGeometry.fromPolyline(stretchIn)
                currentInLength = stretchInGeom.length()
                stretchOut = point[0: beforeVertexIndexpointOut+1]
                stretchOut.append(pointOut)
                stretchOutGeom = QgsGeometry.fromPolyline(stretchOut)
                currentOutLength = stretchOutGeom.length()
                lung = segment.length()
                inCoeff = currentInLength / lung
                outCoeff = currentOutLength / lung
                changedAttributes = {}


                # Insert value to segment field (xyz)
                fieldIdx = grid.dataProvider().fieldNameIndex('xyz')
                changedAttributes[fieldIdx] = ids + 1

                # Insert value to from_lay
                fieldIdx = grid.dataProvider().fieldNameIndex('from_lay')
                changedAttributes[fieldIdx] = 0

                # Insert value to to_lay
                fieldIdx = grid.dataProvider().fieldNameIndex('to_lay')
                changedAttributes[fieldIdx] = 0

                for i in range(1,nsp+1):
                    # Insert value to segment field (segment)
                    fieldIdx = grid.dataProvider().fieldNameIndex('bhead_' +  str(i))
                    changedAttributes[fieldIdx] = 0.1

                    fieldIdx = grid.dataProvider().fieldNameIndex('cond_' +  str(i))
                    changedAttributes[fieldIdx] = 0.001




                if feature.id() not in ft_dict.keys():
                    ft_dict[feature.id()] = [[feature, changedAttributes]]
                else:
                    ft_dict[feature.id()].append([feature, changedAttributes])

                changedAttributesDict[feature.id()] = changedAttributes



    # Now, add (possible) multiple cells.
    for j, f in enumerate(ft_dict.keys()):
        if len(ft_dict[f]) > 1:
            # this is a multi-cell feature - take all the cells except the last one (already considered!)
            for i, ff in enumerate(ft_dict[f][:-1]):
                ## -- Create futures as copies of this feature
                # Values of dictionary are list of 2-sized list: [feature, changedAttributesDictionary]
                feat = QgsFeature(grid.fields())
                feat.setGeometry(ff[0].geometry())
                # Add attributes row and col
                attributesList = [0, 0, ff[0]['ROW'],ff[0]['COL']]
                feat.setAttributes(ff[0].attributes())

                # Add attributes stored in the dictionary
                for fieldIdx in ff[1].keys():
                    feat.setAttribute(fieldIdx, ff[1][fieldIdx])
                # Add feature
                (res, outFeats) = grid.dataProvider().addFeatures( [ feat ] )


    grid.dataProvider().changeAttributeValues(changedAttributesDict)
    grid.commitChanges()

    # Upload the vlayer into DB SQlite
    newName = newName + '_ghb'
    uploadQgisVectorLayer(dbName, grid, newName, srid= grid.crs().postgisSrid())

    # Retrieve the Spatialite layer and add it to mapp
    uri = QgsDataSourceURI()
    uri.setDatabase(dbName)
    schema = ''
    table = newName
    geom_column = "Geometry"
    uri.setDataSource(schema, table, geom_column)
    display_name = table

    wlayer = QgsVectorLayer(uri.uri(), display_name, 'spatialite')

    QgsMapLayerRegistry.instance().addMapLayer(wlayer)

    deselectAll()


# --- Create a EVT Layer
def createEvtLayer(newName, dbName, gridLayer, csvlayer, nsp ):
    # --- Input:
    # newName: name of the new MDO for river
    # dbName: name of model DB (complete path)
    # gridLayer: grid layer
    # csvlayer: csv table with parameters (for each stress period)
    # nsp: stress period number

    ##
    # trasformo csv layer in dictionary
    dp = csvlayer.dataProvider()
    csv_dict = {}
    flag_autocsv = False
    if "ROWNO" in dp.fields()[0].name() :
        flag_autocsv = True

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

    name_fields = []
    type_fields = []
    default_fields = []

    for i in range(1,nsp+1):
        name_fields.append( 'surf_'+ str(i))
        type_fields.append(QVariant.Double)
        default_fields.append( csv_dict[i][0] )

        name_fields.append('evtr_'+ str(i))
        type_fields.append(QVariant.Double)
        default_fields.append(csv_dict[i][1])

        name_fields.append('exdp_'+ str(i))
        type_fields.append(QVariant.Double)
        default_fields.append(csv_dict[i][2])


    for i in range(1,nsp+1):
        name_fields.append( 'ievt_'+ str(i))
        type_fields.append(QVariant.Int)
        default_fields.append(1)


    createMDOLarge(gridLayer, dbName, newName, name_fields, type_fields, default_fields)



def createUzfLayer(newName, dbName, gridLayer, nsp ):
    # --- Input:
    # newName: name of the new MDO for uzf
    # dbName: name of model DB (complete path)
    # gridLayer: grid layer
    # nsp: stress period number

    # Define List of new fields name
    #

    QApplication.processEvents()

    name_fields = []
    type_fields = []
    default_fields = []

    name_fields.append('iuzfbnd')
    type_fields.append(QVariant.Int)
    default_fields.append(1)

    name_fields.append( 'eps')
    type_fields.append(QVariant.Double)
    default_fields.append(3.3)

    name_fields.append( 'thts')
    type_fields.append(QVariant.Double)
    default_fields.append(0.3)

    name_fields.append( 'thti')
    type_fields.append(QVariant.Double)
    default_fields.append(0.2)

    for i in range(1,nsp+1):
        name_fields.append( 'finf_'+ str(i))
        type_fields.append(QVariant.Double)
        default_fields.append(0)

        name_fields.append('pet_'+ str(i))
        type_fields.append(QVariant.Double)
        default_fields.append(0)

        name_fields.append('extdp_'+ str(i))
        type_fields.append(QVariant.Double)
        default_fields.append(0)

        name_fields.append('extwc_'+ str(i))
        type_fields.append(QVariant.Double)
        default_fields.append(0)



    createMDOLarge(gridLayer, dbName, newName, name_fields, type_fields, default_fields)


def createSmlLayer(newName, dbName, gridLayer, nsp):
    # --- Input:
    # newName: name of Surface Model Layer
    # dbName: name of model DB (complete path)
    # gridLayer: name of the grid

    # Define List of new fields name
    #
    name_fields = []
    type_fields = []
    default_fields = []

    # name_fields.append('manning')
    # type_fields.append(QVariant.Double)
    # default_fields.append(0.03)
    #
    # name_fields.append( 'slope')
    # type_fields.append(QVariant.Double)
    # default_fields.append(10.0)
    #
    # name_fields.append( 'aspect')
    # type_fields.append(QVariant.Double)
    # default_fields.append(100.0)

    name_fields.append('irunbnd')
    type_fields.append(QVariant.Int)
    default_fields.append(0)

    for i in range(1,nsp+1):
        name_fields.append( 'et_'+ str(i))
        type_fields.append(QVariant.Double)
        default_fields.append(0)



    createMDO(gridLayer, dbName, newName, name_fields, type_fields, default_fields)
# ---------------------------

def createSfrLayer(newName, dbName, gridLayer, riverLayer, layer, seg_id):
    # --- Input:
    # newName: name of the new MDO for river
    # dbName: name of model DB (complete path)
    # gridLayer: grid layer
    # river: line layer
    # width: river width, given as scalar (constant)
    # layer: layer where SFR is applied  (integer)
    # seg_id: river segment id (integer)

    # Check if river layer has the correct format
    # To do: change check of Wkb into Geometry type

    if riverLayer.wkbType() != 2:
        raise Exception("Input drain layer must by single part linestring")

    if riverLayer.featureCount() > 1:
        raise Exception("Input drain layer must contain only one feature")

    ##
    layerNameList = getVectorLayerNames()
    (modelNameList, pathList) =  getModelsInfoLists(layerNameList)

    #match the exact model trough the list of all models
    db = [os.path.basename(dbName).replace('.sqlite', '')]
    for i in db:
        if i in modelNameList:
            modelName = 'modeltable_' + str(i)

    model = getVectorLayerByName(modelName)
    for ft in model.getFeatures():
        crs = ft['crs']

    # Create a copy of Grid Layer as temporary layer
    #
    provider = gridLayer.dataProvider()
    grid = QgsVectorLayer("Polygon?crs=" + crs, "temporary_layer", "memory")
    pr = grid.dataProvider()

    # fields of original grid layer
    fld = gridLayer.dataProvider().fields()
    for f in fld:
    	pr.addAttributes([f])
    # features from the original grid layer
    for feature in processing.features(gridLayer):
        att = feature.attributes()
        pr.addFeatures([feature])

    grid.startEditing()
    grid.commitChanges()

    # Add the new fields 'layer', 'seg_id',  (Integer)
    newfield = QgsField('layer', QVariant.Int)
    pr.addAttributes( [newfield] )
    newfield = QgsField('seg_id', QVariant.Int)
    pr.addAttributes( [newfield] )
    newfield = QgsField('ireach', QVariant.Int)
    pr.addAttributes( [newfield] )
    newfield = QgsField('length', QVariant.Double)
    pr.addAttributes( [newfield] )






    # Call the method to write new fields to be added
    # createAttributes(grid, fieldsList)

    #  Start editing to write record
    grid.startEditing()

    changedAttributesDict = {}
    #
    ft_sel_id = []
    reach_iter = 0
    for feature in grid.getFeatures():
        QApplication.processEvents()
        gridPoly = feature.geometry()
        inCellRiverLength = 0.0
        inCoeff = 0.0
        outCoeff = 0.0
        # I take the only one river
        riverFeature = [f for f in riverLayer.getFeatures(QgsFeatureRequest(0))][0]
        riverLine = riverFeature.geometry()
        riverLength = riverLine.length()
        riverPoints = riverLine.asPolyline()
        #
        if riverLine.intersects(gridPoly):
            #
            ft_sel_id.append(feature.id())
            #
            reach_iter = reach_iter + 1
            #
            inCellRiver = riverLine.intersection(gridPoly)
            # print "Grid ID %s - River %s - Lenght in cell: %s" % (feature.id(),riverFeature.id(),inCellRiver.length())
            if inCellRiver.wkbType() == QGis.WKBLineString:
                linepoints = inCellRiver.asPolyline()
                # first point
                pointIn = linepoints[0]
                # last point
                pointOut = linepoints[len(linepoints)-1]
                #print "Ingresso: %s %s" % (pointIn.x(),pointIn.y())
                #print "Uscita: %s %s" % (pointOut.x(),pointOut.y())
            elif inCellRiver.wkbType() == QGis.WKBMultiLineString:
                multilinepoints = inCellRiver.asMultiPolyline()
                # first point of first list
                pointIn = multilinepoints[0][0]
                lastGroup = multilinepoints[len(multilinepoints)-1]
                # last point of last list
                pointOut = lastGroup[len(lastGroup)-1]
                #print "Ingresso: %s %s" % (pointIn.x(),pointIn.y())
                #print "Uscita: %s %s" % (pointOut.x(),pointOut.y())

            _, _, beforeVertexIndexpointIn, _, _ = riverLine.closestVertex(QgsPoint(pointIn))
            _, _, beforeVertexIndexpointOut, _, _ = riverLine.closestVertex(QgsPoint(pointOut))
            beforeVertexpointIn = riverPoints[beforeVertexIndexpointIn]
            beforeVertexpointOut = riverPoints[beforeVertexIndexpointOut]
            #print "Indice vertice precedente Ingresso: %s" % (beforeVertexIndexpointIn)
            #print "Indice vertice precedente Uscita: %s" % (beforeVertexIndexpointOut)
            if beforeVertexIndexpointIn == -1:
                beforeVertexpointIn = pointIn
            if beforeVertexIndexpointOut == -1:
                beforeVertexpointOut = pointOut
            #print "Vertice precedente Ingresso: %s %s" % (beforeVertexpointIn.x(),beforeVertexpointIn.y())
            #print "Vertice precedente Uscita: %s %s" % (beforeVertexpointOut.x(),beforeVertexpointOut.y())

            if beforeVertexIndexpointIn == -1:
                beforeVertexIndexpointIn = 0
            if beforeVertexIndexpointOut == -1:
                beforeVertexIndexpointOut = 0

            stretchIn = riverPoints[0: beforeVertexIndexpointIn+1]
            stretchIn.append(pointIn)
            stretchInGeom = QgsGeometry.fromPolyline(stretchIn)
            currentInLength = stretchInGeom.length()
            stretchOut = riverPoints[0: beforeVertexIndexpointOut+1]
            stretchOut.append(pointOut)

##            try:
            stretchOutGeom = QgsGeometry.fromPolyline(stretchOut)
            currentOutLength = stretchOutGeom.length()
            inCoeff = currentInLength / riverLength
            outCoeff = currentOutLength / riverLength
##                print "Parziale ingresso: %s" % (currentInLength)
##                print "Parziale uscita: %s" % (currentOutLength)
##                print "Coefficiente interpolazioni ingresso: %s" % (float(currentInLength)/riverLine.length())
##                print "Coefficiente interpolazioni ingresso: %s" % (float(currentOutLength)/riverLine.length())
##                print "Grid ID: %s" % (feature.id())
            changedAttributes = {}

            # Insert value to layer field (layer)
            fieldIdx = grid.dataProvider().fieldNameIndex('layer')
            changedAttributes[fieldIdx] = layer

            # Insert value to seg_id field (seg_id)
            fieldIdx = grid.dataProvider().fieldNameIndex('seg_id')
            changedAttributes[fieldIdx] = seg_id

            #Insert value to length field
            fieldIdx = grid.dataProvider().fieldNameIndex('length')
            #Lreach = inCellRiver.length()
            changedAttributes[fieldIdx] = inCellRiver.length()

            #Insert value to ireach field
            fieldIdx = grid.dataProvider().fieldNameIndex('ireach')
            #Lreach = inCellRiver.length()
            changedAttributes[fieldIdx] = reach_iter


            changedAttributesDict[feature.id()] = changedAttributes

    grid.dataProvider().changeAttributeValues(changedAttributesDict)
    grid.commitChanges()

    # Select only intersected cells (whose IDs are in ft_sel_id)
    grid.select(ft_sel_id)


    # Upload the vlayer into DB SQlite
    newName = newName + '_sfr'
    uploadQgisVectorLayer(dbName, grid, newName, srid= grid.crs().postgisSrid(), selected = True )

    # Retrieve the Spatialite layer and add it to mapp
    uri = QgsDataSourceURI()
    uri.setDatabase(dbName)
    schema = ''
    table = newName
    geom_column = "Geometry"
    uri.setDataSource(schema, table, geom_column)
    display_name = table

    wlayer = QgsVectorLayer(uri.uri(), display_name, 'spatialite')

    QgsMapLayerRegistry.instance().addMapLayer(wlayer)
# ---------------------------

def createSwiLayer(gridLayer, name, pathFile, modelName, nlay):

    QApplication.processEvents()

    layerName = name + "_swi"
    dbName = os.path.join(pathFile , modelName + '.sqlite')

    name_fields = []
    type_fields = []
    default_fields = []

    for i in range(1, nlay+1):
        name_fields.append('zeta_' + str(i))
        type_fields.append(QVariant.Double)
        default_fields.append(-10.0)
        name_fields.append('ssz_' + str(i))
        type_fields.append(QVariant.Double)
        default_fields.append(0.2)
        name_fields.append('isource_' + str(i))
        type_fields.append(QVariant.Int)
        default_fields.append(1)

    # Create MDO for SWI:
    createMDO(gridLayer, dbName, layerName, name_fields, type_fields, default_fields)
# ---------------------------
def createZoneLayer(gridLayer, name, pathFile, modelName, nlay):

    QApplication.processEvents()

    layerName = name + "_zone"
    dbName = pathFile + '/' + modelName + '.sqlite'

    name_fields = []
    type_fields = []
    default_fields = []

    for i in range(1, nlay+1):
        name_fields.append('zone_' + 'lay_' + str(i))
        type_fields.append(QVariant.Double)
        default_fields.append(1)

    # Create MDO for Zone:
    createMDO(gridLayer, dbName, layerName, name_fields, type_fields, default_fields)
# ---------------------------


# --- Create a Transport Layer
def createTransportLayer(gridLayer, name, pathFile, modelName, lay, nspec):

    layerName = name
    dbName = pathFile + '/' + modelName + '.sqlite'

    name_fields = []
    type_fields = []
    default_fields = []

    # Porosity
    #name_fields.append('PRSITY_' + str(i) )
    #type_fields.append(QVariant.Double)
    #default_fields.append(0.25)
    # Active (ICBUND )
    name_fields.append('ACTIVE' )
    type_fields.append(QVariant.Int)
    default_fields.append(1)
    # Longitudinal Dispersivity aL
    name_fields.append('LONG_D'  )
    type_fields.append(QVariant.Double)
    default_fields.append(10)
    # Ratio aT/aL
    name_fields.append('TRPT'  )
    type_fields.append(QVariant.Double)
    default_fields.append(1.0)
    # Ratio aV/aL
    name_fields.append('TRPV'  )
    type_fields.append(QVariant.Double)
    default_fields.append(0.3)

    for j in range(1,nspec +1):
        # Initial concentration (for each species)
        name_fields.append('SCONC_' + str(j) )
        type_fields.append(QVariant.Double)
        default_fields.append(0.0)
        # Molecular diffusion (for each species)
        name_fields.append('DMCOEF_' +  str(j) )
        type_fields.append(QVariant.Double)
        default_fields.append(0.00001)


    createMDO(gridLayer, dbName, layerName, name_fields, type_fields, default_fields)
# ---------------------------
# --- Create a Transport Layer
def createUztLayer(gridLayer, name, pathFile, modelName, nsp, nspec):

    layerName = name
    dbName = pathFile + '/' + modelName + '.sqlite'

    name_fields = []
    type_fields = []
    default_fields = []

    # Active (IUZFBND = 1 active, IUZFBND = 0 no active )
    name_fields.append('IUZFBND' )
    type_fields.append(QVariant.Int)
    default_fields.append(1)

    name_fields.append('WC' )
    type_fields.append(QVariant.Double)
    default_fields.append(0.2)

    for k in range(1, nsp+1):
        for j in range(1,nspec +1):
            # Concentration of the infiltrating flux (for each species)
            name_fields.append('CUZINF_sp_'+ str(k) + '_spec_' + str(j) )
            type_fields.append(QVariant.Double)
            default_fields.append(0.0)
            # concentration of ET fluxes originating from the unsaturated zone (for each species)
            name_fields.append('CUZET_sp_'+ str(k) + '_spec_' + str(j))
            type_fields.append(QVariant.Double)
            default_fields.append(0.0)

    createMDOLarge(gridLayer, dbName, layerName, name_fields, type_fields, default_fields)
# ---------------------------
#--- Create an SFT Layer

#The createSftLayer function requires the following inputs:
# - sfrLayer (the _sfr MDO (the _sft MDO will actually be a copy of this _sfr MDO))
# - name (the name to be assigned to the new _sft MDO)
# - pathFile (the path of the working directory)
# - modelName (the name of the model)
# - nspec (the number of simulated species)
#The output of the createSftLayer function is the _sft MDO
def createSftLayer(sfrLayer, name, pathFile, modelName, nspec):

    #The new MDO will be renamed as name_sft (name is an input to the createSftLayer function)
    layerName = name + "_sft"
    #dbName is something like pathfile/modelName.sqlite (pathfile and modelName are two inputs to the createSftLayer function)
    dbName = os.path.join(pathFile, modelName + '.sqlite')

    #retreive information of the modeltable and the crs field
    #The list layerNameList is filled with the names of all the vector layers available
    #in the Legend (getVectorLayerNames is a function of the freewat_utils.py module;
    #such function returns a list)
    layerNameList = getVectorLayerNames()
    #The tuple (modelNameList, pathList) is made of two lists:
    # -the list modelNameList is filled with the names of all the flow models available in the Legend
    # -the list pathList is filled with paths of working directories where the flow models DBs are saved
    #(getModelsInfoLists is a function of the freewat_utils.py module; such function
    #takes the list layerNameList as input, looks for layers whose names start with
    #modeltable and reads the content of the name and working_dir fields in the
    #modeltables found; it returns a tuple)
    (modelNameList, pathList) =  getModelsInfoLists(layerNameList)

    #match the exact model trough the list of all models
    #db is a list of the base names of all the models in the working directory
    db = [os.path.basename(dbName).replace('.sqlite', '')]
    #i iterates over the list db
    for i in db:
        #if i is also the name of one of the flow models available in the Legend...
        if i in modelNameList:
            #...then modelName is set to modeltable_i
            modelName = 'modeltable_' + str(i)

    #modelName is just a string like modeltable_i
    #Here, model is the layer object whose name is modelName
    #(getVectorLayerByName is a function of the freewat_utils.py module; such function returns a layer object)
    model = getVectorLayerByName(modelName)

    #SR iterates over the rows of the layer model (i.e., over the rows of the modeltable selected)
    for SR in model.getFeatures():
        #crs is the content of the cell in the modeltable whose heading is 'crs'
        crs = SR['crs']

    # Create a memory copy of the SfrLayer with only ROW and COL fields
    #Provider for the sfrLayer
    fld = sfrLayer.dataProvider().fields()

    #vl is a new temporary polygon layer to which the crs reference system is assigned
    vl = QgsVectorLayer("Polygon?crs=" + crs, "temporary_points", "memory")
    #Provider for the vl layer
    pr = vl.dataProvider()
    #The editing mode for the vl layer is activated
    vl.startEditing()

    #f iterates over the fields of the Attribute Table of the sfrLayer
    for f in fld:
        #All the fields of the sfrLayer are added in the Attribute Table of the vl layer
        pr.addAttributes([f])

    #feature iterates over the rows of the sfrLayer
    for feature in processing.features(sfrLayer):
        att = feature.attributes()
        #All the rows of the sfrLayer are added in the Attribute Table of the vl layer
        pr.addFeatures([feature])

    #Changes are saved in the Attribute Table of the vl layer
    vl.commitChanges()
    # ---
    #campi is the list of field names in the Attribute Table of the vl layer
    campi = pr.fields()
    #idcol is the index of the ROW field in the Attribute Table of the vl layer
    idcol = campi.indexFromName('ireach')
    #idx is an index going from idcol+1 until the total number of the fields in the Attribute Table of the vl layer
    #NOTE: idcol is the index for the layer field, as such, idx starts from the field following the ireach field
    idx = range(idcol+1, campi.count())
    #All fields in the Attribute Table of the vl layer with index equal to idx are deleted
    #(i.e., all fields except ROW, COL, layer, seg_id and ireach and the fields preceeding ROW (i.e., PKUID and ID))
    pr.deleteAttributes(idx)
    #Changes are saved in the Attribute Table of the vl layer
    vl.updateFields()
    #The vl layer is saved permanently
    newlayer = vl

    #The list name_fields will be filled with the names of further fields in the Attribute Table of the newlayer layer
    #ATTENTION: these fields will NOT be added in the Attribute Table of the newlayer layer! This will be done recalling
    #the createMDO function below. Here, we are just creating a list of fields to be added
    name_fields = []
    #The list type_fields will be filled with the type of further fields in the Attribute Table of the newlayer layer
    type_fields = []

    #The first field to be added in the Attribute Table of the newlayer layer is reach_ID (integer),
    #containing a progressive ID of stream reaches, according to the sorting assigned in the sfrLayer
    name_fields.append('reach_ID')
    type_fields.append(QVariant.Int)

    #default_fields is a list of default values which will be used to fill the fields in the
    #Attribute Table of the newlayer layer. Actually the fields PKUID, ID, ROW, COL, layer
    #are already filled with values copied from the Attribute Table of the sfrLayer.
    #The reach_ID field will have 1 as a default value
    default_fields = [1]

    #Other fields to be added in the Attribute Table of the newlayer layer are coldsf and dispsf (double)
    #for each species simulated in the transport model.
    #coldsf_m contains the starting concentration values for the m-th species (0 is the defaulsìt value);
    #dispsf_m contains the dispersion coefficient values for the m-th species (1e-5 is the default value)
    for j in range(1,nspec +1):
        name_fields.append('coldsf_spec_' + str(j))
        type_fields.append(QVariant.Double)
        default_fields.append(0.0)
        name_fields.append('dispsf_spec_' + str(j))
        type_fields.append(QVariant.Double)
        default_fields.append(0.00001)


    # Create MDO for SFT:
    #The createMDO function is recalled. The first argument of the createMDO function should be gridLayer,
    #as the createMDO function assumes that a selection of the grid cells where the new MDO is going to be
    #created is performed before recalling the createMDO function itself. In such case, the _sft MDO will
    #be created over the newlayer layer. As such, through the createMDO function, the Attribute Table of the
    #newlayer layer (which contains so far the PKUID, ID, ROW and COL fields) will be updated with the fields
    #whose names are contained in the list name_fields, whose types are contained in the list type_fields,
    #and whose default values are contained in the list default_fields
    createMDO(newlayer, dbName, layerName, name_fields, type_fields, default_fields)

# ---------------------------
#
# --- Create a Source Sink Layer
def createSSM_RchEvtLayer(gridLayer, name, pathFile, modelName, nsp, nlay, nspec, inrch, inevt):
    # inrch and inevt: list (with length nsp) of "T", true, or "F", false
    # they states if for that SP the package is activated (T) or not (F)

    layerName = name + "_zone_ssm"
    dbName = pathFile + '/' + modelName + '.sqlite'

    name_fields = []
    type_fields = []
    default_fields = []

    for k in range(1,nsp+1):
         if inrch[k-1] == 'T':
            for j in range(1,nspec +1):
                name_fields.append('CRCH_sp_'+ str(k) + '_spec_' + str(j) )
                type_fields.append(QVariant.Double)
                default_fields.append(0.0)
         if inevt[k-1] == 'T':
            for j in range(1,nspec +1):
                name_fields.append('CEVT_sp_'+ str(k) + '_spec_' + str(j) )
                type_fields.append(QVariant.Double)
                default_fields.append(0.0)


    createMDOLarge(gridLayer, dbName, layerName, name_fields, type_fields, default_fields)

# ---------------------------

# ---------------------------
#
# --- Create a Source Sink Layer
def createSSM_PointSourceLayer(gridLayer, name, whichtype, pathFile, modelName, splist, nsp, nlay, nspec):
    # whichtype is the string to set the layer name, depending on type of source/flow relationship
    #           It cna be: "chd", "wel", "riv", "ghb", "mass", "const"
    #

    layerName = name + "_%s_ssm"%whichtype

    dbName = pathFile + '/' + modelName + '.sqlite'

    name_fields = []
    type_fields = []
    default_fields = []

    for k in splist:
        name_fields.append('from_lay_'+ str(k) )
        type_fields.append(QVariant.Int)
        default_fields.append(1)
        name_fields.append('to_lay_'+ str(k) )
        type_fields.append(QVariant.Int)
        default_fields.append(1)

        for j in range(1,nspec +1):
            name_fields.append('CSSMS_sp_'+ str(k) + '_spec_' + str(j) )
            type_fields.append(QVariant.Double)
            default_fields.append(0.0)


    # For Modflow-associated sources, the grid layer is actually the MDO selected
    #      for the associated MFLOW package (chd,riv,wel,ghb).
    #      Thus, I have to take olny fields ROW and COL

    if (whichtype == 'mass' ) or (whichtype == 'const'):
        newlayer = gridLayer
    else:
        # Create a memory copy of the gridLayer with only ROW and COL fields
        fld = gridLayer.dataProvider().fields()

        vl = QgsVectorLayer("Polygon", "temporary_points", "memory")
        pr = vl.dataProvider()
        vl.startEditing()

        for f in fld:
        	pr.addAttributes([f])

        for feature in processing.features(gridLayer):
            att = feature.attributes()
            pr.addFeatures([feature])

        vl.commitChanges()
        # ---
        campi = pr.fields()
        idcol = campi.indexFromName('COL')
        idx = range(idcol+1, campi.count())
        pr.deleteAttributes(  idx   )
        vl.updateFields()
        newlayer = vl

    # Create the MDO
    createMDO(newlayer, dbName, layerName, name_fields, type_fields, default_fields)

# ---------------------------
# ---------------------------
#
# --- Create a Reaction Layer
def createReactionLayer(gridLayer, name, pathFile, modelName, nlay, nspec):

    layerName = name + "_rct"

    dbName = pathFile + '/' + modelName + '.sqlite'

    name_fields = []
    type_fields = []
    default_fields = []

    for i in range(1,nlay+1):
        #if isothm != 5 :
        name_fields.append('RHOB_lay_'+ str(i) )
        type_fields.append(QVariant.Double)
        default_fields.append(1700.0)
        #if isothm == 5 or isothm == 6:
        name_fields.append('PRSITY2_lay_' + str(i))
        type_fields.append(QVariant.Double)
        default_fields.append(0.2)

        for j in range(1,nspec +1):
            #if isothm > 0:
            name_fields.append('SP1_lay_' + str(i) + '_spec_' + str(j) )
            type_fields.append(QVariant.Double)
            default_fields.append(0.2)
            name_fields.append('SP2_lay_' + str(i) + '_spec_' +  str(j) )
            type_fields.append(QVariant.Double)
            default_fields.append(0.7)

            #if ireact > 0:
            name_fields.append('RC1_lay_' + str(i) + '_spec_' + str(j) )
            type_fields.append(QVariant.Double)
            default_fields.append(10.0)

            name_fields.append('RC2_lay_' + str(i) + '_spec_' +  str(j) )
            type_fields.append(QVariant.Double)
            default_fields.append(10.0)


    # Create the MDO
    createMDOLarge(gridLayer, dbName, layerName, name_fields, type_fields, default_fields)

# Farm Process Layers

def createFarmWellLayer(gridLayer, name, pathFile, modelName, nsp):

    layerName = name + "_farm_well"
    dbName = pathFile + '/' + modelName + '.sqlite'

    #retreive information of the modeltable and the crs field
    layerNameList = getVectorLayerNames()
    (modelNameList, pathList) =  getModelsInfoLists(layerNameList)

    #match the exact model trough the list of all models
    db = [os.path.basename(dbName).replace('.sqlite', '')]
    for i in db:
        if i in modelNameList:
            modelName = 'modeltable_' + str(i)

    model = getVectorLayerByName(modelName)

    for SR in model.getFeatures():
        crs = SR['crs']

    # Create a memory copy of the WellLayer with only ROW and COL fields
    fld = gridLayer.dataProvider().fields()

    vl = QgsVectorLayer("Polygon?crs=" + crs, "temporary_points", "memory")
    pr = vl.dataProvider()
    vl.startEditing()

    for f in fld:
        pr.addAttributes([f])

    for feature in processing.features(gridLayer):
        att = feature.attributes()
        pr.addFeatures([feature])

    vl.commitChanges()
    # ---
    campi = pr.fields()
    idcol = campi.indexFromName('ROW')
    idx = range(idcol+2, campi.count())
    pr.deleteAttributes(idx)
    vl.updateFields()
    newlayer = vl

    name_fields = []
    type_fields = []

    name_fields.append('Well_ID')
    type_fields.append(QVariant.Int)
    name_fields.append('Farm_ID')
    type_fields.append(QVariant.Int)
    name_fields.append('from_lay')
    type_fields.append(QVariant.Int)
    name_fields.append('to_lay')
    type_fields.append(QVariant.Int)


    default_fields = [1, 1, 1, 1]

    for i in range(1,nsp+1):
        name_fields.append('Qmax_'+ str(i))
        type_fields.append(QVariant.Double)
        default_fields.append(-100)


    # Create MDO for WEL:
    createMDO(newlayer, dbName, layerName, name_fields, type_fields, default_fields)
# ---------------------------

def createFarmId(gridLayer, name, pathFile, modelName, nlay):

    layerName = name + "_farm_id"
    dbName = pathFile + '/' + modelName + '.sqlite'

    name_fields = []
    type_fields = []
    default_fields = []

    name_fields.append('Farm_ID')
    type_fields.append(QVariant.Int)
    default_fields.append(1)

    # for i in range(1, nlay+1):
    #     name_fields.append('zone_' + 'lay_' + str(i))
    #     type_fields.append(QVariant.Double)
    #     default_fields.append(1)

    # Create MDO for Zone:
    createMDOLarge(gridLayer, dbName, layerName, name_fields, type_fields, default_fields)
# ---------------------------


def createFarmCropSoil(newName, dbName, gridLayer):
    # --- Input:
    # newName: name of Surface Model Layer
    # dbName: name of model DB (complete path)
    # gridLayer: name of the grid

    # Define List of new fields name
    #
    name_fields = []
    type_fields = []
    default_fields = []

    name_fields.append('soil_id')
    type_fields.append(QVariant.Int)
    default_fields.append(1)

    # name_fields.append('cap_fringe')
    # type_fields.append(QVariant.Double)
    # default_fields.append(1.0)

    name_fields.append('crop_id')
    type_fields.append(QVariant.Int)
    default_fields.append(1)


    createMDO(gridLayer, dbName, newName, name_fields, type_fields, default_fields)
# ---------------------------

def createFarmPipelineDiversion(gridLayer, name, pathFile, modelName, nlay):

    layerName = name + "_div"
    dbName = pathFile + '/' + modelName + '.sqlite'

    name_fields = []
    type_fields = []
    default_fields = []

    name_fields.append('farm_id')
    type_fields.append(QVariant.Int)
    default_fields.append(1)

    name_fields.append('segment')
    type_fields.append(QVariant.Int)
    default_fields.append(1)

    name_fields.append('ireach')
    type_fields.append(QVariant.Int)
    default_fields.append(1)

    # for i in range(1, nlay+1):
    #     name_fields.append('zone_' + 'lay_' + str(i))
    #     type_fields.append(QVariant.Double)
    #     default_fields.append(1)

    # Create MDO for Zone:
    createMDO(gridLayer, dbName, layerName, name_fields, type_fields, default_fields)
# ---------------------------

def createHobLayer(newName, dbName, gridLayer, pointLayer, pointFields, use_oat=False):
    # --- Input:
    # newName: name of the new MDO for river
    # dbName: name of model DB (complete path)
    # gridLayer: grid layer
    # pointLayer: point layer of Observations (wells)
    # pointFields: the list of names for fields required by point layer
    #              [Name, ZTop, ZBot ]

    # ft_lst = []
    # --
    # Create a copy of Grid Layer as temporary layer

    # retreive information of the modeltable and the crs field

    layerNameList = getVectorLayerNames()
    (modelNameList, pathList) =  getModelsInfoLists(layerNameList)

    #match the exact model trough the list of all models
    db = [os.path.basename(dbName).replace('.sqlite', '')]
    for i in db:
        if i in modelNameList:
            modelName = 'modeltable_' + str(i)
            timetableName = 'timetable_' + str(i)
    # Get model info from modeltable
    model = getVectorLayerByName(modelName)
    timetable =getVectorLayerByName(timetableName)

    for ft in model.getFeatures():
        crs = ft['crs']
        timeString = ft['time_unit']
        date0str = ft['initial_date']
        t0str = ft['initial_time']


    # Transform initial time and date in datetime objects, according to format
    date0 = datetime.datetime.strptime(date0str, "%Y-%m-%d")
    t0 = datetime.datetime.strptime(t0str,"%H:%M:%S")
    date0 = date0.date()
    t0 = t0.time()
    # Get info on Stress Periods and their lenght (used below to calculate time offset)
    nper = timetable.featureCount()
    perlen = [i for i in range(nper)]
    tid = 0
    for tf in timetable.getFeatures():
        perlen[tid] = float(tf['length'])
        tid += 1

    # Retrieve the List of layers
    myModel =  [os.path.basename(dbName).replace('.sqlite', '')][0]
    lpflayer = getVectorLayerByName('lpf_' + myModel)
    mLayerList = [getVectorLayerByName(f['name']) for f in lpflayer.getFeatures()]
    nlay = len(mLayerList)

    #
    provider = gridLayer.dataProvider()
    grid = QgsVectorLayer("Polygon?crs=" + crs, "temporary_layer", "memory")
    pr = grid.dataProvider()

    # fields of original grid layer
    fld = gridLayer.dataProvider().fields()
    for f in fld:
    	pr.addAttributes([f])
    # features from the original grid layer
    for feature in processing.features(gridLayer):
        att = feature.attributes()
        pr.addFeatures([feature])

    grid.startEditing()
    grid.commitChanges()

    # Add the new field 'NAME'
    newfield = QgsField('WELLNAME', QVariant.String)
    pr.addAttributes( [newfield] )

    newfield = QgsField('OBSNAM', QVariant.String)
    pr.addAttributes( [newfield] )
    newfield = QgsField('WEIGHT', QVariant.Double)
    pr.addAttributes( [newfield] )
    newfield = QgsField('STATISTICS', QVariant.String)
    pr.addAttributes( [newfield] )


    # Add the new field 'from_lay', 'to_lay'
    newfield = QgsField('from_lay', QVariant.Int)
    pr.addAttributes( [newfield] )
    newfield = QgsField('to_lay', QVariant.Int)
    pr.addAttributes( [newfield] )

    # Define List of new fields of type Double
    fieldsList = ['HOBS','ROFF','COFF']
    # Add PR_(n) where n is the number of layer
    for i in range(1,nlay+1):
        fieldsList.append('PR_'+ str(i) )
    #
    fieldsList.append('IREFSP')
    fieldsList.append('TOFFSET')

    # Call the method to write new fields of type Double
    createAttributes(grid, fieldsList)

    #  Start editing to write record
    grid.startEditing()

    changedAttributesDict = {}
    #
    ft_sel_id = []

    ### OAT read poi from sensor_freewat layer
    if use_oat:

        oat_stats = {}

        if pointLayer.selectedFeatureCount() < 1:
            QMessageBox.warning(iface.mainWindow(), 'Warning!', 'Please select at least one sensor')
            return

        grid_crs = grid.crs()
        point_crs = pointLayer.crs()
        transform = QgsCoordinateTransform(point_crs, grid_crs)

    ft_dict = {}
    for ft in grid.getFeatures():
    #for ft in processing.features(gridLayer):
        cellPoly = ft.geometry()
        cellCenter = cellPoly.centroid().asPoint()
        cellRect = cellPoly.boundingBox()
        delC = cellRect.width()
        delR = cellRect.height()
        ### changed poiLayer.getFeatures() with pointLayer.selectedFeatures()

        if use_oat:
            poi_features = pointLayer.selectedFeaturesIterator()
        else:
            poi_features = pointLayer.getFeatures()

        for fp in poi_features:
            pointGeom = fp.geometry()

            ### OAT if OAT source transform coordinates
            if use_oat:
                if fp['use'] == 0:
                    continue
                pointGeom.transform(transform)

            # Call the table of time-variant values, for this observation point
            #                named as the point itself
            # if fp['name'] not in layerNameList:
            #     QMessageBox.warning(self, self.tr('Warning!'), self.tr('Table for observation %s does not exists!' % fp['name']))
            #     return
            # else:
            #     obstable = getVectorLayerByName(fp['name'])

            if pointGeom.intersects(cellPoly):
                # Prepare dictionary to insert values within new fields
                changedAttributes = {}

                #
                myCellId = ft.id()
                ft_sel_id.append(myCellId)

                myPointFeature = pointGeom.intersection(cellPoly)
                myPoint = myPointFeature.asPoint()

                roff = (cellCenter.y() - myPoint.y())/delR
                coff = (myPoint.x() - cellCenter.x())/delC
##                print 'row', ft['ROW'], 'col', ft['COL'], roff, coff

                # Now, run through Model Layer List to calculate from_lay to_lay
                # and percentage of layer intersection (if any)
                # ztop and bottom are stored as 2nd and 3rd item (=[1,2]) in the list:
                ztop = fp[ pointFields[1]]
                zbot = fp[ pointFields[2]]
                #screen = ztop -zbot

                for m, mlayer in enumerate(mLayerList) :
                    lf = [ f for f in mlayer.getFeatures( QgsFeatureRequest().setFilterFid ( myCellId ) )][0]
                    l1, l2 = 0, 0
                    occupied = 0
                    top, bot = lf['TOP'], lf['BOTTOM']
                    thk = top - bot
                    thk_well = ztop - zbot
                    perTemp = 0

                    if ztop >= top and zbot <= bot:
                        perTemp = thk/thk_well
                    if ztop <  top and zbot > bot:
                        perTemp = 1.0
                    if ztop < top and ztop > bot and zbot <= bot:
                        perTemp = (ztop - bot)/thk_well
                    if ztop >= top and zbot < top and zbot > bot:
                        perTemp = (top - zbot)/thk_well

                    # Get from_lay and to_lay
                    if top >= ztop and bot <= ztop:
                        from_lay = m +1
                        occupied = 1
                        l1 = top - ztop
                    if top >= zbot and bot <= zbot:
                        to_lay = m +1
                        occupied = 1
                        l2 = zbot - bot
##                    # Update percentage
##                    perTemp  = occupied*(1.0-(l1+l2)/thk)

                    #perTemp = (l1+l2)/thk

                    # Insert value to PR field
                    fieldIdx = grid.dataProvider().fieldNameIndex('PR_' + str(m+1))
                    changedAttributes[fieldIdx] = perTemp

                # -- Now, compute for all observation in this point, the time offset
                # using initial date and time and the ones stored in observations table
                # It is OAT-compliant, named 'time'
                # for fv in obstable.getFeatures():
                    # timeStr = fv['time']
                    # # Convert in datetime objects
                    # timenowObject = datetime.datetime.strptime(timeStr, "%Y-%m-%d %H:%M:%S")
                    # datenow = timenowObject.date()
                    # timenow = timenowObject.time()
                    #
                    # # Convert in datetime objects
                    # # timenow = timenow.time()
                    # # Call the offset calculator
                    # [iref, offset, tsp] = timeoffsetCalc(t0, date0, timenow, datenow, nper, perlen, timeString)
                    # #
                    # ft_sel_id.append(ft.id())

                fieldIdx = grid.dataProvider().fieldNameIndex('from_lay')
                changedAttributes[fieldIdx] = from_lay
                #
                fieldIdx = grid.dataProvider().fieldNameIndex('to_lay')
                changedAttributes[fieldIdx] = to_lay
                #
                fieldIdx = grid.dataProvider().fieldNameIndex('ROFF')
                changedAttributes[fieldIdx] = roff
                #
                fieldIdx = grid.dataProvider().fieldNameIndex('COFF')
                changedAttributes[fieldIdx] = coff
                    #
                    # fieldIdx = grid.dataProvider().fieldNameIndex('TOFFSET')
                    # changedAttributes[fieldIdx] = offset
                    # #
                    # fieldIdx = grid.dataProvider().fieldNameIndex('IREFSP')
                    # changedAttributes[fieldIdx] = iref
                    # #
                    # fieldIdx = grid.dataProvider().fieldNameIndex('HOBS')
                    # # HOBS is Head, stored in table as 'data':
                    # changedAttributes[fieldIdx] = fv['data']
                    # #
                fieldIdx = grid.dataProvider().fieldNameIndex('WELLNAME')
                # NAME is stored as 1st item (=[0] ) in the list:
                changedAttributes[fieldIdx] = fp[pointFields[0]]

                ### OAT
                if use_oat:

                    oat_stats[fp[pointFields[0]]] = fp['statflag']

                #
                #changedAttributesDict[ft.id()] = changedAttributes
                changedAttributesDict[ft.id()] = changedAttributes

                # Dictionary to store (possible) multiple cells
                if ft.id() not in ft_dict.keys():
                    ft_dict[ft.id()] = [[ft, changedAttributes]]
                else:
                    ft_dict[ft.id()].append([ft, changedAttributes])
##                if ft.id() not in ft_dict.keys():
##                    ft_dict[ft.id()] = [[ft.geometry().asMultiPolygon()[0], ft.attributes(), changedAttributes]]
##                else:
##                    ft_dict[ft.id()].append([ft.geometry().asMultiPolygon()[0], ft.attributes(), changedAttributes])



    grid.dataProvider().changeAttributeValues(changedAttributesDict)
    grid.commitChanges()

    # Select only intersected cells (whose IDs are in ft_sel_id)
    grid.select(ft_sel_id)

    # Upload the vlayer into DB SQlite
    uploadQgisVectorLayer(dbName, grid, newName + 'temp_hob', srid=grid.crs().postgisSrid(), selected=True)

    # Retrieve the Spatialite layer and add it to mapp
    uri = QgsDataSourceURI()
    uri.setDatabase(dbName)
    schema = ''
    table = newName + 'temp_hob'
    geom_column = "Geometry"
    uri.setDataSource(schema, table, geom_column)
    display_name = table

    tmp_table = '{}temp_hob'.format(newName)

    wlayer = QgsVectorLayer(uri.uri(), display_name, 'spatialite')
    #
    del grid
    grid = wlayer
    grid.startEditing()
    # Now, add (possible) multiple cells.
    nfeats = grid.featureCount()
    for j, f in enumerate(ft_dict.keys()):
        if len(ft_dict[f]) > 1:
            # this is a multi-cell feature - take all the cells except the last one (already considered!)
            for i, ff in enumerate(ft_dict[f][:-1]):
                ## -- Create futures as copies of this feature
                # Values of dictionary are list of 2-sized list: [feature, changedAttributesDictionary]
                feat = QgsFeature(grid.fields())

                # Update the features counter
                nfeats += 1
                feat.id = nfeats

                # feat.setGeometry(new_geometry)
                feat.setGeometry(ff[0].geometry())

                # Here we need to change attributes 0, which is the primary key
                attributesList = ff[0].attributes()
                attributesList[0] = attributesList[0] + nfeats
                # Set attributes
                feat.setAttributes(attributesList)

                # Add attributes stored in the dictionary
                for fieldIdx in ff[1].keys():
                    feat.setAttribute(fieldIdx, ff[1][fieldIdx])

                # Add feature
                (res, outFeats) = grid.dataProvider().addFeatures( [ feat ] )
                del feat

    grid.commitChanges()
    # Create another memory copy, to add time-varying features
    del grid
    grid = QgsVectorLayer("Polygon?crs=" + crs, "temporary_layer", "memory")
    pr = grid.dataProvider()

    # fields of original grid layer
    fld = wlayer.dataProvider().fields()
    for f in fld:
        pr.addAttributes([f])

    for f in wlayer.getFeatures():

        # OAT: load data sensor table from db into a qgsLayer
        if use_oat:
            uri = QgsDataSourceURI()
            uri.setDatabase(dbName)
            table = f['WELLNAME']
            uri.setDataSource('', table, '')
            display_name = table

            obstable = QgsVectorLayer(uri.uri(), display_name, 'spatialite')

        else:

            if f['WELLNAME'] not in layerNameList:
                QMessageBox.warning(iface.mainWindow(), 'Table for observation %s does not exists!' % f['WELLNAME'])
                return
            else:
                # Get the table corresponding to this observation
                obstable = getVectorLayerByName(f['WELLNAME'])

        for fv in obstable.getFeatures():

            if use_oat:
                if fv['use'] == 0:
                    continue

            timeStr = fv['time']
            # Convert in datetime objects
            timenowObject = datetime.datetime.strptime(timeStr, "%Y-%m-%d %H:%M:%S")
            datenow = timenowObject.date()
            timenow = timenowObject.time()

            # Call the offset calculator
            if use_oat:
                [iref, offset, tsp] = timeoffsetCalcIsodate(t0, date0, timenow, datenow, nper, perlen, timeString)
            else:
                [iref, offset, tsp] = timeoffsetCalc(t0, date0, timenow, datenow, nper, perlen, timeString)
            #
            # Create a new feature for the new layer
            newfeature = QgsFeature()
            newfeature.setAttributes(f.attributes())
            newfeature.setGeometry(f.geometry())

            ### OAT: load correct data
            if use_oat:
                fieldIdx = grid.dataProvider().fieldNameIndex('OBSNAM')
                newfeature.setAttribute(fieldIdx, fv['obs_index'])
                #
                fieldIdx = wlayer.dataProvider().fieldNameIndex('HOBS')
                newfeature.setAttribute(fieldIdx, float(fv['data']))
                #
                fieldIdx = wlayer.dataProvider().fieldNameIndex('WEIGHT')
                newfeature.setAttribute(fieldIdx, float(fv['quality']))
                #
                fieldIdx = wlayer.dataProvider().fieldNameIndex('STATISTICS')
                newfeature.setAttribute(fieldIdx, oat_stats[f['WELLNAME']])
            else:
                fieldIdx = grid.dataProvider().fieldNameIndex('OBSNAM')
                newfeature.setAttribute(fieldIdx, fv['obsname'])
                #
                fieldIdx = wlayer.dataProvider().fieldNameIndex('HOBS')
                newfeature.setAttribute(fieldIdx, float(fv['data']))
                #
                fieldIdx = wlayer.dataProvider().fieldNameIndex('WEIGHT')
                newfeature.setAttribute(fieldIdx, float(fv['weight']))
                #
                fieldIdx = wlayer.dataProvider().fieldNameIndex('STATISTICS')
                newfeature.setAttribute(fieldIdx, fv['statistics'])
            #
            fieldIdx = wlayer.dataProvider().fieldNameIndex('TOFFSET')
            newfeature.setAttribute(fieldIdx, offset)
            #
            fieldIdx = grid.dataProvider().fieldNameIndex('IREFSP')
            newfeature.setAttribute(fieldIdx, iref)

            # Add feature to new layer
            pr.addFeatures([newfeature])

    # Upload the vlayer into DB SQlite
    newName = newName + '_hob'
    uploadQgisVectorLayer(dbName, grid, newName, srid=grid.crs().postgisSrid(), selected=False)
    # Retrieve the Spatialite layer and add it to mapp
    uri = QgsDataSourceURI()
    uri.setDatabase(dbName)
    schema = ''
    table = newName
    geom_column = "Geometry"
    uri.setDataSource(schema, table, geom_column)
    display_name = table

    final_layer = QgsVectorLayer(uri.uri(), display_name, 'spatialite')

    QgsMapLayerRegistry.instance().addMapLayer(final_layer)

    # remove tmp layer
    con = sqlite3.connect(dbName)
    cur = con.cursor()
    cur.execute("DROP TABLE IF EXISTS {}".format(tmp_table))

    cur.close()
    con.commit()
    con.close()


# ---------------------------
# --- Create a Flow Observations Layer
def createFlowObsLayer(newName, dbName, gridLayer, gageGroupNameField, segmentLayer, package, use_oat=False):
    # --- Input:
    # newName: name of the new MDO for Flow Observations
    # dbName: name of model DB (complete path)
    # gridLayer: layer of the MDO corresponding to Package (RIV, DRN or GHB)
    # gageGroupNameField: name of this Group of Flow Observation
    # segmentLayer: line layer of the gage segment
    # package: string. It can be RIV, DRN or GHB

    # Check if segment layer has the correct format
    # To do: change check of Wkb into Geometry type

    if segmentLayer.wkbType() != 2:
        raise Exception("Input gage segment layer must by single part linestring")

    # if segmentLayer.featureCount() > 1:
    #    raise Exception("Input gage segment layer must contain only one feature")

    # retreive information of the modeltable and the crs field

    layerNameList = getVectorLayerNames()
    (modelNameList, pathList) =  getModelsInfoLists(layerNameList)

    # match the exact model trough the list of all models
    db = [os.path.basename(dbName).replace('.sqlite', '')]
    for i in db:
        if i in modelNameList:
            modelName = 'modeltable_' + str(i)
            timetableName = 'timetable_' + str(i)

    model = getVectorLayerByName(modelName)
    timetable = getVectorLayerByName(timetableName)

    for ft in model.getFeatures():
        crs = ft['crs']
        timeString = ft['time_unit']
        date0str = ft['initial_date']
        t0str = ft['initial_time']

    # Transform initial time and date in datetime objects, according to format
    date0 = datetime.datetime.strptime(date0str, "%Y-%m-%d")
    t0 = datetime.datetime.strptime(t0str, "%H:%M:%S")
    date0 = date0.date()
    t0 = t0.time()

    # Get info on Stress Periods and their lenght (used below to calculate time offset)
    nper = timetable.featureCount()
    perlen = [i for i in range(nper)]
    tid = 0
    for tf in timetable.getFeatures():
        perlen[tid] = float(tf['length'])
        tid += 1

    # Create a memory copy of the MDO Layer with only ROW, COL and LAY info
    fld = gridLayer.dataProvider().fields()
    vl = QgsVectorLayer("Polygon?crs=" + crs, "temporary_layer", "memory")
    pr = vl.dataProvider()
    vl.startEditing()
    for f in fld:
        pr.addAttributes([f])
    for feature in processing.features(gridLayer):
        att = feature.attributes()
        pr.addFeatures([feature])
    vl.commitChanges()
    # ---
    campi = pr.fields()
    idcol = campi.indexFromName('COL')
    # If package = RIV or DRN, take ROW, COl, layer
    idx = range(idcol+2, campi.count())
    # If package = GHB, take ROW, COL, from_lay, to_lay
    if package == 'GHB':
        idx = range(idcol+2, campi.count())
    # delete other fields
    pr.deleteAttributes(idx)
    vl.updateFields()

    # Group ID for the gage segment
    newfield = QgsField('GAGENAME', QVariant.String)
    pr.addAttributes( [newfield] )
    newfield = QgsField('OBSNAM', QVariant.String)
    pr.addAttributes( [newfield] )
    newfield = QgsField('FLOWOBS', QVariant.Double)
    pr.addAttributes( [newfield] )
    newfield = QgsField('WEIGHT', QVariant.Double)
    pr.addAttributes( [newfield] )
    newfield = QgsField('STATISTICS', QVariant.String)
    pr.addAttributes( [newfield] )
    newfield = QgsField('TOFFSET', QVariant.Double)
    pr.addAttributes( [newfield] )
    newfield = QgsField('IREFSP', QVariant.Int)
    pr.addAttributes( [newfield] )

    # Start editing to write record
    vl.startEditing()

    changedAttributesDict = {}
    #
    ft_sel_id = []
    for feature in vl.getFeatures():
        # loop over MDO features (polygons)
        gridPoly = feature.geometry()
        # Loop over Gage segment layer (line)
        for fln in segmentLayer.getFeatures():
            riverFeature = fln
            riverLine = riverFeature.geometry()
            riverLength = riverLine.length()
            riverPoints = riverLine.asPolyline()
        #
            if riverLine.intersects(gridPoly):
                #
                ft_sel_id.append(feature.id())
                #
                inCellRiver = riverLine.intersection(gridPoly)

                changedAttributes = {}

                # -- Now, compute for this observation the time offset
                # using initial date and time and the ones stored in point layer
                # They are stored in field ids: [2, 3]
                # datenowStr = fln[segmentFields[2]]
                # timenowStr = fln[segmentFields[3]]
                # # Cambiare in lettura unica di date e time, in questo modo:
                # timeStr = fln[segmentFields[2]]
                # Convert in datetime objects
                # timenowObject = datetime.datetime.strptime(timeStr, "%d.%m.%Y %H:%M:%S")
                # datenow = timenowObject.date()
                # timenow = timenowObject.date()
                # Convert in datetime objects
                # datenow = datetime.datetime.strptime(datenowStr, "%d.%m.%Y")
                # timenow = datetime.datetime.strptime(timenowStr,"%H:%M:%S")
                # timenow = timenow.time()
                # # Call the offset calculator
                # [iref, offset, tsp] = timeoffsetCalc(t0, date0, timenow, datenow, nper, perlen, timeString)
                # #
                # fieldIdx = vl.dataProvider().fieldNameIndex('TOFFSET')
                # changedAttributes[fieldIdx] = offset
                # #
                # fieldIdx = vl.dataProvider().fieldNameIndex('IREFSP')
                # changedAttributes[fieldIdx] = iref
                # #
                # fieldIdx = vl.dataProvider().fieldNameIndex('ObsName')
                # # Observation NAME is stored as 2nd item (=[0] ) in the list:
                # changedAttributes[fieldIdx] = fln[ segmentFields[0] ]
                # #
                # fieldIdx = vl.dataProvider().fieldNameIndex('FlowObs')
                # # Flow Observation is stored as 3rd item (=[1] ) in the list:
                # changedAttributes[fieldIdx] = fln[ segmentFields[1] ]
                # #
                fieldIdx = vl.dataProvider().fieldNameIndex('GAGENAME')
                changedAttributes[fieldIdx] = fln[gageGroupNameField]

                changedAttributesDict[feature.id()] = changedAttributes

    vl.dataProvider().changeAttributeValues(changedAttributesDict)
    vl.commitChanges()

    # Select only intersected cells (whose IDs are in ft_sel_id)
    vl.select(ft_sel_id)

    # Upload the vlayer into DB SQlite
    uploadQgisVectorLayer(dbName, vl, newName + 'temp_flowob', srid=vl.crs().postgisSrid(), selected=True)

    # Retrieve the Spatialite layer and add it to mapp
    uri = QgsDataSourceURI()
    uri.setDatabase(dbName)
    schema = ''
    table = newName + 'temp_flowob'
    geom_column = "Geometry"
    uri.setDataSource(schema, table, geom_column)
    display_name = table

    wlayer = QgsVectorLayer(uri.uri(), display_name, 'spatialite')

    tmp_table = '{}temp_flowob'.format(newName)

    # Create another memory copy, to add time-varying features
    del vl
    grid = QgsVectorLayer("Polygon?crs=" + crs, "temporary_layer", "memory")
    pr = grid.dataProvider()

    # fields of original grid layer
    fld = wlayer.dataProvider().fields()
    for f in fld:
        pr.addAttributes([f])

    for f in wlayer.getFeatures():
        if f['GAGENAME'] not in layerNameList:
            if use_oat:
                uri = QgsDataSourceURI()
                uri.setDatabase(dbName)
                table = f['GAGENAME']
                uri.setDataSource('', table, '')
                display_name = table
                obstable = QgsVectorLayer(uri.uri(), display_name, 'spatialite')
            else:
                QMessageBox.warning(iface.mainWindow(), 'Warning!', 'Table for observation %s does not exists!' % f['GAGENAME'])
                return
        else:
            # Get the table corresponding to this observation
            obstable = getVectorLayerByName(f['GAGENAME'])

        for fv in obstable.getFeatures():
            timeStr = fv['time']
            # Convert in datetime objects
            timenowObject = datetime.datetime.strptime(timeStr, "%Y-%m-%d %H:%M:%S")
            datenow = timenowObject.date()
            timenow = timenowObject.time()

            # Call the offset calculator
            [iref, offset, tsp] = timeoffsetCalc(t0, date0, timenow, datenow, nper, perlen, timeString)
            #
            # Create a new feature for the new layer
            newfeature = QgsFeature()
            newfeature.setAttributes(f.attributes())
            newfeature.setGeometry(f.geometry())
            #
            fieldIdx = grid.dataProvider().fieldNameIndex('OBSNAM')
            newfeature.setAttribute(fieldIdx, fv['obsname'])
            #
            fieldIdx = wlayer.dataProvider().fieldNameIndex('FLOWOBS')
            newfeature.setAttribute(fieldIdx, float(fv['data']))
            #
            fieldIdx = wlayer.dataProvider().fieldNameIndex('WEIGHT')
            newfeature.setAttribute(fieldIdx, float(fv['weight']))
            #
            fieldIdx = wlayer.dataProvider().fieldNameIndex('STATISTICS')
            newfeature.setAttribute(fieldIdx, fv['statistics'])
            #
            fieldIdx = wlayer.dataProvider().fieldNameIndex('TOFFSET')
            newfeature.setAttribute(fieldIdx, offset)
            #
            fieldIdx = grid.dataProvider().fieldNameIndex('IREFSP')
            newfeature.setAttribute(fieldIdx, iref)

            # Add feature to new layer
            pr.addFeatures( [ newfeature ] )

    # Upload the vlayer into DB SQlite
    # build the name:
    if package == 'RIV':
        typeobs = '_rvob'
    if package == 'DRN':
        typeobs = '_drob'
    if package == 'GHB':
        typeobs = '_gbob'

    newName = newName + typeobs
    #
    uploadQgisVectorLayer(dbName, grid, newName, srid=grid.crs().postgisSrid(), selected=False)
    # Retrieve the Spatialite layer and add it to mapp
    uri = QgsDataSourceURI()
    uri.setDatabase(dbName)
    schema = ''
    table = newName
    geom_column = "Geometry"
    uri.setDataSource(schema, table, geom_column)
    display_name = table

    final_layer = QgsVectorLayer(uri.uri(), display_name, 'spatialite')

    QgsMapLayerRegistry.instance().addMapLayer(final_layer)

    # remove tmp layer

    con = sqlite3.connect(dbName)
    cur = con.cursor()
    cur.execute("DROP TABLE IF EXISTS {}".format(tmp_table))

    cur.close()
    con.commit()
    con.close()

# ---------------------------
# --- Create a USB Layer (as input for unsat solute balance calculator)
def createUsbLayer(gridLayer, name, pathFile, modelName, nsp):

    layerName = name + "_usb"
    dbName = pathFile + '/' + modelName + '.sqlite'

    # Initialize Default values for fields:
    name_fields = []
    type_fields = []
    default_fields = []

    name_fields.append('decay_cnst')
    type_fields.append(QVariant.Double)
    default_fields.append(0.001)

    for i in range(1,nsp+1):
        name_fields.append('init_conc_' + str(i))
        type_fields.append(QVariant.Double)
        default_fields.append(0.0)

    # Create MDO for USB:
    createMDO(gridLayer, dbName, layerName, name_fields, type_fields, default_fields)
