Fix quasi3d plotting (#819)

* Fix plotting methods with quasi3d layers

* Add test

* Rename test

* close figures in test

* fix codacy and travis bugs

* Do not check for quasi3d if nlay==1

And fix vtk-export

* make plot_discharge in map also work without head

* Add laycbd as parameter to StructuredGrid

Make sure zcellcentersis calculated right
fix double line in test (minor change)

* fix xyzcellcenters again

autotest/t070_test_quasi3layers.py

@@ -0,0 +1,98 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Feb 17 12:29:35 2020
+
+@author: Artesia
+"""
+
+
+import os
+import numpy as np
+import flopy
+import matplotlib.pyplot as plt
+
+def test_plotting_with_quasi3d_layers():
+    modelname = 'model_mf'
+    model_ws = os.path.join('.', 'temp', 't069a')
+    exe_name = 'mf2005'
+    mf = flopy.modflow.Modflow(modelname, model_ws=model_ws, exe_name=exe_name)
+
+    # Model domain and grid definition
+    Lx = 1000.
+    Ly = 1000.
+    ztop = 0.
+    zbot = -30.
+    nlay = 3
+    nrow = 10
+    ncol = 10
+    delr = Lx / ncol
+    delc = Ly / nrow
+    laycbd = [0]*(nlay)
+    laycbd[0] = 1
+    botm = np.linspace(ztop, zbot, nlay + np.sum(laycbd) + 1)[1:]
+
+    # Create the discretization object
+    flopy.modflow.ModflowDis(mf, nlay, nrow, ncol, delr=delr, delc=delc,
+                             top=ztop, botm=botm, laycbd=laycbd)
+
+    # Variables for the BAS package
+    ibound = np.ones((nlay, nrow, ncol), dtype=np.int32)
+    ibound[:, :, 0] = -1
+    ibound[:, :, -1] = -1
+    strt = np.ones((nlay, nrow, ncol), dtype=np.float32)
+    strt[:, :, 0] = 10.
+    strt[:, :, -1] = 0.
+    flopy.modflow.ModflowBas(mf, ibound=ibound, strt=strt)
+
+    # Add LPF package to the MODFLOW model
+    flopy.modflow.ModflowLpf(mf, hk=10., vka=10., ipakcb=53, vkcb=10)
+
+    # add a well
+    row = int((nrow-1)/2)
+    col = int((ncol-1)/2)
+    spd = {0:[[1, row, col, -1000]]}
+    flopy.modflow.ModflowWel(mf, stress_period_data=spd)
+
+    # Add OC package to the MODFLOW model
+    spd = {(0, 0): ['save head', 'save budget']}
+    flopy.modflow.ModflowOc(mf, stress_period_data=spd, compact=True)
+
+    # Add PCG package to the MODFLOW model
+    flopy.modflow.ModflowPcg(mf)
+
+    # Write the MODFLOW model input files
+    mf.write_input()
+
+    # Run the MODFLOW model
+    success, buff = mf.run_model()
+
+    # read output
+    hf = flopy.utils.HeadFile(os.path.join(mf.model_ws,'{}.hds'.format(mf.name)))
+    head = hf.get_data(totim=1.0)
+    cbb = flopy.utils.CellBudgetFile(os.path.join(mf.model_ws,'{}.cbc'.format(mf.name)))
+    frf = cbb.get_data(text='FLOW RIGHT FACE', totim=1.0)[0]
+    fff = cbb.get_data(text='FLOW FRONT FACE', totim=1.0)[0]
+    flf = cbb.get_data(text='FLOW LOWER FACE', totim=1.0)[0]
+
+    # plot a map
+    plt.figure()
+    mv = flopy.plot.PlotMapView(model=mf,layer=1)
+    mv.plot_grid()
+    mv.plot_array(head)
+    mv.contour_array(head)
+    mv.plot_ibound()
+    mv.plot_bc('wel')
+    mv.plot_discharge(frf,fff, head=head)
+    plt.close()
+
+    # plot a cross-section
+    plt.figure()
+    cs = flopy.plot.PlotCrossSection(model=mf, line={'row':int((nrow-1)/2)})
+    cs.plot_grid()
+    cs.plot_array(head)
+    cs.contour_array(head)
+    cs.plot_ibound()
+    cs.plot_bc('wel')
+    cs.plot_discharge(frf, fff, flf, head=head)
+    plt.close()
+

flopy/discretization/structuredgrid.py

@@ -37,7 +37,8 @@ class for a structured model grid
     """
     def __init__(self, delc=None, delr=None, top=None, botm=None, idomain=None,
                  lenuni=None, epsg=None, proj4=None, prj=None, xoff=0.0,
-                 yoff=0.0, angrot=0.0, nlay=None, nrow=None, ncol=None):
+                 yoff=0.0, angrot=0.0, nlay=None, nrow=None, ncol=None,
+                 laycbd=None):
         super(StructuredGrid, self).__init__('structured', top, botm, idomain,
                                              lenuni, epsg, proj4, prj, xoff,
                                              yoff, angrot)
