Update(plotting): update plotting examples and fixes for mf6 bugs (#755)

* Update(#741): Remove redundancy from plotting Notebooks

* Removed flopy3_MapExample.ipynb, flopy3_PlotCrossSection_demo.ipynb, flopy3_PlotMapView_demo.ipynb, flopy3_mf6_BasicPlotting.ipynb, flopy3_mf6_plotting_freyberg.ipynb, flopy3_mf6_vertex_plotting.ipynb, flopy3_UnstructuredGridPlotting.ipynb

* Added flopy3.3_PlotCrossSection.ipynb, flopy3.3_PlotMapView.ipynb

* Fixed VertexGrid errors in PlotMapView and PlotCrossSection

* Fix(moflowpy#753): mask no flow and dry cells in PlotUtilities._plot_array_helper for MODFLOW 6 output

*update mfmodel.py hdry and hnoflo to default values

* Cast integer arrays as floating point arrays in PlotUtilities._plot_array_helper

examples/Notebooks/flopy3_PlotArrayExample.ipynb

@@ -4617,7 +4617,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.3"
+   "version": "3.6.5"
   }
  },
  "nbformat": 4,

examples/Notebooks/setup_pmv_demo.py

@@ -6,64 +6,7 @@
 
 
 def run():
-    # run installed version of flopy or add local path
-    try:
-        import flopy
-    except:
-        fpth = os.path.abspath(os.path.join('..', '..'))
-        sys.path.append(fpth)
-        import flopy
-
-    # Set name of MODFLOW exe
-    #  assumes executable is in users path statement
-    version = 'mf2005'
-    exe_name = 'mf2005'
-    exe_mp = 'mp6'
-    if platform.system() == 'Windows':
-        exe_name += '.exe'
-        exe_mp += '.exe'
-    mfexe = exe_name
-
-    # Set the paths
-    loadpth = os.path.join('..', 'data', 'freyberg')
-    modelpth = os.path.join('data')
-
-    # make sure modelpth directory exists
-    if not os.path.exists(modelpth):
-        os.makedirs(modelpth)
-
-    ml = flopy.modflow.Modflow.load('freyberg.nam', model_ws=loadpth,
-                                    exe_name=exe_name, version=version)
-    ml.change_model_ws(new_pth=modelpth)
-    ml.write_input()
-    success, buff = ml.run_model()
-    if not success:
-        print('Something bad happened.')
-    files = ['freyberg.hds', 'freyberg.cbc']
-    for f in files:
-        if os.path.isfile(os.path.join(modelpth, f)):
-            msg = 'Output file located: {}'.format(f)
-            print(msg)
-        else:
-            errmsg = 'Error. Output file cannot be found: {}'.format(f)
-            print(errmsg)
-
-    mp = flopy.modpath.Modpath('freybergmp', exe_name=exe_mp, modflowmodel=ml, model_ws=modelpth)
-    mpbas = flopy.modpath.ModpathBas(mp, hnoflo=ml.bas6.hnoflo, hdry=ml.lpf.hdry,
-                                     ibound=ml.bas6.ibound.array, prsity=0.2, prsityCB=0.2)
-    sim = mp.create_mpsim(trackdir='forward', simtype='endpoint', packages='RCH')
-    mp.write_input()
-    mp.run_model()
-
-    mpp = flopy.modpath.Modpath('freybergmpp', exe_name=exe_mp, modflowmodel=ml, model_ws=modelpth)
-    mpbas = flopy.modpath.ModpathBas(mpp, hnoflo=ml.bas6.hnoflo, hdry=ml.lpf.hdry,
-                                     ibound=ml.bas6.ibound.array, prsity=0.2, prsityCB=0.2)
-    sim = mpp.create_mpsim(trackdir='backward', simtype='pathline', packages='WEL')
-    mpp.write_input()
-    mpp.run_model()
-
-
-    ## load and run second example
+    ## load and run vertex grid example
     # run installed version of flopy or add local path
     try:
         import flopy
@@ -170,6 +113,7 @@ def run():
 
     # node property flow
     npf = flopy.mf6.ModflowGwfnpf(gwf, xt3doptions=[('xt3d')],
+                                  save_specific_discharge=True,
                                   icelltype=[1, 0, 0],
                                   k=[50.0, 0.01, 200.0],
                                   k33=[10., 0.01, 20.])

flopy/mf6/mfmodel.py

@@ -442,11 +442,11 @@ def laytyp(self):
 
     @property
     def hdry(self):
-        return None
+        return -1e30
 
     @property
     def hnoflo(self):
-        return None
+        return 1e30
 
     @property
     def laycbd(self):

flopy/plot/map.py

@@ -128,6 +128,9 @@ def plot_array(self, a, masked_values=None, **kwargs):
                 if a.shape[0] == 1:
                     a = np.squeeze(a, axis=1)
                     plotarray = a[self.layer, :]
+                elif a.shape[1] == 1:
+                    a = np.squeeze(a, axis=1)
+                    plotarray = a[self.layer, :]
                 else:
                     raise Exception("Array must be of dimension 1 or 2")
             elif a.ndim == 2:
@@ -242,6 +245,9 @@ def contour_array(self, a, masked_values=None, **kwargs):
                 if a.shape[0] == 1:
                     a = np.squeeze(a, axis=1)
                     plotarray = a[self.layer, :]
+                elif a.shape[1] == 1:
+                    a = np.squeeze(a, axis=1)
+                    plotarray = a[self.layer, :]
                 else:
                     raise Exception("Array must be of dimension 1 or 2")
             elif a.ndim == 2:

flopy/plot/plotbase.py

@@ -371,10 +371,10 @@ def plot_bc(self, ftype=None, package=None, kper=0, color=None,
         # Plot the list locations, change this to self.mg.shape
         if self.mg.grid_type == "vertex":
             plotarray = np.zeros((self.mg.nlay, self.mg.ncpl), dtype=np.int)
+            plotarray[idx] = 1
         else:
             plotarray = np.zeros((self.mg.nlay, self.mg.nrow, self.mg.ncol), dtype=np.int)
-
-        plotarray[idx] = 1
+            plotarray[idx[0], idx[1], idx[2]] = 1
 
         plotarray = np.ma.masked_equal(plotarray, 0)
         if color is None:
@@ -470,7 +470,7 @@ def plot_specific_discharge(self, spdis, head=None, kstep=1,
         if self.direction == 'x':
             qx[idx] = spdis['qx']
         elif self.direction == 'y':
-            qx[idx] = spdis['qy']
+            qx[idx] = spdis['qy'] * -1
         else:
             err_msg = 'plot_specific_discharge does not ' \
                       'support arbitrary cross-sections'
@@ -663,7 +663,7 @@ def plot_discharge(self, frf, fff, flf=None,
 
             ib = ib.ravel()
             qx = qx.ravel()
-            qy = qy.ravel()
+            qy = qy.ravel() * -1
             qz = qz.ravel()
 
             temp = []

flopy/plot/plotutil.py

@@ -1198,6 +1198,8 @@ def _plot_array_helper(plotarray, model=None, modelgrid=None, axes=None,
             if key in kwargs:
                 defaults[key] = kwargs.pop(key)
 
+        plotarray = plotarray.astype(float)
+
         # test if this is vertex or structured grid
         if model is not None:
             grid_type = model.modelgrid.grid_type
@@ -1223,6 +1225,18 @@ def _plot_array_helper(plotarray, model=None, modelgrid=None, axes=None,
         else:
             grid_type = "structured"
 
+        ib = None
+        if modelgrid is not None:
+            if modelgrid.idomain is not None:
+                ib = modelgrid.idomain
+
+        else:
+            if ib is None:
+                try:
+                    ib = model.modelgrid.idomain
+                except:
+                    pass
+
         # reshape 2d arrays to 3d for convenience
         if len(plotarray.shape) == 2 and grid_type == "structured":
             plotarray = plotarray.reshape((1, plotarray.shape[0],
@@ -1241,7 +1255,7 @@ def _plot_array_helper(plotarray, model=None, modelgrid=None, axes=None,
         for idx, k in enumerate(range(i0, i1)):
             fig = plt.figure(num=fignum[idx])
             pmv = PlotMapView(ax=axes[idx], model=model,
-                             modelgrid=modelgrid, layer=k)
+                              modelgrid=modelgrid, layer=k)
             if defaults['pcolor']:
                 cm = pmv.plot_array(plotarray[k],
                                     masked_values=defaults['masked_values'],
@@ -1266,19 +1280,9 @@ def _plot_array_helper(plotarray, model=None, modelgrid=None, axes=None,
             if defaults['grid']:
                 pmv.plot_grid(ax=axes[idx])
 
-            ib = None
             if defaults['inactive']:
-                if modelgrid is not None:
-                    if modelgrid.idomain is not None:
-                        ib = modelgrid.idomain
-                        pmv.plot_inactive(ibound=ib, ax=axes[idx])
-
-                if ib is None:
-                    try:
-                        ib = model.modelgrid.idomain
-                        pmv.plot_inactive(ibound=ib, ax=axes[idx])
-                    except:
-                        pass
+                if ib is not None:
+                    pmv.plot_inactive(ibound=ib, ax=axes[idx])
 
         if len(axes) == 1:
             axes = axes[0]
@@ -1827,6 +1831,9 @@ def line_intersect_grid(ptsin, xgrid, ygrid):
                             continue
                         elif i in vdict[cell]:
                             continue
+                        elif np.isnan(i[0]) or np.isinf(i[0]) \
+                                or np.isinf(i[1]) or np.isnan(i[1]):
+                            continue
                         else:
                             vdict[cell].append(i)
                 else:
@@ -1840,6 +1847,9 @@ def line_intersect_grid(ptsin, xgrid, ygrid):
                             continue
                         elif i in t:
                             continue
+                        elif np.isnan(i[0]) or np.isinf(i[0]) \
+                                or np.isinf(i[1]) or np.isnan(i[1]):
+                            continue
                         else:
                             t.append(i)