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

#******************************************************************************
#
# Freewat
# ---------------------------------------------------------
#
#
# Copyright (C) 2019 - 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.QtCore import QSettings
from PyQt4.QtGui import *
from qgis.gui import QgsGenericProjectionSelector
import os
from PyQt4 import QtGui, uic
#
import numpy as np
#
# Import from flopy
from flopy.modflow import Modflow
#
# module shapefile is saved under freewat module
from freewat.shapefile import Shape, Writer
from freewat.sqlite_utils import checkIfTableExists, uploadQgisVectorLayer
from freewat.freewat_utils import *
from freewat.mdoCreate_utils import createModel
import datetime
from pyspatialite import dbapi2 as sqlite3

# ---- This is taken from export module of Flopy, which some time seems not working properly.
def write_grid_shapefile(filename, sr, array_dict, nan_val=-1.0e9, selectedCell = None):
    """
    Write a grid shapefile array_dict attributes.

    Parameters
    ----------
    filename : string
        name of the shapefile to write
    sr : spatial reference instance
        spatial reference object for model grid
    array_dict : dict
       Dictionary of name and 2D array pairs.  Additional 2D arrays to add as
       attributes to the grid shapefile.

    selectedCell : if not None, is an array (nrow, ncol) with value = 1 for
       cells to be include, 0 otherwise.

    Returns
    -------
    None

    """

    try:
        import shapefile
    except Exception as e:
        raise Exception("io.to_shapefile(): error " +
                        "importing shapefile - try pip install pyshp")

    wr = shapefile.Writer(filename, shapeType=shapefile.POLYGON)
    wr.field("ROW", "N", 10, 0)
    wr.field("COL", "N", 10, 0)

    # manage selcetion option
    if selectedCell is None:
        selectedCell = np.ones(shape = (sr.nrow, sr.ncol))

    arrays = []
    names = list(array_dict.keys())
    names.sort()
    # for name,array in array_dict.items():
    for name in names:
        array = array_dict[name]
        if array.ndim == 3:
            assert array.shape[0] == 1
            array = array[0, :, :]
        assert array.shape == (sr.nrow, sr.ncol)
        array[np.where(np.isnan(array))] = nan_val
        if array.dtype in [np.int,np.int32,np.int64]:
            wr.field(name, "N", 10, 0)
        else:
            wr.field(name, "N", 20, 12)
        arrays.append(array)

    for i in range(sr.nrow):
        for j in range(sr.ncol):
            if selectedCell[i,j] > 0:
                pts = sr.get_vertices(i, j)
                wr.poly([pts])
                rec = [i + 1, j + 1]
                for array in arrays:
                    rec.append(array[i, j])
                wr.record(*rec)

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

#
FORM_CLASS, _ = uic.loadUiType(os.path.join( os.path.dirname(__file__), 'ui/ui_importModflow.ui') )

#class CreateModelDialog(QDialog, Ui_CreateModelDialog):
class ImportModflowModelDialog(QDialog, FORM_CLASS):
    def __init__(self, iface):
        QDialog.__init__(self)
        self.iface = iface
        self.setupUi(self)

        self.OutFilePath = None
        self.encoding = None
        self.buttonBox.rejected.connect(self.reject)
        self.buttonBox.button(QDialogButtonBox.Ok).clicked.connect(self.importModel)

        self.btnModflowDir.clicked.connect(self.modflowFile)
        self.btnWorkingDir.clicked.connect(self.workingDir)

    def modflowFile(self):
        (self.InputFilePath) = fileDialog(self)
        self.txtModflowFile.setText(str(self.InputFilePath))
        if self.InputFilePath is None or self.encoding is None:
            return

    def workingDir(self):
        (self.OutFilePath, self.encoding) = dirDialog(self)
        #see the current working directory in the text box
        self.txtDirectory.setText(self.OutFilePath)
        if self.OutFilePath is None or self.encoding is None:
            return

##    def restoreGui(self):
##        self.progressBar.setFormat("%p%")
##        self.progressBar.setRange(0, 1)
##        self.progressBar.setValue(0)
##        self.cancelButton.clicked.disconnect(self.stopProcessing)
##        self.okButton.setEnabled(True)

    def importModel(self):
        modelName = str(self.txtModelName.text())
        model_ws = self.txtDirectory.text()
        nam_file = os.path.basename(str(self.txtModflowFile.text()))
        dbname = os.path.join(model_ws, modelName + '.sqlite')
        #
        self.progressBar.setMinimum(0)
        # ------------------
        #
        if os.path.isfile(os.path.join(model_ws,modelName + '.sqlite')):
            pop_message(self.tr('In this Working Folder you halready have a model named:\n {} \n Please change directory or model name'.forma(modelName)), self.tr('warning'))
            return
        #
        m = Modflow()
        #
        ml = Modflow.load(nam_file, model_ws=model_ws, check=False)
        # Here we need to enter CRS and x_upperLeft, y_upperLeft
        # IF you have a rectangle representing the domain (e.g.the grid), you can load it in QGIS to get geo-information
        # as follows:
        ##vl = iface.activeLayer()
        ##xul, yul  = vl.extent().xMinimum(), vl.extent().yMaximum()
        ##proj4_str = vl.crs().toProj4 ()
        if self.anchorGroup.isChecked():
            xul = float(self.txtX.text())
            yul = float(self.txtY.text())
            rotation = float(self.txtAngle.text())
            # get Projection as string from the Qgs object
            proj4_str = str(self.projApp.crs().toProj4())
            crs = str(self.projApp.crs().authid())

        else:
            # set default values
            xul = 0.0
            yul = 0.0
            rotation = 0.0
            proj4_str = '+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs'
            crs = 'EPSG:32632'

        # DIS
        ml.dis.sr.proj4_str = proj4_str
        ml.dis.sr.xul = xul
        ml.dis.sr.yul = yul
        ml.dis.sr.rotation = rotation

        nrow, ncol, nlay = ml.nrow , ml.ncol, ml.nlay
        #
        self.progressBar.setValue(1)

        # --
        # Create Spatialite DB of your model
        lengthInteger = ml.dis.lenuni
        timeInteger = ml.dis.itmuni
        # convert to string, for FREEWAT table
        # length unit
        if lengthInteger == 0:
            lengthString  = 'undefined'
        elif lengthInteger == 1:
            lengthString  = 'ft'
        elif lengthInteger == 2:
            lengthString  = 'm'
        elif lengthInteger == 3:
            lengthString  = 'cm'
        # time unit
        if timeInteger == 0:
            timeString = 'undefined'
        elif timeInteger == 1:
            timeString = 'sec'
        elif timeInteger == 2:
            timeString = 'min'
        elif timeInteger == 3:
            timeString = 'hour'
        elif timeInteger == 4:
            timeString = 'day'
        elif timeInteger == 5:
            timeString = 'year'

        isChild = 1.0
        # info on stress period
        kper = ml.kper
        length = ml.dis.perlen.array
        time_steps = ml.dis.nstp.array
        multiplier = ml.dis.tsmult.array
        state = []
        for s in [ss for i, ss in enumerate(ml.dis.steady.array)]:
            if s :
                state.append('SS')
            else:
                state.append('TR')

        # first stress period info for defining the FREEWAT model
        timeParameters  = [float(length[0]), int(time_steps[0]), float(multiplier[0]), state[0]]

        # Info on initial time: this is fake!
        initTimeString = '12:00'
        initDateString = '2019-01-01'
        #
        self.progressBar.setValue(2)
        #
        createModel(modelName, model_ws, isChild, lengthString, timeString, timeParameters, initDateString, initTimeString,  crs)
        #
        self.progressBar.setValue(3)
        # --
        # start retriving info on the grid
        # grid
        fileshp = os.path.join(model_ws, 'model_shapefiles'  ,  'mygrid_grid.shp')
        wr = write_grid_shapefile(fileshp, ml.dis.sr, {'top': ml.dis.top.array})
        #
        # 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()
        #
        #
        if self.chkSqlite.isChecked():
            # temporarly disable the pop-up asking CRS
            settings = QSettings()
            oldProjValue = settings.value( "/Projections/defaultBehaviour", "prompt")
            settings.setValue( "/Projections/defaultBehaviour", "useProject" )
            layerName = 'model_grid'
            wb = QgsVectorLayer(fileshp, layerName, 'ogr')
            srs = QgsCoordinateReferenceSystem(unicode(crs))
            wb.setCrs(srs)
            srid = wb.crs().postgisSrid()
            uploadQgisVectorLayer(dbname, wb, layerName, srid = srid)
            #
            # Retrieve the Spatialite layer and add it to mapp
            uri = QgsDataSourceURI()
            uri.setDatabase(dbname)
            schema = ''
            table = layerName
            geom_column = "Geometry"
            uri.setDataSource(schema, table, geom_column)
            display_name = table
            wlayer = QgsVectorLayer(uri.uri(), display_name, 'spatialite')
            QgsMapLayerRegistry.instance().addMapLayer(wlayer)
            # revert settings to original
            settings.setValue( "/Projections/defaultBehaviour", oldProjValue )
        #
        # Update TIME TABLE with other stress periods (if any): all info is already uploaded!
        timeTable = 'timetable_' + modelName
        updateList = []
        for i, ln in enumerate(length[1:]):
            updateList.append((i+2,float(ln), int(time_steps[i+1]) , float(multiplier[i+1]), state[i+1] ))

        cur.executemany('INSERT INTO {} (sp, length, ts, multiplier, state) VALUES (?, ?, ?, ?, ?);'.format(timeTable), updateList)

        # Create new LPF table in DB
        nameTable = "lpf_"+ modelName
        ##        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)

        # prepare for LPF table insert
        layType = []
        for i, lt in enumerate(ml.lpf.laytyp.array):
            if lt == 1:
                layType.append('convertible')
            else:
                layType.append('confined')
        #
        layAverage = []
        for i, lv in enumerate(ml.lpf.layavg.array):
            if lv == 0:
                layAverage.append('harmonic')
            elif lv == 1:
                layAverage.append('logarithmic')
            else:
                layAverage.append('arithmetic-mean')
        # Remark: in FREEWAT chani is always 1.0
        chani = [1.0 for i in ml.lpf.chani.array]
        #
        layWet = []
        for i, lw in enumerate(ml.lpf.laywet.array):
            if lw == 1:
                layWet.append('Yes')
            else:
                layWet.append('No')
        #
        #
        self.progressBar.setValue(5)
        # model layers
        for l in range(nlay):
            fileshp = os.path.join(model_ws, 'model_shapefiles'  ,  'model_layer_%i.shp'%(l+1))
            modelLayerDict = {}
            # DIS
            'BORDER' ,'ACTIVE','TOP', 'BOTTOM', 'THICKNESS', 'STRT', 'KX', 'KY', 'KZ', 'SS', 'SY', 'NT', 'NE', ''
            if l == 0:
                modelLayerDict['TOP'] = ml.dis.top.array
            else:
                modelLayerDict['TOP'] = ml.dis.botm.array[l-1]
            modelLayerDict['BOTTOM'] = ml.dis.botm.array[l]
            # BAS
            modelLayerDict['ACTIVE'] = ml.bas6.ibound.array[l]
            modelLayerDict['STRT'] = ml.bas6.strt.array[l]
            modelLayerDict['KX'] = ml.lpf.hk.array[l] # KX = hk
            modelLayerDict['KY'] = ml.lpf.hk.array[l] # KY = KX
            modelLayerDict['KZ'] = ml.lpf.vka.array[l] # KZ = vka
            modelLayerDict['SS'] = ml.lpf.ss.array[l]
            modelLayerDict['SY'] = ml.lpf.sy.array[l]
            modelLayerDict['WETDRY'] = ml.lpf.wetdry.array[l]
            modelLayerDict['NT'] = 0.2*np.ones(shape = (nrow,ncol))
            modelLayerDict['NE'] = 0.2*np.ones(shape = (nrow,ncol))


            wr = write_grid_shapefile(fileshp, ml.dis.sr, modelLayerDict)
            # update LPF table
            # Insert values in LPF Table
            ftinsert = ('model_layer_%i'%(l+1),layType[l],layAverage[l],chani[l],layWet[l])
            sql3 = 'INSERT INTO %s '%nameTable +  '(name,type,layavg,chani,laywet)'
            cur.execute(sql3 + 'VALUES (?, ?, ?, ?, ?);', ftinsert)
            #
            self.progressBar.setValue(5 + l)

        # commit and close connection
        # Close SpatiaLiteDB
        cur.close()
        # Save the changes
        con.commit()
        # Close connection
        con.close()
        # Upload table into map
        # Add the model table into QGis map
        uri = QgsDataSourceURI()
        uri.setDatabase(dbname)
        schema = ''
        table = nameTable
        #geom_column = 'the_geom'
        geom_column = None
        uri.setDataSource(schema, table, geom_column)
        display_name = nameTable
        tableLayer = QgsVectorLayer(uri.uri(), display_name, 'spatialite')
        QgsMapLayerRegistry.instance().addMapLayer(tableLayer)
        #
        val = self.progressBar.value()
        self.progressBar.setValue(val + nlay + 1)

        # Upload the layers in DB:
        # this could be done inside the previous loop, but we need to chenge the method asking
        # only for cur and for dbName
        #
        if self.chkSqlite.isChecked():
            # temporarly disable the pop-up asking CRS
            settings = QSettings()
            oldProjValue = settings.value( "/Projections/defaultBehaviour", "prompt")
            settings.setValue( "/Projections/defaultBehaviour", "useProject" )
            for l in range(nlay):
                layerName = 'model_layer_%i'%(l+1)
                fileshp = os.path.join(model_ws, 'model_shapefiles'  ,  'model_layer_%i.shp'%(l+1))
                wb = QgsVectorLayer(fileshp, layerName, 'ogr')
                srs = QgsCoordinateReferenceSystem(unicode(crs))
                wb.setCrs(srs)
                srid = wb.crs().postgisSrid()
                uploadQgisVectorLayer(dbname, wb, layerName, srid = srid)
                #
                # Retrieve the Spatialite layer and add it to mapp
                uri = QgsDataSourceURI()
                uri.setDatabase(dbname)
                schema = ''
                table = layerName
                geom_column = "Geometry"
                uri.setDataSource(schema, table, geom_column)
                display_name = table
                wlayer = QgsVectorLayer(uri.uri(), display_name, 'spatialite')
                QgsMapLayerRegistry.instance().addMapLayer(wlayer)
                #
                val = self.progressBar.value()
                self.progressBar.setValue(val + 1)

            # revert settings to original
            settings.setValue( "/Projections/defaultBehaviour", oldProjValue )

        # --
        # Scout other packages
        pklist = [p.name[0] for p in ml.packagelist]
        nper  = ml.nper
        # RCH
        if 'RCH' in pklist:
            fileshp = os.path.join(model_ws, 'model_shapefiles'  ,  'model_rch.shp')
            rch_dict = {}
        ##    if ml.rch.irc is None:
        ##        irch = np.ones(shape = (nrow, ncol))
        ##    else:
            for t in range(nper):
                if ml.rch.irch is None:
                    irch = np.ones(shape = (nrow, ncol), dtype = np.int)
                else:
                    irch = ml.rch.irch[t].array
                rch_dict['sp_%i_irch'%(t+1)] = irch
                rch_dict['sp_%i_rech'%(t+1)] = ml.rch.rech[t].array

            wr = write_grid_shapefile(fileshp, ml.dis.sr, rch_dict)
            #
            if self.chkSqlite.isChecked():
                # temporarly disable the pop-up asking CRS
                settings = QSettings()
                oldProjValue = settings.value( "/Projections/defaultBehaviour", "prompt")
                settings.setValue( "/Projections/defaultBehaviour", "useProject" )
                layerName = 'model_rch'
                wb = QgsVectorLayer(fileshp, layerName, 'ogr')
                srs = QgsCoordinateReferenceSystem(unicode(crs))
                wb.setCrs(srs)
                srid = wb.crs().postgisSrid()
                uploadQgisVectorLayer(dbname, wb, layerName, srid = srid)
                #
                # Retrieve the Spatialite layer and add it to mapp
                uri = QgsDataSourceURI()
                uri.setDatabase(dbname)
                schema = ''
                table = layerName
                geom_column = "Geometry"
                uri.setDataSource(schema, table, geom_column)
                display_name = table
                wlayer = QgsVectorLayer(uri.uri(), display_name, 'spatialite')
                QgsMapLayerRegistry.instance().addMapLayer(wlayer)
                # revert settings to original
                settings.setValue( "/Projections/defaultBehaviour", oldProjValue )
            #
            val = self.progressBar.value()
            self.progressBar.setValue(val + 2)


        # WEL
        if 'WEL' in pklist:
            fileshp = os.path.join(model_ws, 'model_shapefiles'  ,  'model_well.shp')
            wel_dict = {}
            welCells = np.zeros(shape = (nrow, ncol))
            fromLayarray = np.ones(shape = (nrow,ncol), dtype = np.int)
            toLayarray = np.ones(shape = (nrow,ncol), dtype = np.int)
            for t in range(nper):
                nwells = len(ml.wel.stress_period_data[t])
                pumpingarray = np.zeros(shape = (nrow,ncol))
                for i in range(len(ml.wel.stress_period_data[t])):
                    (lay, r, c, f )= ml.wel.stress_period_data[t][i]
                    if fromLayarray[r,c] > lay + 1:
                        fromLayarray[r,c] = np.int(lay + 1)
                    if toLayarray[r,c] < lay + 1:
                        toLayarray[r,c] = np.int(lay + 1)
                    pumpingarray[r,c] = f
                    welCells[r,c] = 1
                wel_dict['sp_%i'%(t+1)] = pumpingarray
            wel_dict['from_lay'] = fromLayarray
            wel_dict['to_lay'] = toLayarray
            wr = write_grid_shapefile(fileshp, ml.dis.sr, wel_dict, selectedCell = welCells)
            #
            if self.chkSqlite.isChecked():
                # temporarly disable the pop-up asking CRS
                settings = QSettings()
                oldProjValue = settings.value( "/Projections/defaultBehaviour", "prompt")
                settings.setValue( "/Projections/defaultBehaviour", "useProject" )
                layerName = 'model_well'
                wb = QgsVectorLayer(fileshp, layerName, 'ogr')
                srs = QgsCoordinateReferenceSystem(unicode(crs))
                wb.setCrs(srs)
                srid = wb.crs().postgisSrid()
                uploadQgisVectorLayer(dbname, wb, layerName, srid = srid)
                #
                # Retrieve the Spatialite layer and add it to mapp
                uri = QgsDataSourceURI()
                uri.setDatabase(dbname)
                schema = ''
                table = layerName
                geom_column = "Geometry"
                uri.setDataSource(schema, table, geom_column)
                display_name = table
                wlayer = QgsVectorLayer(uri.uri(), display_name, 'spatialite')
                QgsMapLayerRegistry.instance().addMapLayer(wlayer)
                # revert settings to original
                settings.setValue( "/Projections/defaultBehaviour", oldProjValue )
            #
            #
            val = self.progressBar.value()
            self.progressBar.setValue(val + 2)

        # CHD
        if 'CHD' in pklist:
            fileshp = os.path.join(model_ws, 'model_shapefiles'  ,  'model_chd.shp')
            chd_dict = {}
            chdCells = np.zeros(shape = (nrow, ncol))
            fromLayarray = np.ones(shape = (nrow,ncol), dtype = np.int)
            toLayarray = np.ones(shape = (nrow,ncol), dtype = np.int)
            for t in range(nper):
                nchd = len(ml.chd.stress_period_data[t])
                sheadarray = np.zeros(shape = (nrow,ncol))
                eheadarray = np.zeros(shape = (nrow,ncol))
                for i in range(len(ml.chd.stress_period_data[t])):
                    (lay, r, c, starting, ending )= ml.chd.stress_period_data[t][i]
                    if fromLayarray[r,c] > lay + 1:
                        fromLayarray[r,c] = np.int(lay + 1)
                    if toLayarray[r,c] < lay + 1:
                        toLayarray[r,c] = np.int(lay + 1)
                    sheadarray[r,c] = starting
                    eheadarray[r,c] = ending
                    chdCells[r,c] = 1
                chd_dict['%i_shead'%(t+1)] = sheadarray
                chd_dict['%i_ehead'%(t+1)] = eheadarray
            chd_dict['from_lay'] = fromLayarray
            chd_dict['to_lay'] = toLayarray
            wr = write_grid_shapefile(fileshp, ml.dis.sr, chd_dict, selectedCell = chdCells)
            #
            #
            if self.chkSqlite.isChecked():
                # temporarly disable the pop-up asking CRS
                settings = QSettings()
                oldProjValue = settings.value( "/Projections/defaultBehaviour", "prompt")
                settings.setValue( "/Projections/defaultBehaviour", "useProject" )
                layerName = 'model_chd'
                wb = QgsVectorLayer(fileshp, layerName, 'ogr')
                srs = QgsCoordinateReferenceSystem(unicode(crs))
                wb.setCrs(srs)
                srid = wb.crs().postgisSrid()
                uploadQgisVectorLayer(dbname, wb, layerName, srid = srid)
                #
                # Retrieve the Spatialite layer and add it to mapp
                uri = QgsDataSourceURI()
                uri.setDatabase(dbname)
                schema = ''
                table = layerName
                geom_column = "Geometry"
                uri.setDataSource(schema, table, geom_column)
                display_name = table
                wlayer = QgsVectorLayer(uri.uri(), display_name, 'spatialite')
                QgsMapLayerRegistry.instance().addMapLayer(wlayer)
                # revert settings to original
                settings.setValue( "/Projections/defaultBehaviour", oldProjValue )
            #
            #
            val = self.progressBar.value()
            self.progressBar.setValue(val + 2)

        # GHB
        if 'GHB' in pklist:
            fileshp = os.path.join(model_ws, 'model_shapefiles'  ,  'model_ghb.shp')
            ghb_dict = {}
            ghbCells = np.zeros(shape = (nrow, ncol))
            fromLayarray = np.ones(shape = (nrow,ncol), dtype = np.int)
            toLayarray = np.ones(shape = (nrow,ncol), dtype = np.int)
            for t in range(nper):
                nghb = len(ml.ghb.stress_period_data[t])
                bheadarray = np.zeros(shape = (nrow,ncol))
                condarray = np.zeros(shape = (nrow,ncol))
                xyzarray = np.zeros(shape = (nrow,ncol))
                # first deal with the possibility to manage the existance of (optional) [xyz] input in stress period data
                if len(ml.ghb.stress_period_data[0][0]) == 6:
                    xyzOption = True
                else:
                    xyzOption = False

                for i in range(len(ml.ghb.stress_period_data[t])):
                    if xyzOption:
                        (lay, r, c, bhead, cond , xyz ) = ml.ghb.stress_period_data[t][i]
                    else:
                        (lay, r, c, bhead, cond ) = ml.ghb.stress_period_data[t][i]
                        xyz = 1

                    if fromLayarray[r,c] > lay + 1:
                        fromLayarray[r,c] = np.int(lay + 1)
                    if toLayarray[r,c] < lay + 1:
                        toLayarray[r,c] = np.int(lay + 1)
                    bheadarray[r,c] = bhead
                    condarray[r,c] = cond
                    ghbCells[r,c] = 1
                    xyzarray[r,c] = xyz
                ghb_dict['bhead_%i'%(t+1)] = bheadarray
                ghb_dict['cond_%i'%(t+1)] = condarray
            ghb_dict['from_lay'] = fromLayarray
            ghb_dict['to_lay'] = toLayarray

            wr = write_grid_shapefile(fileshp, ml.dis.sr, ghb_dict, selectedCell = ghbCells)
            #
            if self.chkSqlite.isChecked():
                # temporarly disable the pop-up asking CRS
                settings = QSettings()
                oldProjValue = settings.value( "/Projections/defaultBehaviour", "prompt")
                settings.setValue( "/Projections/defaultBehaviour", "useProject" )
                layerName = 'model_ghb'
                wb = QgsVectorLayer(fileshp, layerName, 'ogr')
                srs = QgsCoordinateReferenceSystem(unicode(crs))
                wb.setCrs(srs)
                srid = wb.crs().postgisSrid()
                uploadQgisVectorLayer(dbname, wb, layerName, srid = srid)
                #
                # Retrieve the Spatialite layer and add it to mapp
                uri = QgsDataSourceURI()
                uri.setDatabase(dbname)
                schema = ''
                table = layerName
                geom_column = "Geometry"
                uri.setDataSource(schema, table, geom_column)
                display_name = table
                wlayer = QgsVectorLayer(uri.uri(), display_name, 'spatialite')
                QgsMapLayerRegistry.instance().addMapLayer(wlayer)
                # revert settings to original
                settings.setValue( "/Projections/defaultBehaviour", oldProjValue )
            #
            val = self.progressBar.value()
            self.progressBar.setValue(val + 2)

        # RIV
        if 'RIV' in pklist:
            fileshp = os.path.join(model_ws, 'model_shapefiles'  ,  'model_riv.shp')
            riv_dict = {}
            rivCells = np.zeros(shape = (nrow, ncol))
            layArray = np.ones(shape = (nrow,ncol), dtype = np.int)
            for t in range(nper):
                nriv = len(ml.riv.stress_period_data[t])
                stagearray = np.zeros(shape = (nrow,ncol))
                rbotarray = np.zeros(shape = (nrow,ncol))
                condarray = np.zeros(shape = (nrow,ncol))
                #
                for i in range(len(ml.riv.stress_period_data[t])):
                    (lay, r, c, stage, cond , rbot ) = ml.riv.stress_period_data[t][i]
                    if lay +1 != 1:
                        layArray[r,c] = np.int(lay + 1)
                    stagearray[r,c] = stage
                    condarray[r,c] = cond
                    rbotarray[r,c] = rbot
                    rivCells[r,c] = 1
                riv_dict['stage_%i'%(t+1)] = stagearray
                riv_dict['cond_%i'%(t+1)] = condarray
                riv_dict['rbot_%i'%(t+1)] = rbotarray

            riv_dict['layer'] = layArray
            wr = write_grid_shapefile(fileshp, ml.dis.sr, riv_dict, selectedCell = rivCells)
            #
            if self.chkSqlite.isChecked():
                # temporarly disable the pop-up asking CRS
                settings = QSettings()
                oldProjValue = settings.value( "/Projections/defaultBehaviour", "prompt")
                settings.setValue( "/Projections/defaultBehaviour", "useProject" )
                layerName = 'model_riv'
                wb = QgsVectorLayer(fileshp, layerName, 'ogr')
                srs = QgsCoordinateReferenceSystem(unicode(crs))
                wb.setCrs(srs)
                srid = wb.crs().postgisSrid()
                uploadQgisVectorLayer(dbname, wb, layerName, srid = srid)
                #
                # Retrieve the Spatialite layer and add it to mapp
                uri = QgsDataSourceURI()
                uri.setDatabase(dbname)
                schema = ''
                table = layerName
                geom_column = "Geometry"
                uri.setDataSource(schema, table, geom_column)
                display_name = table
                wlayer = QgsVectorLayer(uri.uri(), display_name, 'spatialite')
                QgsMapLayerRegistry.instance().addMapLayer(wlayer)
                # revert settings to original
                settings.setValue( "/Projections/defaultBehaviour", oldProjValue )
            #
            val = self.progressBar.value()
            self.progressBar.setValue(val + 2)

        # DRN
        if 'DRN' in pklist:
            fileshp = os.path.join(model_ws, 'model_shapefiles'  ,  'model_drn.shp')
            drn_dict = {}
            drnCells = np.zeros(shape = (nrow, ncol))
            layArray = np.ones(shape = (nrow,ncol), dtype = np.int)
            lengthArray = np.zeros(shape = (nrow,ncol))

            for t in range(nper):
                ndrn = len(ml.drn.stress_period_data[t])
                elevarray = np.zeros(shape = (nrow,ncol))
                condarray = np.zeros(shape = (nrow,ncol))
                #
                for i in range(len(ml.drn.stress_period_data[t])):
                    (lay, r, c, elev, cond ) = ml.drn.stress_period_data[t][i]
                    if lay +1 != 1:
                        layArray[r,c] = np.int(lay + 1)
                    elevarray[r,c] = elev
                    condarray[r,c] = cond
                    drnCells[r,c] = 1
                drn_dict['elev_%i'%(t+1)] = elevarray
                drn_dict['cond_%i'%(t+1)] = condarray

            drn_dict['layer'] = layArray
            drn_dict['length'] = lengthArray
            drn_dict['segment'] = np.ones(shape = (nrow,ncol))
            wr = write_grid_shapefile(fileshp, ml.dis.sr, drn_dict, selectedCell = drnCells)
            #
            if self.chkSqlite.isChecked():
                # temporarly disable the pop-up asking CRS
                settings = QSettings()
                oldProjValue = settings.value( "/Projections/defaultBehaviour", "prompt")
                settings.setValue( "/Projections/defaultBehaviour", "useProject" )
                layerName = 'model_drn'
                wb = QgsVectorLayer(fileshp, layerName, 'ogr')
                srs = QgsCoordinateReferenceSystem(unicode(crs))
                wb.setCrs(srs)
                srid = wb.crs().postgisSrid()
                uploadQgisVectorLayer(dbname, wb, layerName, srid = srid)
                #
                # Retrieve the Spatialite layer and add it to mapp
                uri = QgsDataSourceURI()
                uri.setDatabase(dbname)
                schema = ''
                table = layerName
                geom_column = "Geometry"
                uri.setDataSource(schema, table, geom_column)
                display_name = table
                wlayer = QgsVectorLayer(uri.uri(), display_name, 'spatialite')
                QgsMapLayerRegistry.instance().addMapLayer(wlayer)
                # revert settings to original
                settings.setValue( "/Projections/defaultBehaviour", oldProjValue )
            #
            val = self.progressBar.value()
            self.progressBar.setValue(val + 2)

        # SWI
        print('lista pacchetti', pklist)
        if 'SWI2' in pklist:
            fileshp = os.path.join(model_ws, 'model_shapefiles'  ,  'model_swi.shp')
            swi_dict = {}
            # here we assume that nsfr = 1 (only 1 surface) so that zeta is a listof length1
            # zeta[0] is an array (nlay, nrow, ncol)
            zetaarray = ml.swi2.zeta[0].array
            for l in range(nlay):
                swi_dict['zeta_%i'%(l+1)] = zetaarray[l]
                swi_dict['ssz_%i'%(l+1)] = ml.swi2.ssz.array[l]
                swi_dict['isource_%i'%(l+1)] = ml.swi2.isource.array[l]

            wr = write_grid_shapefile(fileshp, ml.dis.sr, swi_dict)
            #
            if self.chkSqlite.isChecked():
                # temporarly disable the pop-up asking CRS
                settings = QSettings()
                oldProjValue = settings.value( "/Projections/defaultBehaviour", "prompt")
                settings.setValue( "/Projections/defaultBehaviour", "useProject" )
                layerName = 'model_swi'
                wb = QgsVectorLayer(fileshp, layerName, 'ogr')
                srs = QgsCoordinateReferenceSystem(unicode(crs))
                wb.setCrs(srs)
                srid = wb.crs().postgisSrid()
                uploadQgisVectorLayer(dbname, wb, layerName, srid = srid)
                #
                # Retrieve the Spatialite layer and add it to mapp
                uri = QgsDataSourceURI()
                uri.setDatabase(dbname)
                schema = ''
                table = layerName
                geom_column = "Geometry"
                uri.setDataSource(schema, table, geom_column)
                display_name = table
                wlayer = QgsVectorLayer(uri.uri(), display_name, 'spatialite')
                QgsMapLayerRegistry.instance().addMapLayer(wlayer)
                # revert settings to original
                settings.setValue( "/Projections/defaultBehaviour", oldProjValue )
            #
            val = self.progressBar.value()
            self.progressBar.setValue(val + 2)

        self.progressBar.setValue(100)

        messaggio = 'MODFLOW model coorected imported in ' + model_ws
        Logger.information(None, 'Information', '%s' % (messaggio))