@@ -55,15 +56,19 @@ def __init__(self, delc=None, delr=None, top=None, botm=None, idomain=None,
             assert self.__nrow * self.__ncol == len(np.ravel(top))
         if botm is not None:
             assert self.__nrow * self.__ncol == len(np.ravel(botm[0]))
-            if nlay is not None and nlay < len(botm):
-                self.__nlay_nocbd = nlay
+            if nlay is not None:
+                self.__nlay = nlay
             else:
-                self.__nlay_nocbd = len(botm)
-
-            self.__nlay = len(botm)
+                if laycbd is not None:
+                    self.__nlay = len(botm) - np.sum(laycbd>0)
+                else:
+                    self.__nlay = len(botm)
         else:
             self.__nlay = nlay
-            self.__nlay_nocbd = nlay
+        if laycbd is not None:
+            self.__laycbd = laycbd
+        else:
+            self.__laycbd = np.zeros(self.__nlay, dtype=int)
 
     ####################
     # Properties
@@ -81,10 +86,6 @@ def is_complete(self):
             return True
         return False
 
-    @property
-    def nlay_nocbd(self):
-        return self.__nlay_nocbd
-
     @property
     def nlay(self):
         return self.__nlay
@@ -184,9 +185,13 @@ def xyzcellcenters(self):
                 # get z centers
                 z = np.empty((self.__nlay, self.__nrow, self.__ncol))
                 z[0, :, :] = (self._top[:, :] + self._botm[0, :, :]) / 2.
-                for l in range(1, self.__nlay):
-                    z[l, :, :] = (self._botm[l - 1, :, :] +
-                                  self._botm[l, :, :]) / 2.
+                ibs = np.arange(self.__nlay)
+                quasi3d = [cbd !=0 for cbd in self.__laycbd]
+                if np.any(quasi3d):
+                    ibs[1:] = ibs[1:] + np.cumsum(quasi3d)[:self.__nlay - 1]
+                for l, ib in enumerate(ibs[1:], 1):
+                    z[l, :, :] = (self._botm[ib - 1, :, :] +
+                                  self._botm[ib, :, :]) / 2.
             else:
                 z = None
             if self._has_ref_coordinates:

flopy/export/vtk.py

@@ -203,13 +203,13 @@ def __init__(self, model, verbose=None, nanval=-1e+20, smooth=False,
         self.model = model
         self.modelgrid = model.modelgrid
         self.arrays = {}
-        self.shape = (self.modelgrid.nlay, self.modelgrid.nrow,
-                      self.modelgrid.ncol)
-        self.shape2d = (self.shape[1], self.shape[2])
         self.nlay = self.modelgrid.nlay
+        if hasattr(self.model, 'dis') and hasattr(self.model.dis, 'laycbd'):
+            self.nlay = self.nlay + np.sum(self.model.dis.laycbd.array > 0)
         self.nrow = self.modelgrid.nrow
         self.ncol = self.modelgrid.ncol
-        self.ncol = self.modelgrid.ncol
+        self.shape = (self.nlay, self.nrow, self.ncol)
+        self.shape2d = (self.shape[1], self.shape[2])
         self.nanval = nanval
 
         self.cell_type = 11

flopy/modflow/mf.py

@@ -280,7 +280,8 @@ def modelgrid(self):
                                              xoff=self._modelgrid.xoffset,
                                              yoff=self._modelgrid.yoffset,
                                              angrot=self._modelgrid.angrot,
-                                             nlay=self.dis.nlay)
+                                             nlay=self.dis.nlay,
+                                             laycbd=self.dis.laycbd)
 
         # resolve offsets
         xoff = self._modelgrid.xoffset

flopy/plot/crosssection.py

@@ -298,6 +298,8 @@ def get_centergrids(self, xpts, zpts):
                         break
         else:
             for k in range(0, zpts.shape[0] - 1):
+                if not self.active[k]:
+                    continue
                 nz += 1
                 nx = 0
                 for i in range(0, xpts.shape[0], 2):
@@ -781,6 +783,8 @@ def set_zcentergrid(self, vs):
                     zcentergrid.append(zp)
         else:
             for k in range(0, self.zpts.shape[0] - 1):
+                if not self.active[k]==1:
+                    continue
                 nz += 1
                 nx = 0
                 for i in range(0, self.xpts.shape[0], 2):

flopy/plot/map.py

@@ -843,10 +843,10 @@ def plot_discharge(self, frf=None, fff=None,
             delc = self.mg.delc
             top = np.copy(self.mg.top)
             botm = np.copy(self.mg.botm)
-            nlay = botm.shape[0]
             laytyp = None
             hnoflo = 999.
             hdry = 999.
+            laycbd = None
 
             if self.model is not None:
                 if self.model.laytyp is not None:
@@ -858,11 +858,23 @@ def plot_discharge(self, frf=None, fff=None,
                 if self.model.hdry is not None:
                     hdry = self.model.hdry
 
+                if self.model.laycbd is not None:
+                    laycbd = self.model.laycbd
+
+            if laycbd is not None and 1 in laycbd:
+                active = np.ones((botm.shape[0],), dtype=np.int)
+                kon = 0
+                for cbd in laycbd:
+                    if cbd > 0:
+                        kon += 1
+                        active[kon] = 0
+                botm = botm[active==1]
+
             # If no access to head or laytyp, then calculate confined saturated
             # thickness by setting laytyp to zeros
             if head is None or laytyp is None:
                 head = np.zeros(botm.shape, np.float32)
-                laytyp = np.zeros((nlay,), dtype=np.int)
+                laytyp = np.zeros((botm.shape[0],), dtype=np.int)
 
             # calculate the saturated thickness
             sat_thk = plotutil.PlotUtilities. \

flopy/plot/plotbase.py

@@ -650,6 +650,8 @@ def plot_discharge(self, frf, fff, flf=None,
             delc = self.mg.delc
             top = self.mg.top
             botm = self.mg.botm
+            if not np.all(self.active==1):
+                botm = botm[self.active==1]
             nlay = botm.shape[0]
             laytyp = None
             hnoflo = 999.

flopy/utils/check.py

@@ -385,7 +385,7 @@ def get_active(self, include_cbd=False):
         """
         mg = self.model.modelgrid
         if mg.grid_type == 'structured':
-            inds = (mg.nlay_nocbd, mg.nrow, mg.ncol)
+            inds = (mg.nlay, mg.nrow, mg.ncol)
         elif mg.grid_type == 'vertex':
             inds = (mg.nlay, mg.ncpl)
         else: