Update(plotting and associated modules): major plotting library updates (#1078)

* Removed vcrosssection.py and plotbase.py
* PlotCrossSection() and PlotMapView() updates
* Added plot_shapes() to PlotMapView
* Added styles module for report ready plots using matplotlib contexts
* plotutil.py updates
* postprocessing.py updates
* Fix(plot_vector): set xcenters fixed for kstep parameter
* Update(plotting, reference.py) deprecation warnings updated for planned method removal in version 3.3.5

MANIFEST.in

@@ -1,3 +1,4 @@
 include flopy/export/longnames.json
 include flopy/export/unitsformat.json
 include flopy/mf6/data/dfn/*.dfn
+include flopy/plot/mplstyle/*.mplstyle

autotest/t007_test.py

@@ -1019,18 +1019,18 @@ def test_sr_with_Map():
         modelmap = flopy.plot.ModelMap(model=m, xul=xul, yul=yul,
                                        rotation=rotation)
         assert len(w) == 2, len(w)
-        assert w[0].category == PendingDeprecationWarning, w[0]
-        assert 'ModelMap will be replaced by PlotMapView' in str(w[0].message)
+        assert w[0].category == DeprecationWarning, w[0]
+        assert 'ModelMap is deprecated' in str(w[0].message)
         assert w[1].category == DeprecationWarning, w[1]
         assert 'xul/yul have been deprecated' in str(w[1].message)
 
-    lc = modelmap.plot_grid()
+    pc = modelmap.plot_grid()
     xll, yll = modelmap.mg.xoffset, modelmap.mg.yoffset
     plt.close()
 
     def check_vertices():
-        xllp, yllp = lc._paths[0].vertices[0]
-        xulp, yulp = lc._paths[0].vertices[1]
+        xllp, yllp = pc._paths[780].vertices[3]
+        xulp, yulp = pc._paths[0].vertices[0]
         assert np.abs(xllp - xll) < 1e-6
         assert np.abs(yllp - yll) < 1e-6
         assert np.abs(xulp - xul) < 1e-6
@@ -1044,10 +1044,10 @@ def test_sr_with_Map():
         modelmap = flopy.plot.ModelMap(model=m, xll=xll, yll=yll,
                                        rotation=rotation)
         assert len(w) == 1, len(w)
-        assert w[0].category == PendingDeprecationWarning, w[0]
-        assert 'ModelMap will be replaced by PlotMapView' in str(w[0].message)
+        assert w[0].category == DeprecationWarning, w[0]
+        assert 'ModelMap is deprecated' in str(w[0].message)
 
-    lc = modelmap.plot_grid()
+    pc = modelmap.plot_grid()
     check_vertices()
     plt.close()
 
@@ -1062,10 +1062,10 @@ def test_sr_with_Map():
         modelmap = flopy.plot.ModelMap(model=m)
 
         assert len(w) == 1, len(w)
-        assert w[0].category == PendingDeprecationWarning, w[0]
-        assert 'ModelMap will be replaced by PlotMapView' in str(w[0].message)
+        assert w[0].category == DeprecationWarning, w[0]
+        assert 'ModelMap is deprecated' in str(w[0].message)
 
-    lc = modelmap.plot_grid()
+    pc = modelmap.plot_grid()
     check_vertices()
     plt.close()
 
@@ -1078,10 +1078,10 @@ def test_sr_with_Map():
         modelmap = flopy.plot.ModelMap(model=m, sr=sr)
 
         assert len(w) == 1, len(w)
-        assert w[0].category == PendingDeprecationWarning, w[0]
-        assert 'ModelMap will be replaced by PlotMapView' in str(w[0].message)
+        assert w[0].category == DeprecationWarning, w[0]
+        assert 'ModelMap is deprecated' in str(w[0].message)
 
-    lc = modelmap.plot_grid()
+    pc = modelmap.plot_grid()
     check_vertices()
     plt.close()
 
@@ -1111,14 +1111,14 @@ def test_sr_with_Map():
             assert 'SpatialReference has been deprecated' in str(w[0].message)
             assert w[1].category == DeprecationWarning, w[1]
             assert 'SpatialReference has been deprecated' in str(w[1].message)
-        assert w[-3].category == PendingDeprecationWarning, w[-3]
-        assert 'ModelCrossSection will be replaced by' in str(w[-3].message)
+        assert w[-3].category == DeprecationWarning, w[-3]
+        assert 'ModelCrossSection is Deprecated' in str(w[-3].message)
         assert w[-2].category == DeprecationWarning, w[-2]
         assert 'xul/yul have been deprecated' in str(w[-2].message)
         assert w[-1].category == DeprecationWarning, w[-1]
         assert 'xul/yul have been deprecated' in str(w[-1].message)
 
-    linecollection = modelxsect.plot_grid()
+    patchcollection = modelxsect.plot_grid()
     plt.close()
 
 
@@ -1132,23 +1132,18 @@ def test_modelgrid_with_PlotMapView():
     xll, yll, rotation = 500000., 2934000., 45.
 
     def check_vertices():
-        # vertices = modelmap.mg.xyvertices
-        xllp, yllp = lc._paths[0].vertices[0]
-        # xulp, yulp = lc._paths[0].vertices[1]
+        xllp, yllp = pc._paths[780].vertices[3]
         assert np.abs(xllp - xll) < 1e-6
         assert np.abs(yllp - yll) < 1e-6
-        # assert np.abs(xulp - xul) < 1e-6
-        # assert np.abs(yulp - yul) < 1e-6
 
-    #    check_vertices()
     m.modelgrid.set_coord_info(xoff=xll, yoff=yll, angrot=rotation)
     modelmap = flopy.plot.PlotMapView(model=m)
-    lc = modelmap.plot_grid()
+    pc = modelmap.plot_grid()
     check_vertices()
     plt.close()
 
     modelmap = flopy.plot.PlotMapView(modelgrid=m.modelgrid)
-    lc = modelmap.plot_grid()
+    pc = modelmap.plot_grid()
     check_vertices()
     plt.close()
 
@@ -1157,13 +1152,11 @@ def test_modelgrid_with_PlotMapView():
     # Model domain and grid definition
     dis = flopy.modflow.ModflowDis(mf, nlay=1, nrow=10, ncol=20, delr=1.,
                                    delc=1., xul=100, yul=210)
-    # fig, ax = plt.subplots()
+
     verts = [[101., 201.], [119., 209.]]
-    # modelxsect = flopy.plot.ModelCrossSection(model=mf, line={'line': verts},
-    #                                           xul=mf.dis.sr.xul, yul=mf.dis.sr.yul)
     mf.modelgrid.set_coord_info(xoff=mf.dis.sr.xll, yoff=mf.dis.sr.yll)
     modelxsect = flopy.plot.PlotCrossSection(model=mf, line={'line': verts})
-    linecollection = modelxsect.plot_grid()
+    patchcollection = modelxsect.plot_grid()
     plt.close()
 
 
@@ -1206,7 +1199,7 @@ def test_mapview_plot_bc():
         raise AssertionError("Boundary condition was not drawn")
 
     for col in ax.collections:
-        if not isinstance(col, QuadMesh):
+        if not isinstance(col, PatchCollection):
             raise AssertionError("Unexpected collection type")
     plt.close()
 
@@ -1231,7 +1224,7 @@ def test_mapview_plot_bc():
         raise AssertionError("Boundary condition was not drawn")
 
     for col in ax.collections:
-        if not isinstance(col, QuadMesh):
+        if not isinstance(col, PatchCollection):
             raise AssertionError("Unexpected collection type")
     plt.close()
 
@@ -1249,7 +1242,7 @@ def test_mapview_plot_bc():
         raise AssertionError("Boundary condition was not drawn")
 
     for col in ax.collections:
-        if not isinstance(col, QuadMesh):
+        if not isinstance(col, PatchCollection):
             raise AssertionError("Unexpected collection type")
     plt.close()
 
@@ -1421,6 +1414,7 @@ def test_get_lrc_get_node():
         assert node == n, "get_node() returned {}, expecting {}".format(n, node)
     return
 
+
 def test_vertex_model_dot_plot():
     import matplotlib.pyplot as plt
     from matplotlib import rcParams

autotest/t042_test.py

@@ -93,7 +93,7 @@ def test_get_sat_thickness_gradients():
     assert np.nansum(np.abs(dh / dz - grad[:, 1, 0])) < 1e-6
 
     sat_thick = get_saturated_thickness(hds, m, nodata)
-    assert np.abs(np.sum(sat_thick[:, 1, 1] - np.array([0.2, 1., 1.]))) < 1e-6
+    assert np.abs(np.sum(sat_thick[:, 1, 1] - np.array([0.2, 1., 1.6]))) < 1e-6
 
 if __name__ == '__main__':
     #test_get_transmissivities()

autotest/t050_test.py

@@ -355,15 +355,21 @@ def test_vtk_cbc():
 
 def test_vtk_vector():
     from flopy.utils import postprocessing as pp
+    from flopy.utils import HeadFile, CellBudgetFile
     # test mf 2005 freyberg
     mpth = os.path.join('..', 'examples', 'data',
                         'freyberg_multilayer_transient')
     namfile = 'freyberg.nam'
     cbcfile = os.path.join(mpth, 'freyberg.cbc')
     hdsfile = os.path.join(mpth, 'freyberg.hds')
+    cbc = CellBudgetFile(cbcfile)
+    keys = ["FLOW RIGHT FACE", "FLOW FRONT FACE", "FLOW LOWER FACE"]
+    vectors = [cbc.get_data(text=t)[0] for t in keys]
+    hds = HeadFile(hdsfile)
+    head = hds.get_data()
     m = flopy.modflow.Modflow.load(namfile, model_ws=mpth, verbose=False,
                                    load_only=['dis', 'bas6', 'upw'])
-    q = pp.get_specific_discharge(m, cbcfile=cbcfile)
+    q = pp.get_specific_discharge(vectors, m, head)
     output_dir = os.path.join(cpth, 'freyberg_vector')
     filenametocheck = 'discharge.vtu'
 
@@ -386,12 +392,12 @@ def test_vtk_vector():
     assert(os.path.exists(filetocheck))
 
     # with values directly given at vertices
-    q = pp.get_specific_discharge(m, cbcfile=cbcfile, hdsfile=hdsfile,
+    q = pp.get_specific_discharge(vectors, m, head,
                                   position='vertices')
     nancount = np.count_nonzero(np.isnan(q[0]))
-    assert(nancount==308)
+    assert(nancount == 472)
     overall = np.nansum(q[0]) + np.nansum(q[1]) + np.nansum(q[2])
-    assert np.allclose(overall, -15.467904755216372)
+    assert np.allclose(overall, -45.38671967357735)
     output_dir = os.path.join(cpth, 'freyberg_vector')
     filenametocheck = 'discharge_verts.vtu'
     vtk.export_vector(m, q, output_dir, 'discharge_verts')
@@ -655,7 +661,10 @@ def test_vtk_export_true2d_nonregxy():
 
     # export and check specific discharge given at vertices
     cbcfile = os.path.join(output_dir, name + '.cbc')
-    q = pp.get_specific_discharge(m, cbcfile, position='vertices')
+    cbc = bf.CellBudgetFile(cbcfile)
+    keys = ["FLOW RIGHT FACE", "FLOW FRONT FACE"]
+    vectors = [cbc.get_data(text=t)[0] for t in keys]
+    q = pp.get_specific_discharge(vectors, m, position='vertices')
     vtk.export_vector(m, q, output_dir, name + '_q', point_scalars=True,
                       true2d=True)
     filetocheck = os.path.join(output_dir, name + '_q.vtr')
@@ -663,9 +672,9 @@ def test_vtk_export_true2d_nonregxy():
     # assert(totalbytes1==5772)
     nlines1 = count_lines_in_file(filetocheck)
     assert(nlines1==54)
-
     return
 
+
 def test_vtk_export_true2d_nonregxz():
     import flopy.utils.binaryfile as bf
     from flopy.utils import postprocessing as pp
@@ -722,7 +731,11 @@ def test_vtk_export_true2d_nonregxz():
 
     # export and check specific discharge given at vertices
     cbcfile = os.path.join(output_dir, name + '.cbc')
-    q = pp.get_specific_discharge(m, cbcfile, position='vertices')
+    cbc = bf.CellBudgetFile(cbcfile)
+    keys = ["FLOW RIGHT FACE", "FLOW LOWER FACE"]
+    vectors = [cbc.get_data(text=t)[0] for t in keys]
+    vectors.insert(1, None)
+    q = pp.get_specific_discharge(vectors, m, position='vertices')
     vtk.export_vector(m, q, output_dir, name + '_q', point_scalars=True,
                       true2d=True)
     filetocheck = os.path.join(output_dir, name + '_q.vtu')
@@ -730,9 +743,9 @@ def test_vtk_export_true2d_nonregxz():
     # assert(totalbytes2==7036)
     nlines2 = count_lines_in_file(filetocheck)
     assert(nlines2==123)
-
     return
 
+
 def test_vtk_export_true2d_nonregyz():
     import flopy.utils.binaryfile as bf
     from flopy.utils import postprocessing as pp
@@ -789,7 +802,11 @@ def test_vtk_export_true2d_nonregyz():
 
     # export and check specific discharge given at vertices
     cbcfile = os.path.join(output_dir, name + '.cbc')
-    q = pp.get_specific_discharge(m, cbcfile, position='vertices')
+    cbc = bf.CellBudgetFile(cbcfile)
+    keys = ["FLOW FRONT FACE", "FLOW LOWER FACE"]
+    vectors = [cbc.get_data(text=t)[0] for t in keys]
+    vectors.insert(0, None)
+    q = pp.get_specific_discharge(vectors, m, position='vertices')
     vtk.export_vector(m, q, output_dir, name + '_q', point_scalars=True,
                       true2d=True)
     filetocheck = os.path.join(output_dir, name + '_q.vtu')
@@ -797,9 +814,9 @@ def test_vtk_export_true2d_nonregyz():
     # assert(totalbytes2==7032)
     nlines2 = count_lines_in_file(filetocheck)
     assert(nlines2==123)
-
     return
 
+
 if __name__ == '__main__':
     test_vtk_export_array2d()
     test_vtk_export_array3d()

autotest/t065_test_gridintersect.py

@@ -1265,8 +1265,7 @@ def test_rasters():
     if (np.max(data) - 2608.557) > 1e-4:
         raise AssertionError
 
-    data = rio.resample_to_grid(ml.modelgrid.xcellcenters,
-                                ml.modelgrid.ycellcenters,
+    data = rio.resample_to_grid(ml.modelgrid,
                                 band=rio.bands[0],
                                 method="nearest")
     if data.size != 5913:

autotest/t072_test_spedis.py

@@ -333,8 +333,10 @@ def test_extended_budget_comprehensive():
 def test_specific_discharge_default():
     # load and postprocess
     mf = flopy.modflow.Modflow.load(namfile_mf2005, check=False)
-    qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(mf,
-                                            cbcfile_mf2005)
+    cbc = bf.CellBudgetFile(cbcfile_mf2005)
+    keys = ["FLOW RIGHT FACE", "FLOW FRONT FACE", "FLOW LOWER FACE"]
+    vectors = [cbc.get_data(text=t)[0] for t in keys]
+    qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(vectors, mf)
 
     # overall check
     overall = np.sum(qx) + np.sum(qy) + np.sum(qz)
@@ -345,20 +347,31 @@ def test_specific_discharge_comprehensive():
     import matplotlib.pyplot as plt
     from matplotlib.quiver import Quiver
 
+    hds = bf.HeadFile(hdsfile_mf2005)
+    head = hds.get_data()
     # load and postprocess
     mf = flopy.modflow.Modflow.load(namfile_mf2005, check=False)
-    qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(mf,
-                             cbcfile_mf2005,
-                             boundary_ifaces=boundary_ifaces,
-                             hdsfile=hdsfile_mf2005)
+    Qx_ext, Qy_ext, Qz_ext = \
+        flopy.utils.postprocessing.get_extended_budget(
+            cbcfile_mf2005,
+            boundary_ifaces=boundary_ifaces,
+            hdsfile=hdsfile_mf2005,
+            model=mf
+        )
+
+    qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(
+        (Qx_ext, Qy_ext, Qz_ext),
+        mf,
+        head
+    )
 
     # check nan values
     assert np.isnan(qx[0, 0, 2])
     assert np.isnan(qx[1, 0, 1])
 
     # overall check
-    overall = np.nansum(qx) + np.nansum(qy) + np.nansum(qz)
-    assert np.allclose(overall, -0.8558630187423599)
+    overall = np.nansum(qz)  # np.nansum(qx) + np.nansum(qy) + np.nansum(qz)
+    assert np.allclose(overall, -4.43224582939148)
 
     # plot discharge in map view
     lay = 1
@@ -366,7 +379,6 @@ def test_specific_discharge_comprehensive():
     quiver = modelmap.plot_vector(qx, qy, normalize=True,
                                   masked_values=[qx[lay, 0, 0]],
                                   color='orange')
-
     # check plot
     ax = modelmap.ax
     if len(ax.collections) == 0:
@@ -378,7 +390,7 @@ def test_specific_discharge_comprehensive():
     pos = np.sum(quiver.X) + np.sum(quiver.Y)
     assert np.allclose(pos, 1600.)
     val = np.sum(quiver.U) + np.sum(quiver.V)
-    assert np.allclose(val, 10.908548650065649)
+    assert np.allclose(val, 10.11359908150753)
 
     # close figure
     plt.close()
@@ -402,11 +414,11 @@ def test_specific_discharge_comprehensive():
     X = np.ma.masked_where(quiver.Umask, quiver.X)
     Y = np.ma.masked_where(quiver.Umask, quiver.Y)
     pos = X.sum() + Y.sum()
-    assert np.allclose(pos, -153.8352064441874)
+    assert np.allclose(pos, -153.8341064453125)
     U = np.ma.masked_where(quiver.Umask, quiver.U)
     V = np.ma.masked_where(quiver.Umask, quiver.V)
     val = U.sum() + V.sum()
-    assert np.allclose(val, -3.25453417836753)
+    assert np.allclose(val, -5.2501876516091235)
 
     # close figure
     plt.close()
@@ -424,9 +436,13 @@ def test_specific_discharge_mf6():
     sim = MFSimulation.load(sim_name=modelname_mf6, sim_ws=modelws_mf6,
                             verbosity_level=0)
     gwf = sim.get_model(modelname_mf6)
-    qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(gwf,
-                       cbcfile_mf6, precision='double',
-                       hdsfile=hdsfile_mf6)
+    hds = bf.HeadFile(hdsfile_mf6)
+    head = hds.get_data()
+    cbc = bf.CellBudgetFile(cbcfile_mf6)
+    spdis = cbc.get_data(text="SPDIS")[0]
+    qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(
+        spdis, gwf, head
+        )
 
     # check nan values
     assert np.isnan(qx[0, 0, 2])

examples/common/setup_pmv_demo.py

@@ -3,23 +3,17 @@ import numpy as np
 import os
 import platform
 import sys
+try:
+    import flopy
+except:
+    fpth = os.path.abspath(os.path.join("..", "..", ".."))
+    sys.path.append(fpth)
+    import flopy
 
 
 def run():
     ## load and run vertex grid example
     # 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
-
-        print(sys.version)
-    print("numpy version: {}".format(np.__version__))
-    print("matplotlib version: {}".format(mpl.__version__))
-    print("flopy version: {}".format(flopy.__version__))
-
     if not os.path.exists("data"):
         os.mkdir("data")
 
@@ -328,5 +322,287 @@ def run():
     return
 
 
+example_name = "ex-gwt-keating"
+
+# Model units
+
+length_units = "m"
+time_units = "days"
+
+# Table of model parameters
+
+nlay = 80  # Number of layers
+nrow = 1  # Number of rows
+ncol = 400  # Number of columns
+delr = 25.0  # Column width ($m$)
+delc = 1.0  # Row width ($m$)
+delz = 25.0  # Layer thickness ($m$)
+top = 2000.0  # Top of model domain ($m$)
+bottom = 0.0  # Bottom of model domain ($m$)
+hka = 1.0e-12  # Permeability of aquifer ($m^2$)
+hkc = 1.0e-18  # Permeability of aquitard ($m^2$)
+h1 = 800.0  # Head on left side ($m$)
+h2 = 100.0  # Head on right side ($m$)
+recharge = 0.5  # Recharge ($kg/s$)
+recharge_conc = 1.0  # Normalized recharge concentration (unitless)
+alpha_l = 1.0  # Longitudinal dispersivity ($m$)
+alpha_th = 1.0  # Transverse horizontal dispersivity ($m$)
+alpha_tv = 1.0  # Transverse vertical dispersivity ($m$)
+period1 = 730  # Length of first simulation period ($d$)
+period2 = 29270.0  # Length of second simulation period ($d$)
+porosity = 0.1  # Porosity of mobile domain (unitless)
+obs1 = (49, 1, 119)  # Layer, row, and column for observation 1
+obs2 = (77, 1, 359)  # Layer, row, and column for observation 2
+
+obs1 = tuple([i - 1 for i in obs1])
+obs2 = tuple([i - 1 for i in obs2])
+seconds_to_days = 24.0 * 60.0 * 60.0
+permeability_to_conductivity = 1000.0 * 9.81 / 1.0e-3 * seconds_to_days
+hka = hka * permeability_to_conductivity
+hkc = hkc * permeability_to_conductivity
+botm = [top - (k + 1) * delz for k in range(nlay)]
+x = np.arange(0, 10000.0, delr) + delr / 2.0
+plotaspect = 1.0
+
+# Fill hydraulic conductivity array
+hydraulic_conductivity = np.ones((nlay, nrow, ncol), dtype=float) * hka
+for k in range(nlay):
+    if 1000.0 <= botm[k] < 1100.0:
+        for j in range(ncol):
+            if 3000.0 <= x[j] <= 6000.0:
+                hydraulic_conductivity[k, 0, j] = hkc
+
+# Calculate recharge by converting from kg/s to m/d
+rcol = []
+for jcol in range(ncol):
+    if 4200.0 <= x[jcol] <= 4800.0:
+        rcol.append(jcol)
+number_recharge_cells = len(rcol)
+rrate = recharge * seconds_to_days / 1000.0
+cell_area = delr * delc
+rrate = rrate / (float(number_recharge_cells) * cell_area)
+rchspd = {}
+rchspd[0] = [[(0, 0, j), rrate, recharge_conc] for j in rcol]
+rchspd[1] = [[(0, 0, j), rrate, 0.0] for j in rcol]
+
+def build_mf6gwf(sim_folder):
+    print("Building mf6gwf model...{}".format(sim_folder))
+    ws = os.path.join("data", "mf6-gwt-keating")
+    name = "flow"
+    sim_ws = os.path.join(ws, 'mf6gwf')
+    sim = flopy.mf6.MFSimulation(
+        sim_name=name, sim_ws=sim_ws, exe_name="mf6"
+    )
+    tdis_ds = ((period1, 1, 1.0), (period2, 1, 1.0))
+    flopy.mf6.ModflowTdis(
+        sim, nper=len(tdis_ds), perioddata=tdis_ds, time_units=time_units
+    )
+    flopy.mf6.ModflowIms(
+        sim,
+        print_option="summary",
+        complexity="complex",
+        no_ptcrecord="all",
+        outer_dvclose=1.0e-4,
+        outer_maximum=2000,
+        under_relaxation="dbd",
+        linear_acceleration="BICGSTAB",
+        under_relaxation_theta=0.7,
+        under_relaxation_kappa=0.08,
+        under_relaxation_gamma=0.05,
+        under_relaxation_momentum=0.0,
+        backtracking_number=20,
+        backtracking_tolerance=2.0,
+        backtracking_reduction_factor=0.2,
+        backtracking_residual_limit=5.0e-4,
+        inner_dvclose=1.0e-5,
+        rcloserecord=[0.0001, "relative_rclose"],
+        inner_maximum=100,
+        relaxation_factor=0.0,
+        number_orthogonalizations=2,
+        preconditioner_levels=8,
+        preconditioner_drop_tolerance=0.001,
+    )
+    gwf = flopy.mf6.ModflowGwf(
+        sim, modelname=name, save_flows=True, newtonoptions=["newton"]
+    )
+    flopy.mf6.ModflowGwfdis(
+        gwf,
+        length_units=length_units,
+        nlay=nlay,
+        nrow=nrow,
+        ncol=ncol,
+        delr=delr,
+        delc=delc,
+        top=top,
+        botm=botm,
+    )
+    flopy.mf6.ModflowGwfnpf(
+        gwf,
+        save_specific_discharge=True,
+        save_saturation=True,
+        icelltype=1,
+        k=hydraulic_conductivity,
+    )
+    flopy.mf6.ModflowGwfic(gwf, strt=600.0)
+    chdspd = [[(k, 0, 0), h1] for k in range(nlay) if botm[k] < h1]
+    chdspd += [[(k, 0, ncol - 1), h2] for k in range(nlay) if botm[k] < h2]
+    flopy.mf6.ModflowGwfchd(
+        gwf,
+        stress_period_data=chdspd,
+        print_input=True,
+        print_flows=True,
+        save_flows=False,
+        pname="CHD-1",
+    )
+    flopy.mf6.ModflowGwfrch(
+        gwf,
+        stress_period_data=rchspd,
+        auxiliary=["concentration"],
+        pname="RCH-1",
+    )
+
+    head_filerecord = "{}.hds".format(name)
+    budget_filerecord = "{}.bud".format(name)
+    flopy.mf6.ModflowGwfoc(
+        gwf,
+        head_filerecord=head_filerecord,
+        budget_filerecord=budget_filerecord,
+        saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")],
+    )
+    return sim
+
+
+def build_mf6gwt(sim_folder):
+    print("Building mf6gwt model...{}".format(sim_folder))
+    ws = os.path.join("data", "mf6-gwt-keating")
+    name = "trans"
+    sim_ws = os.path.join(ws, 'mf6gwt')
+    sim = flopy.mf6.MFSimulation(
+        sim_name=name,
+        sim_ws=sim_ws,
+        exe_name="mf6",
+        continue_=True,
+    )
+    tdis_ds = ((period1, 73, 1.0), (period2, 2927, 1.0))
+    flopy.mf6.ModflowTdis(
+        sim, nper=len(tdis_ds), perioddata=tdis_ds, time_units=time_units
+    )
+    flopy.mf6.ModflowIms(
+        sim,
+        print_option="summary",
+        outer_dvclose=1.0e-4,
+        outer_maximum=100,
+        under_relaxation="none",
+        linear_acceleration="BICGSTAB",
+        rcloserecord=[1000.0, "strict"],
+        inner_maximum=20,
+        inner_dvclose=1.0e-4,
+        relaxation_factor=0.0,
+        number_orthogonalizations=2,
+        preconditioner_levels=8,
+        preconditioner_drop_tolerance=0.001,
+    )
+    gwt = flopy.mf6.ModflowGwt(sim, modelname=name, save_flows=True)
+    flopy.mf6.ModflowGwtdis(
+        gwt,
+        length_units=length_units,
+        nlay=nlay,
+        nrow=nrow,
+        ncol=ncol,
+        delr=delr,
+        delc=delc,
+        top=top,
+        botm=botm,
+    )
+    flopy.mf6.ModflowGwtic(gwt, strt=0)
+    flopy.mf6.ModflowGwtmst(gwt, porosity=porosity)
+    flopy.mf6.ModflowGwtadv(gwt, scheme="upstream")
+    flopy.mf6.ModflowGwtdsp(
+        gwt, xt3d_off=True, alh=alpha_l, ath1=alpha_th, atv=alpha_tv
+    )
+    pd = [
+        ("GWFHEAD", "../mf6gwf/flow.hds".format(), None),
+        ("GWFBUDGET", "../mf6gwf/flow.bud", None),
+    ]
+    flopy.mf6.ModflowGwtfmi(
+        gwt, flow_imbalance_correction=True, packagedata=pd
+    )
+    sourcerecarray = [
+        ("RCH-1", "AUX", "CONCENTRATION"),
+    ]
+    flopy.mf6.ModflowGwtssm(gwt, sources=sourcerecarray)
+    saverecord = {
+        0: [
+            ("CONCENTRATION", "STEPS", 10),
+            ("CONCENTRATION", "LAST"),
+            ("CONCENTRATION", "FREQUENCY", 10),
+        ],
+        1: [
+            ("CONCENTRATION", "STEPS", 27, 227),
+            ("CONCENTRATION", "LAST"),
+            ("CONCENTRATION", "FREQUENCY", 10),
+        ],
+    }
+    flopy.mf6.ModflowGwtoc(
+        gwt,
+        budget_filerecord="{}.cbc".format(name),
+        concentration_filerecord="{}.ucn".format(name),
+        concentrationprintrecord=[
+            ("COLUMNS", ncol, "WIDTH", 15, "DIGITS", 6, "GENERAL")
+        ],
+        saverecord=saverecord,
+        printrecord=[
+            ("CONCENTRATION", "LAST"),
+            (
+                "BUDGET",
+                "ALL",
+            ),
+        ],
+    )
+    obs_data = {
+        "{}.obs.csv".format(name): [
+            ("obs1", "CONCENTRATION", obs1),
+            ("obs2", "CONCENTRATION", obs2),
+        ],
+    }
+    flopy.mf6.ModflowUtlobs(
+        gwt, digits=10, print_input=True, continuous=obs_data
+    )
+    return sim
+
+
+def build_model(sim_name):
+    sims = None
+    sim_mf6gwf = build_mf6gwf(sim_name)
+    sim_mf6gwt = build_mf6gwt(sim_name)
+    sim_mf2005 = None  # build_mf2005(sim_name)
+    sim_mt3dms = None  # build_mt3dms(sim_name, sim_mf2005)
+    sims = (sim_mf6gwf, sim_mf6gwt, sim_mf2005, sim_mt3dms)
+    return sims
+
+
+def write_model(sims, silent=True):
+    sim_mf6gwf, sim_mf6gwt, sim_mf2005, sim_mt3dms = sims
+    sim_mf6gwf.write_simulation(silent=silent)
+    sim_mf6gwt.write_simulation(silent=silent)
+
+
+def run_keating_model(sims=example_name, silent=True):
+    sim = build_model(sims)
+    write_model(sim, silent=silent)
+    success = False
+    sim_mf6gwf, sim_mf6gwt, sim_mf2005, sim_mt3dms = sim
+    print("Running mf6gwf model...")
+    success, buff = sim_mf6gwf.run_simulation(silent=silent)
+    if not success:
+        print(buff)
+    print("Running mf6gwt model...")
+    success, buff = sim_mf6gwt.run_simulation(silent=silent)
+    if not success:
+        print(buff)
+    return success
+
+
 if __name__ == "__main__":
     run()
+    run_keating_model()

flopy/discretization/grid.py

@@ -169,6 +169,7 @@ class Grid(object):
         if angrot is None:
             angrot = 0.0
         self._angrot = angrot
+        self._polygons = None
         self._cache_dict = {}
         self._copy_cache = True
 
@@ -289,6 +290,10 @@ class Grid(object):
     def idomain(self):
         return copy.deepcopy(self._idomain)
 
+    @property
+    def ncpl(self):
+        raise NotImplementedError("must define ncpl in child class")
+
     @property
     def nnodes(self):
         raise NotImplementedError("must define nnodes in child class")
@@ -372,6 +377,38 @@ class Grid(object):
     #    raise NotImplementedError(
     #        'must define indices in child '
     #        'class to use this base class')
+    @property
+    def cross_section_vertices(self):
+        return self.xyzvertices[0], self.xyzvertices[1]
+
+    @property
+    def map_polygons(self):
+        """
+        Get a list of matplotlib Polygon patches for plotting
+
+        Returns
+        -------
+            list of Polygon objects
+        """
+        try:
+            from matplotlib.patches import Polygon
+        except ImportError:
+            raise ImportError("matplotlib required to use this method")
+        cache_index = "xyzgrid"
+        if (
+            cache_index not in self._cache_dict
+            or self._cache_dict[cache_index].out_of_date
+        ):
+            self.xyzvertices
+            self._polygons = None
+
+        if self._polygons is None:
+            self._polygons = [
+                Polygon(self.get_cell_vertices(nn), closed=True)
+                for nn in range(self.ncpl)
+            ]
+
+        return copy.copy(self._polygons)
 
     def get_plottable_layer_array(self, plotarray, layer):
         raise NotImplementedError(
@@ -430,11 +467,11 @@ class Grid(object):
         if not np.isscalar(x):
             x, y = x.copy(), y.copy()
 
-        x, y = geometry.rotate(
-            x, y, self._xoff, self._yoff, -self.angrot_radians
+        x, y = geometry.transform(
+            x, y, self._xoff, self._yoff, self.angrot_radians, inverse=True
         )
-        x -= self._xoff
-        y -= self._yoff
+        # x -= self._xoff
+        # y -= self._yoff
 
         return x, y
 

flopy/discretization/structuredgrid.py

@@ -228,6 +228,10 @@ class StructuredGrid(Grid):
     def ncol(self):
         return self.__ncol
 
+    @property
+    def ncpl(self):
+        return self.__nrow * self.__ncol
+
     @property
     def nnodes(self):
         return self.__nlay * self.__nrow * self.__ncol
@@ -800,15 +804,33 @@ class StructuredGrid(Grid):
             vrts = np.array(pts).transpose([2, 0, 1])
             return [v.tolist() for v in vrts]
 
-    def get_cell_vertices(self, i, j):
+    def get_cell_vertices(self, *args, **kwargs):
         """
         Method to get a set of cell vertices for a single cell
-            used in the Shapefile export utilities
+            used in the Shapefile export utilities and plotting code
+        :param node: (int) node number
         :param i: (int) cell row number
         :param j: (int) cell column number
         Returns
         ------- list of x,y cell vertices
         """
+        nn = None
+        if kwargs:
+            if "node" in kwargs:
+                nn = kwargs.pop("node")
+            else:
+                i = kwargs.pop("i")
+                j = kwargs.pop("j")
+
+        if len(args) > 0:
+            if len(args) == 1:
+                nn = args[0]
+            else:
+                i, j = args[0:2]
+
+        if nn is not None:
+            k, i, j = self.get_lrc(nn)[0]
+
         self._copy_cache = False
         cell_verts = [
             (self.xvertices[i, j], self.yvertices[i, j]),
@@ -819,6 +841,48 @@ class StructuredGrid(Grid):
         self._copy_cache = True
         return cell_verts
 
+    def get_lrc(self, nodes):
+        """
+        Get layer, row, column from a list of zero based
+        MODFLOW node numbers.
+
+        Returns
+        -------
+        v : list of tuples containing the layer (k), row (i),
+            and column (j) for each node in the input list
+        """
+        if not isinstance(nodes, list):
+            nodes = [nodes]
+        ncpl = self.ncpl
+        v = []
+        for node in nodes:
+            k = int(np.floor(node / ncpl))
+            ij = int((node) - (ncpl * k))
+            i = int(np.floor(ij / self.__ncol))
+            j = int(ij - (i * self.__ncol))
+
+            v.append((k, i, j))
+        return v
+
+    def get_node(self, lrc_list):
+        """
+        Get node number from a list of zero based MODFLOW
+        layer, row, column tuples.
+
+        Returns
+        -------
+        v : list of MODFLOW nodes for each layer (k), row (i),
+            and column (j) tuple in the input list
+        """
+        if not isinstance(lrc_list, list):
+            lrc_list = [lrc_list]
+        nrc = self.__nrow * self.__ncol
+        v = []
+        for [k, i, j] in lrc_list:
+            node = int(((k) * nrc) + ((i) * self.__ncol) + j)
+            v.append(node)
+        return v
+
     def plot(self, **kwargs):
         """
         Plot the grid lines.
@@ -1370,6 +1434,44 @@ class StructuredGrid(Grid):
 
         return afaces
 
+    @property
+    def cross_section_vertices(self):
+        """
+        Get a set of xvertices and yvertices ordered by node
+        for plotting cross sections
+
+        Returns
+        -------
+            xverts, yverts: (np.ndarray, np.ndarray)
+
+        """
+        xv = self.xyzvertices[0]
+        yv = self.xyzvertices[1]
+
+        xverts, yverts = [], []
+        for i in range(self.nrow):
+            for j in range(self.ncol):
+                xverts.append(
+                    [
+                        xv[i, j],
+                        xv[i + 1, j],
+                        xv[i + 1, j + 1],
+                        xv[i, j + 1],
+                        xv[i, j],
+                    ]
+                )
+                yverts.append(
+                    [
+                        yv[i, j],
+                        yv[i + 1, j],
+                        yv[i + 1, j + 1],
+                        yv[i, j + 1],
+                        yv[i, j],
+                    ]
+                )
+
+        return np.array(xverts), np.array(yverts)
+
     def get_number_plottable_layers(self, a):
         """
         Calculate and return the number of 2d plottable arrays that can be
@@ -1414,28 +1516,3 @@ class StructuredGrid(Grid):
         msg = "{} /= {}".format(plotarray.shape, required_shape)
         assert plotarray.shape == required_shape, msg
         return plotarray
-
-
-if __name__ == "__main__":
-    delc = np.ones((10,)) * 1
-    delr = np.ones((20,)) * 1
-
-    top = np.ones((10, 20)) * 2000
-    botm = np.ones((1, 10, 20)) * 1100
-
-    t = StructuredGrid(delc, delr, top, botm, xoff=0, yoff=0, angrot=45)
-
-    t.use_ref_coords = False
-    x = t.xvertices
-    y = t.yvertices
-    xc = t.xcellcenters
-    yc = t.ycellcenters
-    grid = t.grid_lines
-
-    t.use_ref_coords = True
-    sr_x = t.xvertices
-    sr_y = t.yvertices
-    sr_xc = t.xcellcenters
-    sr_yc = t.ycellcenters
-    sr_grid = t.grid_lines
-    print(sr_grid)

flopy/discretization/unstructuredgrid.py

@@ -1,3 +1,4 @@
+import copy
 import numpy as np
 from .grid import Grid, CachedData
 
@@ -301,6 +302,50 @@ class UnstructuredGrid(Grid):
         else:
             return self._cache_dict[cache_index].data_nocopy
 
+    @property
+    def map_polygons(self):
+        """
+        Property to get Matplotlib polygon objects for the modelgrid
+
+        Returns
+        -------
+            list or dict of matplotlib.collections.Polygon
+        """
+        try:
+            from matplotlib.patches import Polygon
+        except ImportError:
+            raise ImportError("matplotlib required to use this method")
+
+        cache_index = "xyzgrid"
+        if (
+            cache_index not in self._cache_dict
+            or self._cache_dict[cache_index].out_of_date
+        ):
+            self.xyzvertices
+            self._polygons = None
+
+        if self._polygons is None:
+            if self.grid_varies_by_layer:
+                self._polygons = {}
+                ilay = 0
+                lay_break = np.cumsum(self.ncpl)
+                for nn in range(self.nnodes):
+                    if nn in lay_break:
+                        ilay += 1
+
+                    if ilay not in self._polygons:
+                        self._polygons[ilay] = []
+
+                    p = Polygon(self.get_cell_vertices(nn), closed=True)
+                    self._polygons[ilay].append(p)
+            else:
+                self._polygons = [
+                    Polygon(self.get_cell_vertices(nn), closed=True)
+                    for nn in range(self.ncpl[0])
+                ]
+
+        return copy.copy(self._polygons)
+
     def intersect(self, x, y, local=False, forgive=False):
         x, y = super(UnstructuredGrid, self).intersect(x, y, local, forgive)
         raise Exception("Not implemented yet")

flopy/discretization/vertexgrid.py

@@ -1,3 +1,4 @@
+import copy
 import numpy as np
 
 try:
@@ -111,6 +112,8 @@ class VertexGrid(Grid):
             return len(self._cell1d)
         if self._botm is not None:
             return len(self._botm[0])
+        if self._cell2d is not None and self._nlay is None:
+            return len(self._cell2d)
         else:
             return self._ncpl
 
@@ -262,6 +265,15 @@ class VertexGrid(Grid):
         Returns
         ------- list of x,y cell vertices
         """
+        while cellid >= self.ncpl:
+            if cellid > self.nnodes:
+                err = "cellid {} out of index for size {}".format(
+                    cellid, self.nnodes
+                )
+                raise IndexError(err)
+
+            cellid -= self.ncpl
+
         self._copy_cache = False
         cell_verts = list(zip(self.xvertices[cellid], self.yvertices[cellid]))
         self._copy_cache = True
@@ -361,7 +373,7 @@ class VertexGrid(Grid):
             yvertices = yvertxform
 
         self._cache_dict[cache_index_cc] = CachedData(
-            [xcenters, ycenters, zcenters]
+            [np.array(xcenters), np.array(ycenters), np.array(zcenters)]
         )
         self._cache_dict[cache_index_vert] = CachedData(
             [xvertices, yvertices, zvertices]

flopy/mf6/utils/reference.py

@@ -206,7 +206,6 @@ class StructuredSpatialReference(object):
     @classmethod
     def from_gridspec(cls, gridspec_file, lenuni=0):
         f = open(gridspec_file, "r")
-        lines = f.readlines()
         raw = f.readline().strip().split()
         nrow = int(raw[0])
         ncol = int(raw[1])
@@ -220,10 +219,10 @@ class StructuredSpatialReference(object):
                 if "*" in r:
                     rraw = r.split("*")
                     for n in range(int(rraw[0])):
-                        delr.append(int(rraw[1]))
+                        delr.append(float(rraw[1]))
                         j += 1
                 else:
-                    delr.append(int(r))
+                    delr.append(float(r))
                     j += 1
         delc = []
         i = 0
@@ -233,10 +232,10 @@ class StructuredSpatialReference(object):
                 if "*" in r:
                     rraw = r.split("*")
                     for n in range(int(rraw[0])):
-                        delc.append(int(rraw[1]))
+                        delc.append(float(rraw[1]))
                         i += 1
                 else:
-                    delc.append(int(r))
+                    delc.append(float(r))
                     i += 1
         f.close()
         return cls(

flopy/modflow/mfdis.py

@@ -559,17 +559,7 @@ class ModflowDis(Package):
         v : list of tuples containing the layer (k), row (i),
             and column (j) for each node in the input list
         """
-        if not isinstance(nodes, list):
-            nodes = [nodes]
-        nrc = self.nrow * self.ncol
-        v = []
-        for node in nodes:
-            k = int(node / nrc)
-            ij = node - k * nrc
-            i = int(ij / self.ncol)
-            j = ij - i * self.ncol
-            v.append((k, i, j))
-        return v
+        return self.parent.modelgrid.get_lrc(nodes)
 
     def get_node(self, lrc_list):
         """
@@ -581,14 +571,7 @@ class ModflowDis(Package):
         v : list of MODFLOW nodes for each layer (k), row (i),
             and column (j) tuple in the input list
         """
-        if not isinstance(lrc_list, list):
-            lrc_list = [lrc_list]
-        nrc = self.nrow * self.ncol
-        v = []
-        for [k, i, j] in lrc_list:
-            node = int((k * nrc) + (i * self.ncol) + j)
-            v.append(node)
-        return v
+        return self.parent.modelgrid.get_node(lrc_list)
 
     def get_layer(self, i, j, elev):
         """Return the layer for an elevation at an i, j location.

flopy/modflow/mfupw.py

@@ -342,7 +342,7 @@ class ModflowUpw(Package):
         transient = not self.parent.get_package("DIS").steady.all()
         for k in range(nlay):
             f_upw.write(self.hk[k].get_file_entry())
-            if self.chani[k] < 1:
+            if self.chani[k] < 0:
                 f_upw.write(self.hani[k].get_file_entry())
             f_upw.write(self.vka[k].get_file_entry())
             if transient:
@@ -499,7 +499,7 @@ class ModflowUpw(Package):
             hk[k] = t
 
             # hani
-            if chani[k] < 1:
+            if chani[k] < 0:
                 if model.verbose:
                     print("   loading hani layer {0:3d}...".format(k + 1))
                 if "hani" not in par_types:

flopy/plot/__init__.py

@@ -27,5 +27,5 @@ from .plotutil import (
     PlotUtilities,
 )
 from .map import ModelMap, PlotMapView
-from .crosssection import ModelCrossSection
-from .plotbase import PlotCrossSection
+from .crosssection import ModelCrossSection, PlotCrossSection
+from .styles import styles

flopy/plot/map.py

@@ -7,7 +7,7 @@ try:
     import matplotlib.colors
     from matplotlib.collections import PatchCollection
     from matplotlib.patches import Polygon
-except ImportError:
+except (ImportError, ModuleNotFoundError):
     plt = None
 
 from . import plotutil
@@ -58,23 +58,19 @@ class PlotMapView(object):
         self.layer = layer
         self.mg = None
 
-        if model is not None:
-            self.mg = model.modelgrid
-        elif modelgrid is not None:
+        if modelgrid is not None:
             self.mg = modelgrid
+        elif model is not None:
+            self.mg = model.modelgrid
         else:
             err_msg = "A model grid instance must be provided to PlotMapView"
             raise AssertionError(err_msg)
 
-        if self.mg.grid_type not in ("structured", "vertex", "unstructured"):
-            err_msg = "Unrecognized modelgrid type {}"
-            raise TypeError(err_msg.format(self.mg.grid_type))
-
         if ax is None:
             try:
                 self.ax = plt.gca()
                 self.ax.set_aspect("equal")
-            except:
+            except (AttributeError, ValueError):
                 self.ax = plt.subplot(1, 1, 1, aspect="equal", axisbg="white")
         else:
             self.ax = ax
@@ -111,16 +107,12 @@ class PlotMapView(object):
 
         """
 
-        if self.mg.grid_type not in ("structured", "vertex", "unstructured"):
-            raise TypeError(
-                "Unrecognized grid type {}".format(self.mg.grid_type)
-            )
-
         if not isinstance(a, np.ndarray):
             a = np.array(a)
 
         # Use the model grid to pass back an array of the correct shape
         plotarray = self.mg.get_plottable_layer_array(a, self.layer)
+        plotarray = plotarray.ravel()
 
         # if masked_values are provided mask the plotting array
         if masked_values is not None:
@@ -130,50 +122,33 @@ class PlotMapView(object):
         # add NaN values to mask
         plotarray = np.ma.masked_where(np.isnan(plotarray), plotarray)
 
-        if "ax" in kwargs:
-            ax = kwargs.pop("ax")
-        else:
-            ax = self.ax
+        ax = kwargs.pop("ax", self.ax)
 
-        # Get vertices for the selected layer
-        xgrid = self.mg.get_xvertices_for_layer(self.layer)
-        ygrid = self.mg.get_yvertices_for_layer(self.layer)
+        # use cached patch collection for plotting
+        polygons = self.mg.map_polygons
+        if isinstance(polygons, dict):
+            polygons = polygons[self.layer]
 
-        if self.mg.grid_type == "structured":
-            quadmesh = ax.pcolormesh(xgrid, ygrid, plotarray)
-        else:
-            # use patch collection for vertex and unstructured
-            patches = [
-                Polygon(list(zip(xgrid[i], ygrid[i])), closed=True)
-                for i in range(xgrid.shape[0])
-            ]
-            quadmesh = PatchCollection(patches)
-            quadmesh.set_array(plotarray)
+        collection = PatchCollection(polygons)
+        collection.set_array(plotarray)
 
         # set max and min
-        if "vmin" in kwargs:
-            vmin = kwargs.pop("vmin")
-        else:
-            vmin = None
-
-        if "vmax" in kwargs:
-            vmax = kwargs.pop("vmax")
-        else:
-            vmax = None
+        vmin = kwargs.pop("vmin", None)
+        vmax = kwargs.pop("vmax", None)
 
         # limit the color range
-        quadmesh.set_clim(vmin=vmin, vmax=vmax)
+        collection.set_clim(vmin=vmin, vmax=vmax)
 
         # send rest of kwargs to quadmesh
-        quadmesh.set(**kwargs)
+        collection.set(**kwargs)
 
         # add collection to axis
-        ax.add_collection(quadmesh)
+        ax.add_collection(collection)
 
         # set limits
         ax.set_xlim(self.extent[0], self.extent[1])
         ax.set_ylim(self.extent[2], self.extent[3])
-        return quadmesh
+        return collection
 
     def contour_array(self, a, masked_values=None, **kwargs):
         """
@@ -210,15 +185,8 @@ class PlotMapView(object):
         # work around for tri-contour ignore vmin & vmax
         # necessary block for tri-contour NaN issue
         if "levels" not in kwargs:
-            if "vmin" not in kwargs:
-                vmin = np.nanmin(plotarray)
-            else:
-                vmin = kwargs.pop("vmin")
-            if "vmax" not in kwargs:
-                vmax = np.nanmax(plotarray)
-            else:
-                vmax = kwargs.pop("vmax")
-
+            vmin = kwargs.pop("vmin", np.nanmin(plotarray))
+            vmax = kwargs.pop("vmax", np.nanmax(plotarray))
             levels = np.linspace(vmin, vmax, 7)
             kwargs["levels"] = levels
 
@@ -241,10 +209,7 @@ class PlotMapView(object):
                     t = np.isclose(plotarray, mval)
                     ismasked += t
 
-        if "ax" in kwargs:
-            ax = kwargs.pop("ax")
-        else:
-            ax = self.ax
+        ax = kwargs.pop("ax", self.ax)
 
         if "colors" in kwargs.keys():
             if "cmap" in kwargs.keys():
@@ -311,9 +276,6 @@ class PlotMapView(object):
         quadmesh : matplotlib.collections.QuadMesh
 
         """
-        if plt is None:
-            err_msg = "matplotlib must be installed to use plot_inactive()"
-            raise ImportError(err_msg)
 
         if ibound is None:
             if self.mg.idomain is None:
@@ -359,9 +321,6 @@ class PlotMapView(object):
         quadmesh : matplotlib.collections.QuadMesh
 
         """
-        if plt is None:
-            err_msg = "matplotlib must be installed to use plot_ibound()"
-            raise ImportError(err_msg)
 
         if ibound is None:
             if self.model is not None:
@@ -399,32 +358,34 @@ class PlotMapView(object):
         lc : matplotlib.collections.LineCollection
 
         """
-        if plt is None:
-            err_msg = "matplotlib must be installed to use plot_grid()"
-            raise ImportError(err_msg)
+
+        from matplotlib.collections import PatchCollection
+
+        ax = kwargs.pop("ax", self.ax)
+        edgecolor = kwargs.pop("colors", "grey")
+        edgecolor = kwargs.pop("color", edgecolor)
+        edgecolor = kwargs.pop("ec", edgecolor)
+        edgecolor = kwargs.pop("edgecolor", edgecolor)
+        facecolor = kwargs.pop("facecolor", "none")
+        facecolor = kwargs.pop("fc", facecolor)
+
+        # use cached patch collection for plotting
+        polygons = self.mg.map_polygons
+        if isinstance(polygons, dict):
+            polygons = polygons[self.layer]
+
+        if len(polygons) > 0:
+            patches = PatchCollection(
+                polygons, edgecolor=edgecolor, facecolor=facecolor, **kwargs
+            )
         else:
-            from matplotlib.collections import LineCollection
+            patches = None
 
-        if "ax" in kwargs:
-            ax = kwargs.pop("ax")
-        else:
-            ax = self.ax
-
-        if "colors" not in kwargs:
-            kwargs["colors"] = "0.5"
-
-        grid_lines = self.mg.grid_lines
-        if isinstance(grid_lines, dict):
-            # grid_lines are passed back as a dictionary with keys equal to
-            # layers for an UnstructuredGrid
-            grid_lines = grid_lines[self.layer]
-        lc = LineCollection(grid_lines, **kwargs)
-
-        ax.add_collection(lc)
+        ax.add_collection(patches)
         ax.set_xlim(self.extent[0], self.extent[1])
         ax.set_ylim(self.extent[2], self.extent[3])
 
-        return lc
+        return patches
 
     def plot_bc(
         self,
@@ -522,15 +483,11 @@ class PlotMapView(object):
                 else:
                     idx = mflist["node"]
 
-        if plotAll and self.mg.grid_type == "unstructured":
-            raise Exception("plotAll cannot be used with unstructured grid.")
-        else:
-            nlay = self.mg.nlay
+        nlay = self.mg.nlay
 
-        # Plot the list locations
-        plotarray = np.zeros(self.mg.shape, dtype=int)
-        if plotAll:
-            pa = np.zeros(self.mg.shape[1:], dtype=int)
+        plotarray = np.zeros(self.mg.shape, dtype=np.int)
+        if plotAll and len(self.mg.shape) > 1:
+            pa = np.zeros(self.mg.shape[1:], dtype=np.int)
             pa[tuple(idx[1:])] = 1
             for k in range(nlay):
                 plotarray[k] = pa.copy()
@@ -574,18 +531,43 @@ class PlotMapView(object):
             Keyword arguments passed to plotutil.plot_shapefile()
 
         """
-        if "ax" in kwargs:
-            ax = kwargs.pop("ax")
-        else:
-            ax = self.ax
-        patch_collection = plotutil.plot_shapefile(shp, ax, **kwargs)
+        return self.plot_shapes(shp, **kwargs)
 
+    def plot_shapes(self, obj, **kwargs):
+        """
+        Plot shapes is a method that facilitates plotting a collection
+        of geospatial objects
+
+        Parameters
+        ----------
+        obj : collection object
+            obj can accept the following types
+
+            str : shapefile name
+            shapefile.Reader object
+            list of [shapefile.Shape, shapefile.Shape,]
+            shapefile.Shapes object
+            flopy.utils.geometry.Collection object
+            list of [flopy.utils.geometry, ...] objects
+            geojson.GeometryCollection object
+            geojson.FeatureCollection object
+            shapely.GeometryCollection object
+            list of [[vertices], ...]
+        kwargs : dictionary
+            keyword arguments passed to plotutil.plot_shapefile()
+
+        Returns
+        -------
+            matplotlib.Collection object
+        """
+        ax = kwargs.pop("ax", self.ax)
+        patch_collection = plotutil.plot_shapefile(obj, ax, **kwargs)
         return patch_collection
 
     def plot_cvfd(self, verts, iverts, **kwargs):
         """
-        Plot a cvfd grid.  The vertices must be in the same coordinates as
-        the rotated and offset grid.
+        Plot a cvfd grid.  The vertices must be in the same
+        coordinates as the rotated and offset grid.
 
         Parameters
         ----------
@@ -598,20 +580,23 @@ class PlotMapView(object):
             Keyword arguments passed to plotutil.plot_cvfd()
 
         """
-        if "ax" in kwargs:
-            ax = kwargs.pop("ax")
-        else:
-            ax = self.ax
-        patch_collection = plotutil.plot_cvfd(
-            verts, iverts, ax, self.layer, **kwargs
+        warnings.warn(
+            "plot_cvfd will be deprecated and will be removed in version "
+            "3.3.5. Use plot_grid or plot_array",
+            PendingDeprecationWarning,
         )
-        return patch_collection
+        a = kwargs.pop("a", None)
+        if a is None:
+            return self.plot_grid(**kwargs)
+        else:
+            return self.plot_array(a, **kwargs)
 
     def contour_array_cvfd(self, vertc, a, masked_values=None, **kwargs):
         """
-        Contour a cvfd array.  If the array is three-dimensional, then the method
-        will contour the layer tied to this class (self.layer). The vertices
-        must be in the same coordinates as the rotated and offset grid.
+        Contour a cvfd array.  If the array is three-dimensional,
+        then the method will contour the layer tied to this class (self.layer).
+        The vertices must be in the same coordinates as the rotated and
+        offset grid.
 
         Parameters
         ----------
@@ -629,69 +614,13 @@ class PlotMapView(object):
         contour_set : matplotlib.pyplot.contour
 
         """
-        try:
-            import matplotlib.tri as tri
-        except ImportError:
-            err_msg = "matplotlib must be updated to use contour_array()"
-            raise ImportError(err_msg)
+        warnings.warn(
+            "contour_cvfd will be deprecated and removed in version 3.3.5. "
+            " Use contour_array",
+            PendingDeprecationWarning,
+        )
 
-        if "ncpl" in kwargs:
-            nlay = self.layer + 1
-            ncpl = kwargs.pop("ncpl")
-            if isinstance(ncpl, int):
-                i = int(ncpl)
-                ncpl = np.ones((nlay,), dtype=int) * i
-            elif isinstance(ncpl, list) or isinstance(ncpl, tuple):
-                ncpl = np.array(ncpl)
-            i0 = 0
-            i1 = 0
-            for k in range(nlay):
-                i0 = i1
-                i1 = i0 + ncpl[k]
-            # retain vertc in selected layer
-            vertc = vertc[i0:i1, :]
-        else:
-            i0 = 0
-            i1 = vertc.shape[0]
-
-        plotarray = a[i0:i1]
-
-        ismasked = None
-        if masked_values is not None:
-            for mval in masked_values:
-                if ismasked is None:
-                    ismasked = np.isclose(plotarray, mval)
-                else:
-                    t = np.isclose(plotarray, mval)
-                    ismasked += t
-
-        # add NaN values to mask
-        if ismasked is None:
-            ismasked = np.isnan(plotarray)
-        else:
-            ismasked += np.isnan(plotarray)
-
-        if "ax" in kwargs:
-            ax = kwargs.pop("ax")
-        else:
-            ax = self.ax
-
-        if "colors" in kwargs.keys():
-            if "cmap" in kwargs.keys():
-                kwargs.pop("cmap")
-
-        triang = tri.Triangulation(vertc[:, 0], vertc[:, 1])
-
-        if ismasked is not None:
-            ismasked = ismasked.flatten()
-            mask = np.any(
-                np.where(ismasked[triang.triangles], True, False), axis=1
-            )
-            triang.set_mask(mask)
-
-        contour_set = ax.tricontour(triang, plotarray, **kwargs)
-
-        return contour_set
+        return self.contour_array(a, masked_values=masked_values, **kwargs)
 
     def plot_vector(
         self,
@@ -734,34 +663,37 @@ class PlotMapView(object):
             result of the quiver function
 
         """
-        if "pivot" in kwargs:
-            pivot = kwargs.pop("pivot")
-        else:
-            pivot = "middle"
+        pivot = kwargs.pop("pivot", "middle")
+        ax = kwargs.pop("ax", self.ax)
 
-        if "ax" in kwargs:
-            ax = kwargs.pop("ax")
-        else:
-            ax = self.ax
+        # get ibound array to mask inactive cells
+        ib = np.ones((self.mg.nnodes,), dtype=int)
+        if self.mg.idomain is not None:
+            ib = self.mg.idomain.ravel()
 
-        # get actual values to plot
-        if self.mg.grid_type == "structured":
-            x = self.mg.xcellcenters[::istep, ::jstep]
-            y = self.mg.ycellcenters[::istep, ::jstep]
-            u = vx[self.layer, ::istep, ::jstep]
-            v = vy[self.layer, ::istep, ::jstep]
-        else:
-            x = self.mg.xcellcenters[::istep]
-            y = self.mg.ycellcenters[::istep]
-            u = vx[self.layer, ::istep]
-            v = vy[self.layer, ::istep]
+        xcentergrid = self.mg.get_xcellcenters_for_layer(self.layer)
+        ycentergrid = self.mg.get_ycellcenters_for_layer(self.layer)
+        vx = self.mg.get_plottable_layer_array(vx, self.layer)
+        vy = self.mg.get_plottable_layer_array(vy, self.layer)
+        ib = self.mg.get_plottable_layer_array(ib, self.layer)
+
+        try:
+            x = xcentergrid[::istep, ::jstep]
+            y = ycentergrid[::istep, ::jstep]
+            u = vx[::istep, ::jstep]
+            v = vy[::istep, ::jstep]
+            ib = ib[::istep, ::jstep]
+        except IndexError:
+            x = xcentergrid[::jstep]
+            y = ycentergrid[::jstep]
+            u = vx[::jstep]
+            v = vy[::jstep]
+            ib = ib[::jstep]
 
         # if necessary, copy to avoid changing the passed values
         if masked_values is not None or normalize:
-            import copy
-
-            u = copy.copy(u)
-            v = copy.copy(v)
+            u = np.copy(u)
+            v = np.copy(v)
 
         # mask values
         if masked_values is not None:
@@ -777,13 +709,13 @@ class PlotMapView(object):
             u[idx] /= vmag[idx]
             v[idx] /= vmag[idx]
 
+        u[ib == 0] = np.nan
+        v[ib == 0] = np.nan
+
         # rotate and plot, offsets must be zero since
         # these are vectors not locations
         urot, vrot = geometry.rotate(u, v, 0.0, 0.0, self.mg.angrot_radians)
-
-        # plot with quiver
         quiver = ax.quiver(x, y, urot, vrot, pivot=pivot, **kwargs)
-
         return quiver
 
     def plot_specific_discharge(
@@ -806,6 +738,10 @@ class PlotMapView(object):
             row frequency to plot. (Default is 1.)
         jstep : int
             column frequency to plot. (Default is 1.)
+        normalize : bool
+            boolean flag used to determine if discharge vectors should
+            be normalized using the magnitude of the specific discharge in each
+            cell. (default is False)
         kwargs : matplotlib.pyplot keyword arguments for the
             plt.quiver method.
 
@@ -816,22 +752,12 @@ class PlotMapView(object):
 
         """
         warnings.warn(
-            "plot_specific_discharge() has been deprecated. Use "
-            "plot_vector() instead, which should follow after "
-            "postprocessing.get_specific_discharge()",
+            "plot_specific_discharge() has been deprecated and will be "
+            "removed in version 3.3.5. Use plot_vector() instead, which "
+            "should follow after postprocessing.get_specific_discharge()",
             DeprecationWarning,
         )
 
-        if "pivot" in kwargs:
-            pivot = kwargs.pop("pivot")
-        else:
-            pivot = "middle"
-
-        if "ax" in kwargs:
-            ax = kwargs.pop("ax")
-        else:
-            ax = self.ax
-
         if isinstance(spdis, list):
             print(
                 "Warning: Selecting the final stress period from Specific"
@@ -848,39 +774,7 @@ class PlotMapView(object):
         qx[idx] = spdis["qx"]
         qy[idx] = spdis["qy"]
 
-        qx = self.mg.get_plottable_layer_array(qx, self.layer)
-        qy = self.mg.get_plottable_layer_array(qy, self.layer)
-
-        # Get vertices for the selected layer
-        xcentergrid = self.mg.get_xcellcenters_for_layer(self.layer)
-        ycentergrid = self.mg.get_ycellcenters_for_layer(self.layer)
-
-        if self.mg.grid_type == "structured":
-            x = xcentergrid[::istep, ::jstep]
-            y = ycentergrid[::istep, ::jstep]
-            u = qx[::istep, ::jstep]
-            v = qy[::istep, ::jstep]
-        else:
-            x = xcentergrid[::istep]
-            y = ycentergrid[::istep]
-            u = qx[::istep]
-            v = qy[::istep]
-
-        # normalize
-        if normalize:
-            vmag = np.sqrt(u ** 2.0 + v ** 2.0)
-            idx = vmag > 0.0
-            u[idx] /= vmag[idx]
-            v[idx] /= vmag[idx]
-
-        u[u == 0] = np.nan
-        v[v == 0] = np.nan
-
-        # Rotate and plot, offsets must be zero since
-        # these are vectors not locations
-        urot, vrot = geometry.rotate(u, v, 0.0, 0.0, self.mg.angrot_radians)
-        quiver = ax.quiver(x, y, urot, vrot, pivot=pivot, **kwargs)
-        return quiver
+        return self.plot_vector(qx, qy, istep, jstep, normalize, **kwargs)
 
     def plot_discharge(
         self,
@@ -928,9 +822,9 @@ class PlotMapView(object):
 
         """
         warnings.warn(
-            "plot_discharge() has been deprecated. Use "
-            "plot_vector() instead, which should follow after "
-            "postprocessing.get_specific_discharge()",
+            "plot_discharge() has been deprecated and will be replaced "
+            "in version 3.3.5. Use plot_vector() instead, which should "
+            "follow after postprocessing.get_specific_discharge()",
             DeprecationWarning,
         )
 
@@ -948,10 +842,6 @@ class PlotMapView(object):
                 )
                 raise AssertionError(err)
 
-            ib = np.ones((self.mg.nlay, self.mg.nrow, self.mg.ncol))
-            if self.mg.idomain is not None:
-                ib = self.mg.idomain
-
             delr = self.mg.delr
             delc = self.mg.delc
             top = np.copy(self.mg.top)
@@ -998,26 +888,8 @@ class PlotMapView(object):
             qx, qy, qz = plotutil.PlotUtilities.centered_specific_discharge(
                 frf, fff, flf, delr, delc, sat_thk
             )
-            ib = ib.ravel()
-            qx = qx.ravel()
-            qy = qy.ravel()
-            del qz
 
-            temp = []
-            for ix, val in enumerate(ib):
-                if val != 0:
-                    temp.append((ix + 1, qx[ix], qy[ix]))
-
-            spdis = np.recarray(
-                (len(temp),),
-                dtype=[("node", int), ("qx", float), ("qy", float)],
-            )
-            for ix, tup in enumerate(temp):
-                spdis[ix] = tup
-
-            return self.plot_specific_discharge(
-                spdis, istep=istep, jstep=jstep, normalize=normalize, **kwargs
-            )
+            return self.plot_vector(qx, qy, istep, jstep, normalize, **kwargs)
 
     def plot_pathline(self, pl, travel_time=None, **kwargs):
         """
@@ -1049,11 +921,8 @@ class PlotMapView(object):
         lc : matplotlib.collections.LineCollection
 
         """
-        if plt is None:
-            err_msg = "matplotlib must be installed to use plot_pathline()"
-            raise ImportError(err_msg)
-        else:
-            from matplotlib.collections import LineCollection
+
+        from matplotlib.collections import LineCollection
 
         # make sure pathlines is a list
         if not isinstance(pl, list):
@@ -1071,32 +940,12 @@ class PlotMapView(object):
         else:
             kon = self.layer
 
-        if "marker" in kwargs:
-            marker = kwargs.pop("marker")
-        else:
-            marker = None
-
-        if "markersize" in kwargs:
-            markersize = kwargs.pop("markersize")
-        elif "ms" in kwargs:
-            markersize = kwargs.pop("ms")
-        else:
-            markersize = None
-
-        if "markercolor" in kwargs:
-            markercolor = kwargs.pop("markercolor")
-        else:
-            markercolor = None
-
-        if "markerevery" in kwargs:
-            markerevery = kwargs.pop("markerevery")
-        else:
-            markerevery = 1
-
-        if "ax" in kwargs:
-            ax = kwargs.pop("ax")
-        else:
-            ax = self.ax
+        marker = kwargs.pop("marker", None)
+        markersize = kwargs.pop("markersize", None)
+        markersize = kwargs.pop("ms", markersize)
+        markercolor = kwargs.pop("markercolor", None)
+        markerevery = kwargs.pop("markerevery", 1)
+        ax = kwargs.pop("ax", self.ax)
 
         if "colors" not in kwargs:
             kwargs["colors"] = "0.5"
@@ -1104,39 +953,7 @@ class PlotMapView(object):
         linecol = []
         markers = []
         for p in pl:
-            if travel_time is None:
-                tp = p.copy()
-            else:
-                if isinstance(travel_time, str):
-                    if "<=" in travel_time:
-                        time = float(travel_time.replace("<=", ""))
-                        idx = p["time"] <= time
-                    elif "<" in travel_time:
-                        time = float(travel_time.replace("<", ""))
-                        idx = p["time"] < time
-                    elif ">=" in travel_time:
-                        time = float(travel_time.replace(">=", ""))
-                        idx = p["time"] >= time
-                    elif "<" in travel_time:
-                        time = float(travel_time.replace(">", ""))
-                        idx = p["time"] > time
-                    else:
-                        try:
-                            time = float(travel_time)
-                            idx = p["time"] <= time
-                        except:
-                            errmsg = (
-                                "flopy.map.plot_pathline travel_time "
-                                + "variable cannot be parsed. "
-                                + "Acceptable logical variables are , "
-                                + "<=, <, >=, and >. "
-                                + "You passed {}".format(travel_time)
-                            )
-                            raise Exception(errmsg)
-                else:
-                    time = float(travel_time)
-                    idx = p["time"] <= time
-                tp = p[idx]
+            tp = plotutil.filter_modpath_by_travel_time(p, travel_time)
 
             # transform data!
             x0r, y0r = geometry.transform(
@@ -1162,7 +979,7 @@ class PlotMapView(object):
                     linecol.append(arr)
                     if marker is not None:
                         for xy in arr[::markerevery]:
-                            if not xy.mask:
+                            if not np.all(xy.mask):
                                 markers.append(xy)
         # create line collection
         lc = None
@@ -1210,104 +1027,10 @@ class PlotMapView(object):
         -------
             lo : list of Line2D objects
         """
-        if plt is None:
-            err_msg = "matplotlib must be installed to use plot_timeseries()"
-            raise ImportError(err_msg)
+        if "color" in kwargs:
+            kwargs["markercolor"] = kwargs["color"]
 
-        # make sure timeseries is a list
-        if not isinstance(ts, list):
-            ts = [ts]
-
-        if "layer" in kwargs:
-            kon = kwargs.pop("layer")
-
-            if isinstance(kon, bytes):
-                kon = kon.decode()
-
-            if isinstance(kon, str):
-                if kon.lower() == "all":
-                    kon = -1
-                else:
-                    kon = self.layer
-        else:
-            kon = self.layer
-
-        if "ax" in kwargs:
-            ax = kwargs.pop("ax")
-
-        else:
-            ax = self.ax
-
-        if "color" not in kwargs:
-            kwargs["color"] = "red"
-
-        linecol = []
-        for t in ts:
-            if travel_time is None:
-                tp = t.copy()
-
-            else:
-                if isinstance(travel_time, str):
-                    if "<=" in travel_time:
-                        time = float(travel_time.replace("<=", ""))
-                        idx = t["time"] <= time
-                    elif "<" in travel_time:
-                        time = float(travel_time.replace("<", ""))
-                        idx = t["time"] < time
-                    elif ">=" in travel_time:
-                        time = float(travel_time.replace(">=", ""))
-                        idx = t["time"] >= time
-                    elif "<" in travel_time:
-                        time = float(travel_time.replace(">", ""))
-                        idx = t["time"] > time
-                    else:
-                        try:
-                            time = float(travel_time)
-                            idx = t["time"] <= time
-                        except:
-                            errmsg = (
-                                "flopy.map.plot_pathline travel_time "
-                                + "variable cannot be parsed. "
-                                + "Acceptable logical variables are , "
-                                + "<=, <, >=, and >. "
-                                + "You passed {}".format(travel_time)
-                            )
-                            raise Exception(errmsg)
-                else:
-                    time = float(travel_time)
-                    idx = t["time"] <= time
-                tp = ts[idx]
-
-            x0r, y0r = geometry.transform(
-                tp["x"],
-                tp["y"],
-                self.mg.xoffset,
-                self.mg.yoffset,
-                self.mg.angrot_radians,
-            )
-
-            # build polyline array
-            arr = np.vstack((x0r, y0r)).T
-            # select based on layer
-            if kon >= 0:
-                kk = t["k"].copy().reshape(t.shape[0], 1)
-                kk = np.repeat(kk, 2, axis=1)
-                arr = np.ma.masked_where((kk != kon), arr)
-
-            else:
-                arr = np.ma.asarray(arr)
-
-            # append line to linecol if there is some unmasked segment
-            if not arr.mask.all():
-                linecol.append(arr)
-
-        # plot timeseries data
-        lo = []
-        for lc in linecol:
-            if not lc.mask.all():
-                lo += ax.plot(lc[:, 0], lc[:, 1], **kwargs)
-
-        return lo
+        return self.plot_pathline(ts, travel_time=travel_time, **kwargs)
 
     def plot_endpoint(
         self,
@@ -1353,89 +1076,12 @@ class PlotMapView(object):
         sp : matplotlib.pyplot.scatter
 
         """
-        if plt is None:
-            err_msg = "matplotlib must be installed to use plot_endpoint()"
-            raise ImportError(err_msg)
 
-        ep = ep.copy()
-        direction = direction.lower()
-        if direction == "starting":
-            xp, yp = "x0", "y0"
-
-        elif direction == "ending":
-            xp, yp = "x", "y"
-
-        else:
-            errmsg = (
-                'flopy.map.plot_endpoint direction must be "ending" '
-                + 'or "starting".'
-            )
-            raise Exception(errmsg)
-
-        if selection_direction is not None:
-            if (
-                selection_direction.lower() != "starting"
-                and selection_direction.lower() != "ending"
-            ):
-                errmsg = (
-                    "flopy.map.plot_endpoint selection_direction "
-                    + 'must be "ending" or "starting".'
-                )
-                raise Exception(errmsg)
-        else:
-            if direction.lower() == "starting":
-                selection_direction = "ending"
-            elif direction.lower() == "ending":
-                selection_direction = "starting"
-
-        # selection of endpoints
-        if selection is not None:
-            if isinstance(selection, int):
-                selection = tuple((selection,))
-            try:
-                if len(selection) == 1:
-                    node = selection[0]
-                    if selection_direction.lower() == "starting":
-                        nsel = "node0"
-                    else:
-                        nsel = "node"
-                    # make selection
-                    idx = ep[nsel] == node
-                    tep = ep[idx]
-                elif len(selection) == 3:
-                    k, i, j = selection[0], selection[1], selection[2]
-                    if selection_direction.lower() == "starting":
-                        ksel, isel, jsel = "k0", "i0", "j0"
-                    else:
-                        ksel, isel, jsel = "k", "i", "j"
-                    # make selection
-                    idx = (ep[ksel] == k) & (ep[isel] == i) & (ep[jsel] == j)
-                    tep = ep[idx]
-                else:
-                    errmsg = (
-                        "flopy.map.plot_endpoint selection must be "
-                        + "a zero-based layer, row, column tuple "
-                        + "(l, r, c) or node number (MODPATH 7) of "
-                        + "the location to evaluate (i.e., well location)."
-                    )
-                    raise Exception(errmsg)
-            except:
-                errmsg = (
-                    "flopy.map.plot_endpoint selection must be a "
-                    + "zero-based layer, row, column tuple (l, r, c) "
-                    + "or node number (MODPATH 7) of the location "
-                    + "to evaluate (i.e., well location)."
-                )
-                raise Exception(errmsg)
-        # all endpoints
-        else:
-            tep = ep.copy()
-
-        if "ax" in kwargs:
-            ax = kwargs.pop("ax")
-        else:
-            ax = self.ax
+        ax = kwargs.pop("ax", self.ax)
 
+        tep, _, xp, yp = plotutil.parse_modpath_selection_options(
+            ep, direction, selection, selection_direction
+        )
         # scatter kwargs that users may redefine
         if "c" not in kwargs:
             c = tep["time"] - tep["time0"]
@@ -1443,24 +1089,13 @@ class PlotMapView(object):
             c = np.empty((tep.shape[0]), dtype="S30")
             c.fill(kwargs.pop("c"))
 
-        s = 50
-        if "s" in kwargs:
-            s = float(kwargs.pop("s")) ** 2.0
-        elif "size" in kwargs:
-            s = float(kwargs.pop("size")) ** 2.0
+        s = kwargs.pop("s", np.sqrt(50))
+        s = float(kwargs.pop("size", s)) ** 2.0
 
         # colorbar kwargs
-        createcb = False
-        if "colorbar" in kwargs:
-            createcb = kwargs.pop("colorbar")
-
-        colorbar_label = "Endpoint Time"
-        if "colorbar_label" in kwargs:
-            colorbar_label = kwargs.pop("colorbar_label")
-
-        shrink = 1.0
-        if "shrink" in kwargs:
-            shrink = float(kwargs.pop("shrink"))
+        createcb = kwargs.pop("colorbar", False)
+        colorbar_label = kwargs.pop("colorbar_label", "Endpoint Time")
+        shrink = float(kwargs.pop("shrink", 1.0))
 
         # transform data!
         x0r, y0r = geometry.transform(
@@ -1573,7 +1208,7 @@ class DeprecatedMapView(PlotMapView):
 
 class ModelMap(object):
     """
-    Pending Depreciation: ModelMap acts as a PlotMapView factory
+    DEPRECATED. ModelMap acts as a PlotMapView factory
     object. Please migrate to PlotMapView for plotting
     functionality and future code compatibility
 
@@ -1634,13 +1269,12 @@ class ModelMap(object):
 
         from ..utils.reference import SpatialReferenceUnstructured
 
-        # from ..plot.plotbase import DeprecatedMapView
-
         err_msg = (
-            "ModelMap will be replaced by "
-            "PlotMapView(); Calling PlotMapView()"
+            "ModelMap is deprecated and has been replaced by "
+            "PlotMapView(). ModelMap will be removed in version 3.3.5; "
+            "Calling PlotMapView()"
         )
-        warnings.warn(err_msg, PendingDeprecationWarning)
+        warnings.warn(err_msg, DeprecationWarning)
 
         modelgrid = None
         if model is not None:

flopy/plot/mplstyle/usgsmap.mplstyle

@@ -0,0 +1,28 @@
+font.family: sans-serif
+font.sans-serif: Arial
+font.size: 7
+axes.labelsize: 9
+axes.titlesize: 9
+axes.linewidth: 0.5
+xtick.labelsize: 7
+xtick.top: True
+xtick.bottom: True
+xtick.major.size: 7.2
+xtick.minor.size: 3.6
+xtick.major.width: 0.5
+xtick.minor.width: 0.5
+xtick.direction: in
+ytick.labelsize: 7
+ytick.left: True
+ytick.right: True
+ytick.major.size: 7.2
+ytick.minor.size: 3.6
+ytick.major.width: 0.5
+ytick.minor.width: 0.5
+ytick.direction: in
+pdf.fonttype: 42
+savefig.dpi: 300
+savefig.transparent: True
+legend.fontsize: 9
+legend.frameon: False
+legend.markerscale: 1.0

flopy/plot/mplstyle/usgsmap_linux.mplstyle

@@ -0,0 +1,28 @@
+font.family: sans-serif
+font.sans-serif: Liberation Sans
+font.size: 7
+axes.labelsize: 9
+axes.titlesize: 9
+axes.linewidth: 0.5
+xtick.labelsize: 7
+xtick.top: True
+xtick.bottom: True
+xtick.major.size: 7.2
+xtick.minor.size: 3.6
+xtick.major.width: 0.5
+xtick.minor.width: 0.5
+xtick.direction: in
+ytick.labelsize: 7
+ytick.left: True
+ytick.right: True
+ytick.major.size: 7.2
+ytick.minor.size: 3.6
+ytick.major.width: 0.5
+ytick.minor.width: 0.5
+ytick.direction: in
+pdf.fonttype: 42
+savefig.dpi: 300
+savefig.transparent: True
+legend.fontsize: 9
+legend.frameon: False
+legend.markerscale: 1.0

flopy/plot/mplstyle/usgsplot.mplstyle

@@ -0,0 +1,28 @@
+font.family: sans-serif
+font.sans-serif: Arial
+font.size: 7
+axes.labelsize: 9
+axes.titlesize: 9
+axes.linewidth: 0.5
+xtick.labelsize: 8
+xtick.top: True
+xtick.bottom: True
+xtick.major.size: 7.2
+xtick.minor.size: 3.6
+xtick.major.width: 0.5
+xtick.minor.width: 0.5
+xtick.direction: in
+ytick.labelsize: 8
+ytick.left: True
+ytick.right: True
+ytick.major.size: 7.2
+ytick.minor.size: 3.6
+ytick.major.width: 0.5
+ytick.minor.width: 0.5
+ytick.direction: in
+pdf.fonttype: 42
+savefig.dpi: 300
+savefig.transparent: True
+legend.fontsize: 9
+legend.frameon: False
+legend.markerscale: 1.0

flopy/plot/mplstyle/usgsplot_linux.mplstyle

@@ -0,0 +1,28 @@
+font.family: sans-serif
+font.sans-serif: Liberation Sans
+font.size: 7
+axes.labelsize: 9
+axes.titlesize: 9
+axes.linewidth: 0.5
+xtick.labelsize: 8
+xtick.top: True
+xtick.bottom: True
+xtick.major.size: 7.2
+xtick.minor.size: 3.6
+xtick.major.width: 0.5
+xtick.minor.width: 0.5
+xtick.direction: in
+ytick.labelsize: 8
+ytick.left: True
+ytick.right: True
+ytick.major.size: 7.2
+ytick.minor.size: 3.6
+ytick.major.width: 0.5
+ytick.minor.width: 0.5
+ytick.direction: in
+pdf.fonttype: 42
+savefig.dpi: 300
+savefig.transparent: True
+legend.fontsize: 9
+legend.frameon: False
+legend.markerscale: 1.0

flopy/plot/plotutil.py

@@ -9,6 +9,7 @@ import os
 import sys
 import math
 import numpy as np
+import warnings
 from ..utils import Util3d
 from ..datbase import DataType, DataInterface
 
@@ -22,6 +23,8 @@ try:
 except ImportError:
     plt = None
 
+warnings.simplefilter("ignore", RuntimeWarning)
+
 bc_color_dict = {
     "default": "black",
     "WEL": "red",
@@ -211,6 +214,7 @@ class PlotUtilities(object):
                     mflay=defaults["mflay"],
                     key=defaults["key"],
                     model_name=defaults["model_name"],
+                    model_grid=model.modelgrid,
                 )
                 # unroll nested lists of axes into a single list of axes
                 if isinstance(caxs, list):
@@ -236,9 +240,11 @@ class PlotUtilities(object):
                             mflay=defaults["mflay"],
                             key=defaults["key"],
                             model_name=defaults["model_name"],
+                            modelgrid=model.modelgrid,
                         )
 
-                        # unroll nested lists of axes into a single list of axes
+                        # unroll nested lists of axes into a single list
+                        # of axes
                         if isinstance(caxs, list):
                             for c in caxs:
                                 axes.append(c)
@@ -294,6 +300,7 @@ class PlotUtilities(object):
             "key": None,
             "initial_fig": 0,
             "model_name": "",
+            "modelgrid": None,
         }
 
         for key in defaults:
@@ -344,6 +351,7 @@ class PlotUtilities(object):
                                 fignum=fignum,
                                 model_name=model_name,
                                 colorbar=True,
+                                modelgrid=defaults["modelgrid"],
                             )
                         )
 
@@ -388,7 +396,8 @@ class PlotUtilities(object):
                         )
                     )
                     defaults["initial_fig"] = fignum[-1] + 1
-                    # need to keep this as value.plot() because of mf6 datatype issues
+                    # need to keep this as value.plot() because of
+                    # mf6 datatype issues
                     ax = value.plot(
                         defaults["key"],
                         names,
@@ -398,6 +407,7 @@ class PlotUtilities(object):
                         mflay=defaults["mflay"],
                         fignum=fignum,
                         colorbar=colorbar,
+                        modelgrid=defaults["modelgrid"],
                         **kwargs
                     )
 
@@ -433,6 +443,7 @@ class PlotUtilities(object):
                                 fignum=fignum,
                                 model_name=model_name,
                                 colorbar=True,
+                                modelgrid=defaults["modelgrid"],
                             )
                         )
 
@@ -463,6 +474,7 @@ class PlotUtilities(object):
                                     fignum=fignum,
                                     model_name=model_name,
                                     colorbar=True,
+                                    modelgrid=defaults["modelgrid"],
                                 )
                             )
 
@@ -492,6 +504,7 @@ class PlotUtilities(object):
                                 kper=defaults["kper"],
                                 fignum=fignum,
                                 colorbar=True,
+                                modelgrid=defaults["modelgrid"],
                             )
                         )
 
@@ -592,6 +605,10 @@ class PlotUtilities(object):
         if "model_name" in kwargs:
             model_name = kwargs.pop("model_name") + " "
 
+        modelgrid = None
+        if "modelgrid" in kwargs:
+            modelgrid = kwargs.pop("modelgrid")
+
         filenames = None
         if filename_base is not None:
             if mflay is not None:
@@ -638,6 +655,7 @@ class PlotUtilities(object):
                 names=names,
                 filenames=filenames,
                 mflay=mflay,
+                modelgrid=modelgrid,
                 **kwargs
             )
         else:
@@ -659,6 +677,7 @@ class PlotUtilities(object):
                 names=names,
                 filenames=filenames,
                 mflay=mflay,
+                modelgrid=modelgrid,
                 **kwargs
             )
         return axes
@@ -728,6 +747,10 @@ class PlotUtilities(object):
         if "model_name" in kwargs:
             model_name = kwargs.pop("model_name") + " "
 
+        modelgrid = None
+        if "modelgrid" in kwargs:
+            modelgrid = kwargs.pop("modelgrid")
+
         if title is None:
             title = "{}{}".format(model_name, util2d.name)
 
@@ -746,6 +769,7 @@ class PlotUtilities(object):
             names=title,
             filenames=filename,
             fignum=fignum,
+            modelgrid=modelgrid,
             **kwargs
         )
         return axes
@@ -815,6 +839,10 @@ class PlotUtilities(object):
         if "model_name" in kwargs:
             model_name = kwargs.pop("model_name")
 
+        modelgrid = None
+        if "modelgrid" in kwargs:
+            modelgrid = kwargs.pop("modelgrid")
+
         if file_extension is not None:
             fext = file_extension
         else:
@@ -859,6 +887,7 @@ class PlotUtilities(object):
             filenames=filenames,
             mflay=mflay,
             fignum=fignum,
+            modelgrid=modelgrid,
             **kwargs
         )
         return axes
@@ -932,6 +961,10 @@ class PlotUtilities(object):
         else:
             fext = "png"
 
+        modelgrid = None
+        if "modelgrid" in kwargs:
+            modelgrid = kwargs.pop("modelgrid")
+
         if isinstance(kper, int):
             k0 = kper
             k1 = kper + 1
@@ -976,6 +1009,7 @@ class PlotUtilities(object):
                     names=title,
                     filenames=filename,
                     fignum=fignum[idx],
+                    modelgrid=modelgrid,
                     **kwargs
                 )
             )
@@ -1012,6 +1046,10 @@ class PlotUtilities(object):
         if "mflay" in kwargs:
             kwargs.pop("mflay")
 
+        modelgrid = None
+        if "modelgrid" in kwargs:
+            modelgrid = kwargs.pop("modelgrid")
+
         title = "{}".format(scalar.name.replace("_", "").upper())
 
         if filename_base is not None:
@@ -1024,6 +1062,7 @@ class PlotUtilities(object):
             scalar.model,
             names=title,
             filenames=filename,
+            modelgrid=modelgrid,
             **kwargs
         )
         return axes
@@ -1569,8 +1608,9 @@ class PlotUtilities(object):
         """
         DEPRECATED. Use postprocessing.get_specific_discharge() instead.
 
-        Using the MODFLOW discharge, calculate the cell centered specific discharge
-        by dividing by the flow width and then averaging to the cell center.
+        Using the MODFLOW discharge, calculate the cell centered specific
+        discharge by dividing by the flow width and then averaging
+        to the cell center.
 
         Parameters
         ----------
@@ -1580,8 +1620,9 @@ class PlotUtilities(object):
             MODFLOW 'flow front face'.  The sign on this array will be flipped
             by this function so that the y axis is positive to north.
         Qz : numpy.ndarray
-            MODFLOW 'flow lower face'.  The sign on this array will be flipped by
-            this function so that the z axis is positive in the upward direction.
+            MODFLOW 'flow lower face'.  The sign on this array will be
+            flipped by this function so that the z axis is positive
+            in the upward direction.
         delr : numpy.ndarray
             MODFLOW delr array
         delc : numpy.ndarray
@@ -1598,7 +1639,8 @@ class PlotUtilities(object):
         import warnings
 
         warnings.warn(
-            "centered_specific_discharge() has been deprecated. Use "
+            "centered_specific_discharge() has been deprecated and will be "
+            "removed in version 3.3.5. Use "
             "postprocessing.get_specific_discharge() instead.",
             DeprecationWarning,
         )
@@ -1733,13 +1775,13 @@ class UnstructuredPlotUtilities(object):
                     # only cycle through the cells that intersect
                     # the infinite line
                     cvert_ix = []
-                    for ix in range(len(cpv)):
-                        if cpv[ix - 1] < 0 and cpv[ix] > 0:
-                            cvert_ix.append(ix - 1)
-                        elif cpv[ix - 1] > 0 and cpv[ix] < 0:
-                            cvert_ix.append(ix - 1)
-                        elif cpv[ix - 1] == 0 and cpv[ix] == 0:
-                            cvert_ix += [ix - 1, ix]
+                    for vx in range(len(cpv)):
+                        if cpv[vx - 1] < 0 and cpv[vx] > 0:
+                            cvert_ix.append(vx - 1)
+                        elif cpv[vx - 1] > 0 and cpv[vx] < 0:
+                            cvert_ix.append(vx - 1)
+                        elif cpv[vx - 1] == 0 and cpv[vx] == 0:
+                            cvert_ix += [vx - 1, vx]
                         else:
                             pass
 
@@ -1760,13 +1802,13 @@ class UnstructuredPlotUtilities(object):
             x = x1 + ua * (x2 - x1)
             y = y1 + ua * (y2 - y1)
 
-            for ix, cell in enumerate(cells):
+            for iix, cell in enumerate(cells):
                 xc = x[cell]
                 yc = y[cell]
                 verts = [
                     (xt, yt)
                     for xt, yt in zip(
-                        xc[cell_vertex_ix[ix]], yc[cell_vertex_ix[ix]]
+                        xc[cell_vertex_ix[iix]], yc[cell_vertex_ix[iix]]
                     )
                 ]
 
@@ -1866,9 +1908,10 @@ class UnstructuredPlotUtilities(object):
         return xverts, yverts
 
     @staticmethod
-    def arctan2(verts):
+    def arctan2(verts, reverse=False):
         """
-        Reads 2 dimensional set of verts and orders them using the arctan 2 method
+        Reads 2 dimensional set of verts and orders them using the
+        arctan 2 method
 
         Parameters
         ----------
@@ -1889,6 +1932,8 @@ class UnstructuredPlotUtilities(object):
         angleidx = angles.argsort()
 
         verts = verts[angleidx]
+        if reverse:
+            return verts[::-1]
         return verts
 
 
@@ -2096,7 +2141,10 @@ def shapefile_to_patch_collection(shp, radius=500.0, idx=None):
 
     """
     if shapefile is None:
-        s = "Could not import shapefile.  Must install pyshp in order to plot shapefiles."
+        s = (
+            "Could not import shapefile.  Must install pyshp "
+            "in order to plot shapefiles."
+        )
         raise PlotException(s)
     if plt is None:
         err_msg = (
@@ -2105,13 +2153,15 @@ def shapefile_to_patch_collection(shp, radius=500.0, idx=None):
         )
         raise ImportError(err_msg)
     else:
-        from matplotlib.patches import Polygon, Circle, Path, PathPatch
+        from matplotlib.patches import Polygon, Circle, PathPatch
+        import matplotlib.path as MPath
         from matplotlib.collections import PatchCollection
-    if isinstance(shp, str):
-        sf = shapefile.Reader(shp)
-    else:
-        sf = shp
-    shapes = sf.shapes()
+        from ..utils.geospatial_utils import GeoSpatialCollection
+        from ..utils.geometry import point_in_polygon
+
+    geofeats = GeoSpatialCollection(shp)
+    shapes = geofeats.shape
+
     nshp = len(shapes)
     ptchs = []
     if idx is None:
@@ -2127,16 +2177,66 @@ def shapefile_to_patch_collection(shp, radius=500.0, idx=None):
             vertices = []
             for p in shapes[n].points:
                 vertices.append([p[0], p[1]])
+            vertices += vertices[::-1]
             vertices = np.array(vertices)
-            path = Path(vertices)
-            ptchs.append(PathPatch(path, fill=False))
+            ptchs.append(Polygon(vertices))
         elif st in [5, 25, 31]:
             # polygons
             pts = np.array(shapes[n].points)
             prt = shapes[n].parts
             par = list(prt) + [pts.shape[0]]
+            polys = []
             for pij in range(len(prt)):
-                ptchs.append(Polygon(pts[par[pij] : par[pij + 1]]))
+                poly = np.array(pts[par[pij] : par[pij + 1]])
+                if not polys:
+                    polys.append(poly)
+                else:
+                    temp = []
+                    for ix, p in enumerate(polys):
+                        # check multipolygons for holes!
+                        mask = point_in_polygon(
+                            poly.T[0].reshape(1, -1),
+                            poly.T[1].reshape(1, -1),
+                            p,
+                        )
+
+                        if np.all(mask):
+                            temp.append((poly, ix))
+                        else:
+                            temp.append((poly, -1))
+
+                    for p, flag in temp:
+                        if flag < 0:
+                            polys.append(p)
+                        else:
+                            # hole in polygon
+                            if isinstance(polys[flag], list):
+                                polys[flag].append(p)
+                            else:
+                                polys[flag] = [polys[flag], p]
+
+            for poly in polys:
+                if isinstance(poly, list):
+                    codes = []
+                    for path in poly:
+                        c = (
+                            np.ones(len(path), dtype=MPath.Path.code_type)
+                            * MPath.Path.LINETO
+                        )
+                        c[0] = MPath.Path.MOVETO
+                        if len(codes) == 0:
+                            codes = c
+                            verts = path
+                        else:
+                            codes = np.concatenate((codes, c))
+                            verts = np.concatenate((verts, path))
+
+                    mplpath = MPath.Path(verts, codes)
+                    ptchs.append(PathPatch(mplpath))
+
+                else:
+                    ptchs.append(Polygon(poly))
+
     pc = PatchCollection(ptchs)
     return pc
 
@@ -2191,18 +2291,14 @@ def plot_shapefile(
     """
 
     if shapefile is None:
-        s = "Could not import shapefile.  Must install pyshp in order to plot shapefiles."
+        s = (
+            "Could not import shapefile.  Must install pyshp in "
+            "order to plot shapefiles."
+        )
         raise PlotException(s)
 
-    if "vmin" in kwargs:
-        vmin = kwargs.pop("vmin")
-    else:
-        vmin = None
-
-    if "vmax" in kwargs:
-        vmax = kwargs.pop("vmax")
-    else:
-        vmax = None
+    vmin = kwargs.pop("vmin", None)
+    vmax = kwargs.pop("vmax", None)
 
     if ax is None:
         ax = plt.gca()
@@ -2248,6 +2344,12 @@ def cvfd_to_patch_collection(verts, iverts):
         should be of len(ncells) with a list of vertex numbers for each cell
 
     """
+    warnings.warn(
+        "cvfd_to_patch_collection is deprecated and will be removed in "
+        "version 3.3.5. Use PlotMapView for plotting",
+        DeprecationWarning,
+    )
+
     if plt is None:
         err_msg = (
             "matplotlib must be installed to "
@@ -2321,6 +2423,11 @@ def plot_cvfd(
     --------
 
     """
+    warnings.warn(
+        "plot_cvfd is deprecated and will be removed in version 3.3.5. "
+        "Use PlotMapView for plotting",
+        DeprecationWarning,
+    )
     if plt is None:
         err_msg = "matplotlib must be installed to use plot_cvfd()"
         raise ImportError(err_msg)
@@ -2411,252 +2518,6 @@ def plot_cvfd(
     return pc
 
 
-def findrowcolumn(pt, xedge, yedge):
-    """
-    Find the MODFLOW cell containing the x- and y- point provided.
-
-    Parameters
-    ----------
-    pt : list or tuple
-        A list or tuple containing a x- and y- coordinate
-    xedge : numpy.ndarray
-        x-coordinate of the edge of each MODFLOW column. xedge is dimensioned
-        to NCOL + 1. If xedge is not a numpy.ndarray it is converted to a
-        numpy.ndarray.
-    yedge : numpy.ndarray
-        y-coordinate of the edge of each MODFLOW row. yedge is dimensioned
-        to NROW + 1. If yedge is not a numpy.ndarray it is converted to a
-        numpy.ndarray.
-
-    Returns
-    -------
-    irow, jcol : int
-        Row and column location containing x- and y- point passed to function.
-
-    Examples
-    --------
-    >>> import flopy
-    >>> irow, jcol = flopy.plotutil.findrowcolumn(pt, xedge, yedge)
-
-    """
-
-    # make sure xedge and yedge are numpy arrays
-    if not isinstance(xedge, np.ndarray):
-        xedge = np.array(xedge)
-    if not isinstance(yedge, np.ndarray):
-        yedge = np.array(yedge)
-
-    # find column
-    jcol = -100
-    for jdx, xmf in enumerate(xedge):
-        if xmf > pt[0]:
-            jcol = jdx - 1
-            break
-
-    # find row
-    irow = -100
-    for jdx, ymf in enumerate(yedge):
-        if ymf < pt[1]:
-            irow = jdx - 1
-            break
-    return irow, jcol
-
-
-def line_intersect_grid(ptsin, xedge, yedge, returnvertices=False):
-    """
-    Intersect a list of polyline vertices with a rectilinear MODFLOW
-    grid. Vertices at the intersection of the polyline with the grid
-    cell edges is returned. Optionally the original polyline vertices
-    are returned.
-
-    Parameters
-    ----------
-    ptsin : list
-        A list of x, y points defining the vertices of a polyline that will be
-        intersected with the rectilinear MODFLOW grid
-    xedge : numpy.ndarray
-        x-coordinate of the edge of each MODFLOW column. xedge is dimensioned
-        to NCOL + 1. If xedge is not a numpy.ndarray it is converted to a
-        numpy.ndarray.
-    yedge : numpy.ndarray
-        y-coordinate of the edge of each MODFLOW row. yedge is dimensioned
-        to NROW + 1. If yedge is not a numpy.ndarray it is converted to a
-        numpy.ndarray.
-    returnvertices: bool
-        Return the original polyline vertices in the list of numpy.ndarray
-        containing vertices resulting from intersection of the provided
-        polygon and the MODFLOW model grid if returnvertices=True.
-        (default is False).
-
-    Returns
-    -------
-    (x, y, dlen) : numpy.ndarray of tuples
-        numpy.ndarray of tuples containing the x, y, and segment length of the
-        intersection of the provided polyline with the rectilinear MODFLOW
-        grid.
-
-    Examples
-    --------
-    >>> import flopy
-    >>> ptsout = flopy.plotutil.line_intersect_grid(ptsin, xedge, yedge)
-
-    """
-
-    small_value = 1.0e-4
-
-    # make sure xedge and yedge are numpy arrays
-    if not isinstance(xedge, np.ndarray):
-        xedge = np.array(xedge)
-    if not isinstance(yedge, np.ndarray):
-        yedge = np.array(yedge)
-
-    # build list of points along current line
-    pts = []
-    npts = len(ptsin)
-    dlen = 0.0
-    for idx in range(1, npts):
-        x0 = ptsin[idx - 1][0]
-        x1 = ptsin[idx][0]
-        y0 = ptsin[idx - 1][1]
-        y1 = ptsin[idx][1]
-        a = x1 - x0
-        b = y1 - y0
-        c = math.sqrt(math.pow(a, 2.0) + math.pow(b, 2.0))
-        # find cells with (x0, y0) and (x1, y1)
-        irow0, jcol0 = findrowcolumn((x0, y0), xedge, yedge)
-        irow1, jcol1 = findrowcolumn((x1, y1), xedge, yedge)
-        # determine direction to go in the x- and y-directions
-        jx = 0
-        incx = abs(small_value * a / c)
-        iy = 0
-        incy = -abs(small_value * b / c)
-        if a == 0.0:
-            incx = 0.0
-        # go to the right
-        elif a > 0.0:
-            jx = 1
-            incx *= -1.0
-        if b == 0.0:
-            incy = 0.0
-        # go down
-        elif b < 0.0:
-            iy = 1
-            incy *= -1.0
-        # process data
-        if irow0 >= 0 and jcol0 >= 0:
-            iadd = True
-            if idx > 1 and returnvertices:
-                iadd = False
-            if iadd:
-                pts.append((x0, y0, dlen))
-        icnt = 0
-        while True:
-            icnt += 1
-            dx = xedge[jcol0 + jx] - x0
-            dlx = 0.0
-            if a != 0.0:
-                dlx = c * dx / a
-            dy = yedge[irow0 + iy] - y0
-            dly = 0.0
-            if b != 0.0:
-                dly = c * dy / b
-            if dlx != 0.0 and dly != 0.0:
-                if abs(dlx) < abs(dly):
-                    dy = dx * b / a
-                else:
-                    dx = dy * a / b
-            xt = x0 + dx + incx
-            yt = y0 + dy + incy
-            dl = math.sqrt(math.pow((xt - x0), 2.0) + math.pow((yt - y0), 2.0))
-            dlen += dl
-            if not returnvertices:
-                pts.append((xt, yt, dlen))
-            x0, y0 = xt, yt
-            xt = x0 - 2.0 * incx
-            yt = y0 - 2.0 * incy
-            dl = math.sqrt(math.pow((xt - x0), 2.0) + math.pow((yt - y0), 2.0))
-            dlen += dl
-            x0, y0 = xt, yt
-            irow0, jcol0 = findrowcolumn((x0, y0), xedge, yedge)
-            if irow0 >= 0 and jcol0 >= 0:
-                if not returnvertices:
-                    pts.append((xt, yt, dlen))
-            elif irow1 < 0 or jcol1 < 0:
-                dl = math.sqrt(
-                    math.pow((x1 - x0), 2.0) + math.pow((y1 - y0), 2.0)
-                )
-                dlen += dl
-                break
-            if irow0 == irow1 and jcol0 == jcol1:
-                dl = math.sqrt(
-                    math.pow((x1 - x0), 2.0) + math.pow((y1 - y0), 2.0)
-                )
-                dlen += dl
-                pts.append((x1, y1, dlen))
-                break
-    return np.array(pts)
-
-
-def cell_value_points(pts, xedge, yedge, vdata):
-    """
-    Intersect a list of polyline vertices with a rectilinear MODFLOW
-    grid. Vertices at the intersection of the polyline with the grid
-    cell edges is returned. Optionally the original polyline vertices
-    are returned.
-
-    Parameters
-    ----------
-    pts : list
-        A list of x, y points and polyline length to extract defining the
-        vertices of a polyline that
-    xedge : numpy.ndarray
-        x-coordinate of the edge of each MODFLOW column. The shape of xedge is
-        (NCOL + 1). If xedge is not a numpy.ndarray it is converted to a
-        numpy.ndarray.
-    yedge : numpy.ndarray
-        y-coordinate of the edge of each MODFLOW row. The shape of yedge is
-        (NROW + 1). If yedge is not a numpy.ndarray it is converted to a
-        numpy.ndarray.
-    vdata : numpy.ndarray
-        Data (i.e., head, hk, etc.) for a rectilinear MODFLOW model grid. The
-        shape of vdata is (NROW, NCOL). If vdata is not a numpy.ndarray it is
-        converted to a numpy.ndarray.
-
-    Returns
-    -------
-    vcell : numpy.ndarray
-        numpy.ndarray of of data values from the vdata numpy.ndarray at x- and
-        y-coordinate locations in pts.
-
-    Examples
-    --------
-    >>> import flopy
-    >>> vcell = flopy.plotutil.cell_value_points(xpts, xedge, yedge, head[0, :, :])
-
-    """
-
-    # make sure xedge and yedge are numpy arrays
-    if not isinstance(xedge, np.ndarray):
-        xedge = np.array(xedge)
-    if not isinstance(yedge, np.ndarray):
-        yedge = np.array(yedge)
-    if not isinstance(vdata, np.ndarray):
-        vdata = np.array(vdata)
-
-    vcell = []
-    for (xt, yt, _) in pts:
-        # find the modflow cell containing point
-        irow, jcol = findrowcolumn((xt, yt), xedge, yedge)
-        if irow >= 0 and jcol >= 0:
-            if np.isnan(vdata[irow, jcol]):
-                vcell.append(np.nan)
-            else:
-                v = np.asarray(vdata[irow, jcol])
-                vcell.append(v)
-
-    return np.array(vcell)
-
-
 def _set_coord_info(mg, xul, yul, xll, yll, rotation):
     """
 
@@ -2722,7 +2583,9 @@ def _depreciated_dis_handler(modelgrid, dis):
     import warnings
 
     warnings.warn(
-        "the dis parameter has been depreciated.", PendingDeprecationWarning
+        "the dis parameter has been depreciated and will be removed in "
+        "version 3.3.5.",
+        PendingDeprecationWarning,
     )
     if modelgrid.grid_type == "vertex":
         modelgrid = VertexGrid(
@@ -2795,3 +2658,303 @@ def advanced_package_bc_helper(pkg, modelgrid, kper):
             "Pkg {} not implemented for bc plotting".format(pkg.package_type)
         )
     return idx
+
+
+def filter_modpath_by_travel_time(recarray, travel_time):
+    """
+
+    :param recarray:
+    :param travel_time:
+    :return:
+    """
+    if travel_time is None:
+        tp = recarray.copy()
+    else:
+        if isinstance(travel_time, str):
+            funcs = {
+                "<=": lambda a, b: a["time"] <= b,
+                ">=": lambda a, b: a["time"] >= b,
+                "<": lambda a, b: a["time"] < b,
+                ">": lambda a, b: a["time"] > b,
+            }
+            idx = None
+            for k, func in sorted(funcs.items())[::-1]:
+                if k in travel_time:
+                    time = float(travel_time.replace(k, ""))
+                    idx = func(recarray, time)
+                    break
+            if idx is None:
+                try:
+                    time = float(travel_time)
+                    idx = recarray["time"] <= time
+                except (ValueError, KeyError):
+                    errmsg = (
+                        "flopy.map.plot_pathline travel_time "
+                        + "variable cannot be parsed. "
+                        + "Acceptable logical variables are , "
+                        + "<=, <, >=, and >. "
+                        + "You passed {}".format(travel_time)
+                    )
+                    raise Exception(errmsg)
+        else:
+            time = float(travel_time)
+            idx = recarray["time"] <= time
+        tp = recarray[idx]
+
+    return tp
+
+
+def intersect_modpath_with_crosssection(
+    recarrays,
+    projpts,
+    xvertices,
+    yvertices,
+    projection,
+    method="cell",
+    starting=False,
+):
+    """
+
+    :param recarrays:
+    :param projpts:
+    :param xvertices:
+    :param yvertices:
+    :param projection:
+    :param method:
+    :param starting:
+    :return:
+    """
+
+    from ..utils.geometry import point_in_polygon
+
+    xp, yp, zp = "x", "y", "z"
+    if starting:
+        xp, yp, zp = "x0", "y0", "z0"
+
+    if not isinstance(recarrays, list):
+        recarrays = [
+            recarrays,
+        ]
+
+    if projection == "x":
+        v_opp = yvertices
+        v_norm = xvertices
+        oprj = yp
+        prj = xp
+    else:
+        v_opp = xvertices
+        v_norm = yvertices
+        oprj = xp
+        prj = yp
+
+    # set points opposite projection direction
+    oppts = {}
+    nppts = {}
+
+    for cell, verts in projpts.items():
+        zmin = np.min(np.array(verts)[:, 1])
+        zmax = np.max(np.array(verts)[:, 1])
+        nmin = np.min(v_norm[cell])
+        nmax = np.max(v_norm[cell])
+        omin = np.min(v_opp[cell])
+        omax = np.max(v_opp[cell])
+        oppts[cell] = np.array(
+            [
+                [omin, zmax],
+                [omax, zmax],
+                [omax, zmin],
+                [omin, zmin],
+                [omin, zmax],
+            ]
+        )
+
+        # intersects w/actual...
+        nppts[cell] = np.array(
+            [
+                [nmin, zmax],
+                [nmax, zmax],
+                [nmax, zmin],
+                [nmin, zmin],
+                [nmin, zmax],
+            ]
+        )
+
+    idict = {}
+    for recarray in recarrays:
+        for cell, _ in projpts.items():
+            m0 = point_in_polygon(
+                recarray[prj].reshape(1, -1),
+                recarray[zp].reshape(1, -1),
+                nppts[cell],
+            )
+            if method == "cell":
+                m1 = point_in_polygon(
+                    recarray[oprj].reshape(1, -1),
+                    recarray[zp].reshape(1, -1),
+                    oppts[cell],
+                )
+                idx = [
+                    i
+                    for i, (x, y) in enumerate(zip(m0[0], m1[0]))
+                    if x == y == True
+                ]
+            else:
+                idx = [i for i, x in enumerate(m0[0]) if x == True]
+
+            if idx:
+                if cell not in idict:
+                    idict[cell] = [recarray[idx]]
+                else:
+                    idict[cell].append(recarray[idx])
+
+    return idict
+
+
+def reproject_modpath_to_crosssection(
+    idict,
+    projpts,
+    xypts,
+    projection,
+    modelgrid,
+    geographic_coords,
+    starting=False,
+):
+    """
+
+    :param idict:
+    :param projpts:
+    :param xypts:
+    :param projection:
+    :param modelgrid:
+    :param geographic_coords:
+    :param starting:
+    :return:
+    """
+    from ..utils import geometry
+
+    xp, yp, zp = "x", "y", "z"
+    if starting:
+        xp, yp, zp = "x0", "y0", "z0"
+
+    proj = xp
+    if projection == "y":
+        proj = yp
+
+    ptdict = {}
+    if not geographic_coords:
+        for cell, recarrays in idict.items():
+            line = xypts[cell]
+            if projection == "x":
+                d0 = np.min([i[0] for i in projpts[cell]])
+            else:
+                d0 = np.max([i[0] for i in projpts[cell]])
+            for rec in recarrays:
+                pts = list(zip(rec[xp], rec[yp]))
+                x, y = geometry.project_point_onto_xc_line(
+                    line, pts, d0, projection
+                )
+                rec[xp] = x
+                rec[yp] = y
+                pid = rec["particleid"][0]
+                pline = list(zip(rec[proj], rec[zp]))
+                if pid not in ptdict:
+                    ptdict[pid] = pline
+                else:
+                    ptdict[pid] += pline
+    else:
+        for cell, recarrays in idict.items():
+            for rec in recarrays:
+                x, y = geometry.transform(
+                    rec[xp],
+                    rec[yp],
+                    modelgrid.xoffset,
+                    modelgrid.yoffset,
+                    modelgrid.angrot_radians,
+                )
+                rec[xp] = x
+                rec[yp] = y
+                pid = rec["particleid"][0]
+                pline = list(zip(rec[proj], rec[zp]))
+                if pid not in ptdict:
+                    ptdict[pid] = pline
+                else:
+                    ptdict[pid] += pline
+
+    return ptdict
+
+
+def parse_modpath_selection_options(
+    ep,
+    direction,
+    selection,
+    selection_direction,
+):
+    """
+
+    :return:
+    """
+    ep = ep.copy()
+    direction = direction.lower()
+    if direction == "starting":
+        istart = True
+        xp, yp = "x0", "y0"
+
+    else:
+        istart = False
+        direction = "ending"
+        xp, yp = "x", "y"
+
+    if selection_direction is not None:
+        selection_direction = selection_direction.lower()
+        if selection_direction != "starting":
+            selection_direction = "ending"
+
+    else:
+        if direction.lower() == "starting":
+            selection_direction = "ending"
+        elif direction.lower() == "ending":
+            selection_direction = "starting"
+
+    # selection of endpoints
+    if selection is not None:
+        if isinstance(selection, int):
+            selection = tuple((selection,))
+        try:
+            if len(selection) == 1:
+                node = selection[0]
+                if selection_direction.lower() == "starting":
+                    nsel = "node0"
+                else:
+                    nsel = "node"
+                # make selection
+                idx = ep[nsel] == node
+                tep = ep[idx]
+            elif len(selection) == 3:
+                k, i, j = selection[0], selection[1], selection[2]
+                if selection_direction.lower() == "starting":
+                    ksel, isel, jsel = "k0", "i0", "j0"
+                else:
+                    ksel, isel, jsel = "k", "i", "j"
+                # make selection
+                idx = (ep[ksel] == k) & (ep[isel] == i) & (ep[jsel] == j)
+                tep = ep[idx]
+            else:
+                errmsg = (
+                    "plot_endpoint selection must be "
+                    + "a zero-based layer, row, column tuple "
+                    + "(l, r, c) or node number (MODPATH 7) of "
+                    + "the location to evaluate (i.e., well location)."
+                )
+                raise Exception(errmsg)
+        except (ValueError, KeyError, IndexError):
+            errmsg = (
+                "plot_endpoint selection must be a "
+                + "zero-based layer, row, column tuple (l, r, c) "
+                + "or node number (MODPATH 7) of the location "
+                + "to evaluate (i.e., well location)."
+            )
+            raise Exception(errmsg)
+    else:
+        tep = ep.copy()
+
+    return tep, istart, xp, yp

flopy/plot/styles.py

@@ -0,0 +1,467 @@
+try:
+    import matplotlib.pyplot as plt
+    import matplotlib as mpl
+except (ImportError, ModuleNotFoundError):
+    plt = None
+
+import os
+import platform
+
+
+class styles(object):
+    """Styles class for custom matplotlib styling
+
+    The class contains both custom styles and plotting methods
+    for custom formatting using a specific matplotlib style
+
+    Additional styles can be easily added to the mplstyle folder and
+    accessed using the plt.style.context() method.
+
+    """
+
+    _ws = os.path.abspath(os.path.dirname(__file__))
+    _map_style = os.path.join(_ws, "mplstyle", "usgsmap.mplstyle")
+    _plot_style = os.path.join(_ws, "mplstyle", "usgsplot.mplstyle")
+    if platform.system() == "linux":
+        _map_style = os.path.join(_ws, "mplstyle", "usgsmap_linux.mplstyle")
+        _plot_style = os.path.join(_ws, "mplstyle", "usgsplot_linux.mplstyle")
+
+    @classmethod
+    def USGSMap(cls):
+        return plt.style.context(styles._map_style)
+
+    @classmethod
+    def USGSPlot(cls):
+        return plt.style.context(styles._map_style)
+
+    @classmethod
+    def set_font_type(cls, family, fontname):
+        """
+        Method to set the matplotlib font type for the current style
+
+        Note: this method only works when adding text using the styles
+        methods.
+
+        Parameters
+        ----------
+        family : str
+            matplotlib.rcparams font.family
+        font : str
+            matplotlib.rcparams font.fontname
+
+        Returns
+        -------
+            None
+        """
+        mpl.rcParams["font.family"] = family
+        mpl.rcParams["font." + family] = fontname
+        return mpl.rcParams
+
+    @classmethod
+    def heading(
+        self,
+        ax=None,
+        letter=None,
+        heading=None,
+        x=0.00,
+        y=1.01,
+        idx=None,
+        fontsize=9,
+    ):
+        """Add a USGS-style heading to a matplotlib axis object
+
+        Parameters
+        ----------
+        ax : axis object
+            matplotlib axis object (default is None)
+        letter : str
+            string that defines the subplot (A, B, C, etc.)
+        heading : str
+            text string
+        x : float
+            location of the heading in the x-direction in normalized plot
+            dimensions ranging from 0 to 1 (default is 0.00)
+        y : float
+            location of the heading in the y-direction in normalized plot
+            dimensions ranging from 0 to 1 (default is 1.01)
+        idx : int
+            index for programatically generating the heading letter when letter
+            is None and idx is not None. idx = 0 will generate A
+            (default is None)
+
+        Returns
+        -------
+        text : object
+            matplotlib text object
+
+        """
+        if ax is None:
+            ax = plt.gca()
+
+        if letter is None and idx is not None:
+            letter = chr(ord("A") + idx)
+
+        font = styles.__set_fontspec(
+            bold=True, italic=False, fontsize=fontsize
+        )
+
+        if letter is not None:
+            if heading is None:
+                text = letter.replace(".", "")
+            else:
+                letter = letter.rstrip()
+                if not letter.endswith("."):
+                    letter += "."
+                text = letter + " " + heading
+        else:
+            text = heading
+
+        if text is None:
+            return
+
+        text = ax.text(
+            x,
+            y,
+            text,
+            va="bottom",
+            ha="left",
+            fontdict=font,
+            transform=ax.transAxes,
+        )
+        return text
+
+    @classmethod
+    def xlabel(cls, ax=None, label="", bold=False, italic=False, **kwargs):
+        """Method to set the xlabel using the styled fontdict
+
+        Parameters
+        ----------
+        ax : axis object
+            matplotlib axis object (default is None)
+        label : str
+            axis label for the chart
+        bold : bool
+            flag to switch to boldface test
+        italic : bool
+            flag to use italic text
+        kwargs : dict
+            keyword arguments for the matplotlib set_xlabel method
+
+        Returns
+        -------
+            None
+        """
+        if ax is None:
+            ax = plt.gca()
+        fontsize = kwargs.pop("fontsize", 9)
+        fontspec = styles.__set_fontspec(
+            bold=bold, italic=italic, fontsize=fontsize
+        )
+        ax.set_xlabel(label, fontdict=fontspec, **kwargs)
+
+    @classmethod
+    def ylabel(cls, ax=None, label="", bold=False, italic=False, **kwargs):
+        """Method to set the ylabel using the styled fontdict
+
+        Parameters
+        ----------
+        ax : axis object
+            matplotlib axis object (default is None)
+        label : str
+            axis label for the chart
+        bold : bool
+            flag to switch to boldface test
+        italic : bool
+            flag to use italic text
+        kwargs : dict
+            keyword arguments for the matplotlib set_xlabel method
+
+        Returns
+        -------
+            None
+        """
+        if ax is None:
+            ax = plt.gca()
+
+        fontsize = kwargs.pop("fontsize", 9)
+        fontspec = styles.__set_fontspec(
+            bold=bold, italic=italic, fontsize=fontsize
+        )
+        ax.set_ylabel(label, fontdict=fontspec, **kwargs)
+
+    @classmethod
+    def graph_legend(cls, ax=None, handles=None, labels=None, **kwargs):
+        """Add a USGS-style legend to a matplotlib axis object
+
+        Parameters
+        ----------
+        ax : axis object
+            matplotlib axis object (default is None)
+        handles : list
+            list of legend handles
+        labels : list
+            list of labels for legend handles
+        kwargs : kwargs
+            matplotlib legend kwargs
+
+        Returns
+        -------
+        leg : object
+            matplotlib legend object
+
+        """
+        if ax is None:
+            ax = plt.gca()
+
+        fontspec = styles.__set_fontspec(bold=True, italic=False, family=True)
+
+        if handles is None or labels is None:
+            handles, labels = ax.get_legend_handles_labels()
+        leg = ax.legend(handles, labels, prop=fontspec, **kwargs)
+
+        # add title to legend
+        if "title" in kwargs:
+            title = kwargs.pop("title")
+        else:
+            title = None
+        leg = styles.graph_legend_title(leg, title=title)
+        return leg
+
+    @classmethod
+    def graph_legend_title(cls, leg, title=None):
+        """Set the legend title for a matplotlib legend object
+
+        Parameters
+        ----------
+        leg : legend object
+            matplotlib legend object
+        title : str
+            title for legend
+
+        Returns
+        -------
+        leg : object
+            matplotlib legend object
+
+        """
+        if title is None:
+            title = "EXPLANATION"
+        elif title.lower() == "none":
+            title = None
+
+        fontspec = styles.__set_fontspec(bold=True, italic=False, family=True)
+
+        leg.set_title(title, prop=fontspec)
+        return leg
+
+    @classmethod
+    def add_text(
+        cls,
+        ax=None,
+        text="",
+        x=0.0,
+        y=0.0,
+        transform=True,
+        bold=True,
+        italic=True,
+        fontsize=9,
+        ha="left",
+        va="bottom",
+        **kwargs
+    ):
+        """Add USGS-style text to a axis object
+
+        Parameters
+        ----------
+        ax : axis object
+            matplotlib axis object (default is None)
+        text : str
+            text string
+        x : float
+            x-location of text string (default is 0.)
+        y : float
+            y-location of text string (default is 0.)
+        transform : bool
+            boolean that determines if a transformed (True) or data (False)
+            coordinate system is used to define the (x, y) location of the
+            text string (default is True)
+        bold : bool
+            boolean indicating if bold font (default is True)
+        italic : bool
+            boolean indicating if italic font (default is True)
+        fontsize : int
+            font size (default is 9 points)
+        ha : str
+            matplotlib horizontal alignment keyword (default is left)
+        va : str
+            matplotlib vertical alignment keyword (default is bottom)
+        kwargs : dict
+            dictionary with valid matplotlib text object keywords
+
+        Returns
+        -------
+        text_obj : object
+            matplotlib text object
+
+        """
+        if ax is None:
+            ax = plt.gca()
+
+        if transform:
+            transform = ax.transAxes
+        else:
+            transform = ax.transData
+
+        font = styles.__set_fontspec(
+            bold=bold, italic=italic, fontsize=fontsize
+        )
+
+        text_obj = ax.text(
+            x,
+            y,
+            text,
+            va=va,
+            ha=ha,
+            fontdict=font,
+            transform=transform,
+            **kwargs
+        )
+        return text_obj
+
+    @classmethod
+    def add_annotation(
+        cls,
+        ax=None,
+        text="",
+        xy=None,
+        xytext=None,
+        bold=True,
+        italic=True,
+        fontsize=9,
+        ha="left",
+        va="bottom",
+        **kwargs
+    ):
+        """Add an annotation to a axis object
+
+        Parameters
+        ----------
+        ax : axis object
+            matplotlib axis object (default is None)
+        text : str
+            text string
+        xy : tuple
+            tuple with the location of the annotation (default is None)
+        xytext : tuple
+            tuple with the location of the text
+        bold : bool
+            boolean indicating if bold font (default is True)
+        italic : bool
+            boolean indicating if italic font (default is True)
+        fontsize : int
+            font size (default is 9 points)
+        ha : str
+            matplotlib horizontal alignment keyword (default is left)
+        va : str
+            matplotlib vertical alignment keyword (default is bottom)
+        kwargs : dict
+            dictionary with valid matplotlib annotation object keywords
+
+        Returns
+        -------
+        ann_obj : object
+            matplotlib annotation object
+
+        """
+        if ax is None:
+            ax = plt.gca()
+
+        if xy is None:
+            xy = (0.0, 0.0)
+
+        if xytext is None:
+            xytext = (0.0, 0.0)
+
+        fontspec = styles.__set_fontspec(
+            bold=bold, italic=italic, fontsize=fontsize
+        )
+        # add font information to kwargs
+        if kwargs is None:
+            kwargs = fontspec
+        else:
+            for key, value in fontspec.items():
+                kwargs[key] = value
+
+        # create annotation
+        ann_obj = ax.annotate(text, xy, xytext, va=va, ha=ha, **kwargs)
+
+        return ann_obj
+
+    @classmethod
+    def remove_edge_ticks(cls, ax=None):
+        """Remove unnecessary ticks on the edges of the plot
+
+        Parameters
+        ----------
+        ax : axis object
+            matplotlib axis object (default is None)
+
+        Returns
+        -------
+        ax : axis object
+            matplotlib axis object
+
+        """
+        if ax is None:
+            ax = plt.gca()
+
+        ax.tick_params(axis="both", which="both", length=0)
+
+    @classmethod
+    def __set_fontspec(cls, bold=True, italic=True, fontsize=9, family=False):
+        """Create fontspec dictionary for matplotlib pyplot objects
+
+        Parameters
+        ----------
+        bold : bool
+            boolean indicating if font is bold (default is True)
+        italic : bool
+            boolean indicating if font is italic (default is True)
+        fontsize : int
+            font size (default is 9 point)
+
+
+        Returns
+        -------
+            dict
+        """
+        family = mpl.rcParams["font.family"][0]
+        font = mpl.rcParams["font." + family][0]
+
+        if bold:
+            weight = "bold"
+        else:
+            weight = "normal"
+
+        if italic:
+            style = "italic"
+        else:
+            style = "normal"
+
+        # define fontspec dictionary
+        fontspec = {
+            "fontname": font,
+            "size": fontsize,
+            "weight": weight,
+            "style": style,
+        }
+
+        if family:
+            fontspec.pop("fontname")
+            fontspec["family"] = family
+
+        return fontspec
+
+
+if plt is None:
+    styles = None

flopy/plot/vcrosssection.py

@@ -1,817 +0,0 @@
-import numpy as np
-
-try:
-    import matplotlib.pyplot as plt
-except:
-    plt = None
-from flopy.plot import plotutil
-from flopy.utils import geometry
-from flopy.plot.crosssection import _CrossSection
-import warnings
-
-warnings.filterwarnings("ignore", category=RuntimeWarning)
-
-
-class _VertexCrossSection(_CrossSection):
-    """
-    Class to create a cross section of the model from a vertex
-    discretization.
-
-    Class is not to be instantiated by the user!
-
-    Parameters
-    ----------
-    ax : matplotlib.pyplot axis
-        The plot axis.  If not provided it, plt.gca() will be used.
-    model : flopy.modflow object
-        flopy model object. (Default is None)
-    modelgrid : flopy.discretization.VertexGrid
-        Vertex model grid object
-    line : dict
-        Dictionary with either "row", "column", or "line" key. If key
-        is "row" or "column" key value should be the zero-based row or
-        column index for cross-section. If key is "line" value should
-        be an array of (x, y) tuples with vertices of cross-section.
-        Vertices should be in map coordinates consistent with xul,
-        yul, and rotation.
-    extent : tuple of floats
-        (xmin, xmax, ymin, ymax) will be used to specify axes limits.  If None
-        then these will be calculated based on grid, coordinates, and rotation.
-    geographic_coords : bool
-        boolean flag to allow the user to plot cross section lines in
-        geographic coordinates. If False (default), cross section is plotted
-        as the distance along the cross section line.
-
-    """
-
-    def __init__(
-        self,
-        ax=None,
-        model=None,
-        modelgrid=None,
-        line=None,
-        extent=None,
-        geographic_coords=False,
-    ):
-        super(_VertexCrossSection, self).__init__(
-            ax=ax,
-            model=model,
-            modelgrid=modelgrid,
-            geographic_coords=geographic_coords,
-        )
-
-        if line is None:
-            err_msg = "line must be specified."
-            raise Exception(err_msg)
-
-        linekeys = [linekeys.lower() for linekeys in list(line.keys())]
-
-        if len(linekeys) != 1:
-            err_msg = (
-                "Either row, column, or line must be specified "
-                "in line dictionary.\nkeys specified: "
-            )
-            for k in linekeys:
-                err_msg += "{} ".format(k)
-            raise Exception(err_msg)
-
-        elif "line" not in linekeys:
-            err_msg = (
-                "only line can be specified in line dictionary "
-                "for vertex Discretization"
-            )
-            raise AssertionError(err_msg)
-
-        onkey = linekeys[0]
-
-        if ax is None:
-            self.ax = plt.gca()
-        else:
-            self.ax = ax
-
-        self.direction = "xy"
-        # convert pts list to a numpy array
-        verts = line[onkey]
-        xp = []
-        yp = []
-        for [v1, v2] in verts:
-            xp.append(v1)
-            yp.append(v2)
-
-        # unrotate and untransform modelgrid into modflow coordinates!
-        xp, yp = geometry.transform(
-            xp,
-            yp,
-            self.mg.xoffset,
-            self.mg.yoffset,
-            self.mg.angrot_radians,
-            inverse=True,
-        )
-
-        self.xcellcenters, self.ycellcenters = geometry.transform(
-            self.mg.xcellcenters,
-            self.mg.ycellcenters,
-            self.mg.xoffset,
-            self.mg.yoffset,
-            self.mg.angrot_radians,
-            inverse=True,
-        )
-
-        try:
-            self.xvertices, self.yvertices = geometry.transform(
-                self.mg.xvertices,
-                self.mg.yvertices,
-                self.mg.xoffset,
-                self.mg.yoffset,
-                self.mg.angrot_radians,
-                inverse=True,
-            )
-        except ValueError:
-            # irregular shapes in vertex grid ie. squares and triangles
-            (
-                xverts,
-                yverts,
-            ) = plotutil.UnstructuredPlotUtilities.irregular_shape_patch(
-                self.mg.xvertices, self.mg.yvertices
-            )
-
-            self.xvertices, self.yvertices = geometry.transform(
-                xverts,
-                yverts,
-                self.mg.xoffset,
-                self.mg.yoffset,
-                self.mg.angrot_radians,
-                inverse=True,
-            )
-
-        pts = [(xt, yt) for xt, yt in zip(xp, yp)]
-        self.pts = np.array(pts)
-
-        # get points along the line
-
-        self.xypts = plotutil.UnstructuredPlotUtilities.line_intersect_grid(
-            self.pts, self.xvertices, self.yvertices
-        )
-
-        if len(self.xypts) < 2:
-            s = "cross-section cannot be created\n."
-            s += "   less than 2 points intersect the model grid\n"
-            s += "   {} points intersect the grid.".format(len(self.xypts))
-            raise Exception(s)
-
-        if self.geographic_coords:
-            # transform back to geographic coordinates
-            xypts = {}
-            for nn, pt in self.xypts.items():
-                xp = [t[0] for t in pt]
-                yp = [t[1] for t in pt]
-                xp, yp = geometry.transform(
-                    xp,
-                    yp,
-                    self.mg.xoffset,
-                    self.mg.yoffset,
-                    self.mg.angrot_radians,
-                )
-                xypts[nn] = [(xt, yt) for xt, yt in zip(xp, yp)]
-
-            self.xypts = xypts
-
-        top = self.mg.top
-        top.shape = (1, -1)
-        botm = self.mg.botm
-        nlay = len(botm)
-        ncpl = self.mg.ncpl
-
-        elev = list(top.copy())
-        for k in range(nlay):
-            elev.append(botm[k, :])
-
-        self.elev = np.array(elev)
-
-        self.idomain = self.mg.idomain
-        if self.mg.idomain is None:
-            self.idomain = np.ones((nlay, ncpl), dtype=int)
-
-        # choose a projection direction based on maximum information
-        xpts = []
-        ypts = []
-        for nn, verts in self.xypts.items():
-            for v in verts:
-                xpts.append(v[0])
-                ypts.append(v[1])
-
-        if np.max(xpts) - np.min(xpts) > np.max(ypts) - np.min(ypts):
-            self.direction = "x"
-        else:
-            self.direction = "y"
-
-        # make vertex array based on projection direction
-        self.projpts = self.set_zpts(None)
-
-        # Create cross-section extent
-        if extent is None:
-            self.extent = self.get_extent()
-        else:
-            self.extent = extent
-
-        self.layer0 = None
-        self.layer1 = None
-
-        self.d = {
-            i: (np.min(np.array(v).T[0]), np.max(np.array(v).T[0]))
-            for i, v in sorted(self.projpts.items())
-        }
-
-        self.xpts = None
-        self.active = None
-        self.ncb = None
-        self.laycbd = None
-        self.zpts = None
-        self.xcentergrid = None
-        self.zcentergrid = None
-        self.geographic_xcentergrid = None
-        self.geographic_xpts = None
-
-        # Set axis limits
-        self.ax.set_xlim(self.extent[0], self.extent[1])
-        self.ax.set_ylim(self.extent[2], self.extent[3])
-
-    def plot_array(self, a, masked_values=None, head=None, **kwargs):
-        """
-        Plot a three-dimensional array as a patch collection.
-
-        Parameters
-        ----------
-        a : numpy.ndarray
-            Three-dimensional array to plot.
-        masked_values : iterable of floats, ints
-            Values to mask.
-        head : numpy.ndarray
-            Three-dimensional array to set top of patches to the minimum
-            of the top of a layer or the head value. Used to create
-            patches that conform to water-level elevations.
-        **kwargs : dictionary
-            keyword arguments passed to matplotlib.collections.PatchCollection
-
-        Returns
-        -------
-        patches : matplotlib.collections.PatchCollection
-
-        """
-        if "ax" in kwargs:
-            ax = kwargs.pop("ax")
-        else:
-            ax = self.ax
-
-        if not isinstance(a, np.ndarray):
-            a = np.array(a)
-
-        if a.ndim > 1:
-            a = np.ravel(a)
-
-        if masked_values is not None:
-            for mval in masked_values:
-                a = np.ma.masked_values(a, mval)
-
-        if isinstance(head, np.ndarray):
-            projpts = self.set_zpts(np.ravel(head))
-        else:
-            projpts = self.projpts
-
-        pc = self.get_grid_patch_collection(projpts, a, **kwargs)
-        if pc is not None:
-            ax.add_collection(pc)
-            ax.set_xlim(self.extent[0], self.extent[1])
-            ax.set_ylim(self.extent[2], self.extent[3])
-
-        return pc
-
-    def plot_surface(self, a, masked_values=None, **kwargs):
-        """
-        Plot a two- or three-dimensional array as line(s).
-
-        Parameters
-        ----------
-        a : numpy.ndarray
-            Two- or three-dimensional array to plot.
-        masked_values : iterable of floats, ints
-            Values to mask.
-        **kwargs : dictionary
-            keyword arguments passed to matplotlib.pyplot.plot
-
-        Returns
-        -------
-        plot : list containing matplotlib.plot objects
-        """
-        if "ax" in kwargs:
-            ax = kwargs.pop("ax")
-        else:
-            ax = self.ax
-
-        if "color" in kwargs:
-            color = kwargs.pop("color")
-        elif "c" in kwargs:
-            color = kwargs.pop("c")
-        else:
-            color = "b"
-
-        if not isinstance(a, np.ndarray):
-            a = np.array(a)
-
-        if a.ndim > 1:
-            a = np.ravel(a)
-
-        if a.size % self.mg.ncpl != 0:
-            raise AssertionError("Array size must be a multiple of ncpl")
-
-        if masked_values is not None:
-            for mval in masked_values:
-                a = np.ma.masked_values(a, mval)
-
-        data = []
-        lay_data = []
-        d = []
-        lay_d = []
-        dim = self.mg.ncpl
-        for cell, verts in sorted(self.projpts.items()):
-
-            if cell >= a.size:
-                continue
-            elif np.isnan(a[cell]):
-                continue
-            elif a[cell] is np.ma.masked:
-                continue
-
-            if cell >= dim:
-                data.append(lay_data)
-                d.append(lay_d)
-                dim += self.mg.ncpl
-                lay_data = [(a[cell], a[cell])]
-                lay_d = [self.d[cell]]
-            else:
-                lay_data.append((a[cell], a[cell]))
-                lay_d.append(self.d[cell])
-
-        if lay_data:
-            data.append(lay_data)
-            d.append(lay_d)
-
-        data = np.array(data)
-        d = np.array(d)
-
-        plot = []
-        for k in range(data.shape[0]):
-            if ax is None:
-                ax = plt.gca()
-            for ix, _ in enumerate(data[k]):
-                ax.plot(d[k, ix], data[k, ix], color=color, **kwargs)
-
-            ax.set_xlim(self.extent[0], self.extent[1])
-            ax.set_ylim(self.extent[2], self.extent[3])
-            plot.append(ax)
-
-        return plot
-
-    def plot_fill_between(
-        self,
-        a,
-        colors=("blue", "red"),
-        masked_values=None,
-        head=None,
-        **kwargs
-    ):
-        """
-        Plot a three-dimensional array as lines.
-
-        Parameters
-        ----------
-        a : numpy.ndarray
-            Three-dimensional array to plot.
-        colors: list
-            matplotlib fill colors, two required
-        masked_values : iterable of floats, ints
-            Values to mask.
-        head : numpy.ndarray
-            Three-dimensional array to set top of patches to the minimum
-            of the top of a layer or the head value. Used to create
-            patches that conform to water-level elevations.
-        **kwargs : dictionary
-            keyword arguments passed to matplotlib.pyplot.plot
-
-        Returns
-        -------
-        plot : list containing matplotlib.fillbetween objects
-
-        """
-        if "ax" in kwargs:
-            ax = kwargs.pop("ax")
-        else:
-            ax = self.ax
-
-        if not isinstance(a, np.ndarray):
-            a = np.array(a)
-
-        a = np.ravel(a)
-
-        if masked_values is not None:
-            for mval in masked_values:
-                a = np.ma.masked_values(a, mval)
-
-        if isinstance(head, np.ndarray):
-            projpts = self.set_zpts(head)
-        else:
-            projpts = self.projpts
-
-        plot = []
-        for cell, verts in sorted(projpts.items()):
-            if cell >= a.size:
-                continue
-            elif np.isnan(a[cell]):
-                continue
-            elif a[cell] is np.ma.masked:
-                continue
-
-            x = list(set(np.array(verts.T[0])))
-            y1 = np.max(np.array(verts.T[1]))
-            y2 = np.min(np.array(verts.T[1]))
-            v = a[cell]
-
-            if v > y1:
-                v = y1
-
-            elif v < y2:
-                v = y2
-
-            v = [v] * len(x)
-
-            x = np.array(x)
-            plot.append(ax.fill_between(x, y1, v, color=colors[0], **kwargs))
-            plot.append(ax.fill_between(x, v, y2, color=colors[1], **kwargs))
-
-        return plot
-
-    def contour_array(self, a, masked_values=None, head=None, **kwargs):
-        """
-        Contour a two-dimensional array.
-
-        Parameters
-        ----------
-        a : numpy.ndarray
-            Three-dimensional array to plot.
-        masked_values : iterable of floats, ints
-            Values to mask.
-        head : numpy.ndarray
-            Three-dimensional array to set top of patches to the minimum
-            of the top of a layer or the head value. Used to create
-            patches that conform to water-level elevations.
-        **kwargs : dictionary
-            keyword arguments passed to matplotlib.pyplot.contour
-
-        Returns
-        -------
-        contour_set : matplotlib.pyplot.contour
-
-        """
-        if plt is None:
-            err_msg = (
-                "matplotlib must be installed to " + "use contour_array()"
-            )
-            raise ImportError(err_msg)
-        else:
-            import matplotlib.tri as tri
-
-        if not isinstance(a, np.ndarray):
-            a = np.array(a)
-
-        if a.ndim > 1:
-            a = np.ravel(a)
-
-        if "ax" in kwargs:
-            ax = kwargs.pop("ax")
-        else:
-            ax = self.ax
-
-        xcenters = [
-            np.mean(np.array(v).T[0]) for i, v in sorted(self.projpts.items())
-        ]
-
-        plotarray = np.array([a[cell] for cell in sorted(self.projpts)])
-
-        # work around for tri-contour ignore vmin & vmax
-        # necessary for the tri-contour NaN issue fix
-        if "levels" not in kwargs:
-            if "vmin" not in kwargs:
-                vmin = np.nanmin(plotarray)
-            else:
-                vmin = kwargs.pop("vmin")
-            if "vmax" not in kwargs:
-                vmax = np.nanmax(plotarray)
-            else:
-                vmax = kwargs.pop("vmax")
-
-            levels = np.linspace(vmin, vmax, 7)
-            kwargs["levels"] = levels
-
-        # workaround for tri-contour nan issue
-        plotarray[np.isnan(plotarray)] = -(2 ** 31)
-        if masked_values is None:
-            masked_values = [-(2 ** 31)]
-        else:
-            masked_values = list(masked_values)
-            if -(2 ** 31) not in masked_values:
-                masked_values.append(-(2 ** 31))
-
-        ismasked = None
-        if masked_values is not None:
-            for mval in masked_values:
-                if ismasked is None:
-                    ismasked = np.isclose(plotarray, mval)
-                else:
-                    t = np.isclose(plotarray, mval)
-                    ismasked += t
-
-        if isinstance(head, np.ndarray):
-            zcenters = self.set_zcentergrid(np.ravel(head))
-        else:
-            zcenters = [
-                np.mean(np.array(v).T[1])
-                for i, v in sorted(self.projpts.items())
-            ]
-
-        plot_triplot = False
-        if "plot_triplot" in kwargs:
-            plot_triplot = kwargs.pop("plot_triplot")
-
-        if "extent" in kwargs:
-            extent = kwargs.pop("extent")
-
-            idx = (
-                (xcenters >= extent[0])
-                & (xcenters <= extent[1])
-                & (zcenters >= extent[2])
-                & (zcenters <= extent[3])
-            )
-            plotarray = plotarray[idx].flatten()
-            xcenters = xcenters[idx].flatten()
-            zcenters = zcenters[idx].flatten()
-
-        triang = tri.Triangulation(xcenters, zcenters)
-
-        if ismasked is not None:
-            ismasked = ismasked.flatten()
-            mask = np.any(
-                np.where(ismasked[triang.triangles], True, False), axis=1
-            )
-            triang.set_mask(mask)
-
-        contour_set = ax.tricontour(triang, plotarray, **kwargs)
-
-        if plot_triplot:
-            ax.triplot(triang, color="black", marker="o", lw=0.75)
-
-        ax.set_xlim(self.extent[0], self.extent[1])
-        ax.set_ylim(self.extent[2], self.extent[3])
-
-        return contour_set
-
-    def plot_inactive(self):
-        raise NotImplementedError(
-            "Function must be called in PlotCrossSection"
-        )
-
-    def plot_ibound(self):
-        raise NotImplementedError(
-            "Function must be called in PlotCrossSection"
-        )
-
-    def plot_grid(self):
-        raise NotImplementedError(
-            "Function must be called in PlotCrossSection"
-        )
-
-    def plot_bc(self):
-        raise NotImplementedError(
-            "Function must be called in PlotCrossSection"
-        )
-
-    def plot_specific_discharge(self):
-        raise NotImplementedError(
-            "Function must be called in PlotCrossSection"
-        )
-
-    def plot_discharge(self):
-        raise NotImplementedError(
-            "plot_specific_discharge must be " "used for VertexGrid models"
-        )
-
-    @classmethod
-    def get_grid_patch_collection(cls, projpts, plotarray, **kwargs):
-        """
-        Get a PatchCollection of plotarray in unmasked cells
-
-        Parameters
-        ----------
-        projpts : dict
-            dictionary defined by node number which contains model patch vertices.
-        plotarray : numpy.ndarray
-            One-dimensional array to attach to the Patch Collection.
-        **kwargs : dictionary
-            keyword arguments passed to matplotlib.collections.PatchCollection
-
-        Returns
-        -------
-        patches : matplotlib.collections.PatchCollection
-
-        """
-        if plt is None:
-            err_msg = (
-                "matplotlib must be installed to "
-                + "use get_grid_patch_collection()"
-            )
-            raise ImportError(err_msg)
-        else:
-            from matplotlib.patches import Polygon
-            from matplotlib.collections import PatchCollection
-
-        if "vmin" in kwargs:
-            vmin = kwargs.pop("vmin")
-        else:
-            vmin = None
-        if "vmax" in kwargs:
-            vmax = kwargs.pop("vmax")
-        else:
-            vmax = None
-
-        rectcol = []
-        data = []
-        for cell, verts in sorted(projpts.items()):
-            verts = plotutil.UnstructuredPlotUtilities.arctan2(np.array(verts))
-
-            if np.isnan(plotarray[cell]):
-                continue
-            elif plotarray[cell] is np.ma.masked:
-                continue
-
-            rectcol.append(Polygon(verts, closed=True))
-            data.append(plotarray[cell])
-
-        if len(rectcol) > 0:
-            patches = PatchCollection(rectcol, **kwargs)
-            patches.set_array(np.array(data))
-            patches.set_clim(vmin, vmax)
-
-        else:
-            patches = None
-
-        return patches
-
-    def get_grid_line_collection(self, **kwargs):
-        """
-        Get a LineCollection of the grid
-
-        Parameters
-        ----------
-        **kwargs : dictionary
-            keyword arguments passed to matplotlib.collections.LineCollection
-
-        Returns
-        -------
-        linecollection : matplotlib.collections.LineCollection
-        """
-        if plt is None:
-            err_msg = (
-                "matplotlib must be installed to "
-                + "use get_grid_line_collection()"
-            )
-            raise ImportError(err_msg)
-        else:
-            from matplotlib.patches import Polygon
-            from matplotlib.collections import PatchCollection
-
-        color = "grey"
-        if "ec" in kwargs:
-            color = kwargs.pop("ec")
-        if color in kwargs:
-            color = kwargs.pop("color")
-
-        rectcol = []
-        for _, verts in sorted(self.projpts.items()):
-            verts = plotutil.UnstructuredPlotUtilities.arctan2(np.array(verts))
-
-            rectcol.append(Polygon(verts, closed=True))
-
-        if len(rectcol) > 0:
-            patches = PatchCollection(
-                rectcol, edgecolor=color, facecolor="none", **kwargs
-            )
-        else:
-            patches = None
-
-        return patches
-
-    def set_zpts(self, vs):
-        """
-        Get an array of projection vertices corrected for
-         elevations based on minimum of cell elevation
-        (self.elev) or passed vs numpy.ndarray
-
-        Parameters
-        ----------
-        vs : numpy.ndarray
-            Two-dimensional array to plot.
-
-        Returns
-        -------
-        zpts : dict
-
-        """
-        # make vertex array based on projection direction
-        if vs is not None:
-            if not isinstance(vs, np.ndarray):
-                vs = np.array(vs)
-
-        if self.direction == "x":
-            xyix = 0
-        else:
-            xyix = -1
-
-        projpts = {}
-        for k in range(1, self.mg.nlay + 1):
-            top = self.elev[k - 1, :]
-            botm = self.elev[k, :]
-            adjnn = (k - 1) * self.mg.ncpl
-            d0 = 0
-            for nn, verts in sorted(
-                self.xypts.items(), key=lambda q: q[-1][xyix][xyix]
-            ):
-                if vs is None:
-                    t = top[nn]
-                else:
-                    t = vs[nn]
-                    if top[nn] < vs[nn]:
-                        t = top[nn]
-                b = botm[nn]
-                if self.geographic_coords:
-                    if self.direction == "x":
-                        projt = [(v[0], t) for v in verts]
-                        projb = [(v[0], b) for v in verts]
-                    else:
-                        projt = [(v[1], t) for v in verts]
-                        projb = [(v[1], b) for v in verts]
-                else:
-                    verts = np.array(verts).T
-                    a2 = (np.max(verts[0]) - np.min(verts[0])) ** 2
-                    b2 = (np.max(verts[1]) - np.min(verts[1])) ** 2
-                    c = np.sqrt(a2 + b2)
-                    d1 = d0 + c
-                    projt = [(d0, t), (d1, t)]
-                    projb = [(d0, b), (d1, b)]
-                    d0 += c
-
-                projpts[nn + adjnn] = projt + projb
-
-        return projpts
-
-    def set_zcentergrid(self, vs, kstep=1):
-        """
-        Get an array of z elevations at the center of a cell that is based
-        on minimum of cell top elevation (self.elev) or passed vs numpy.ndarray
-
-        Parameters
-        ----------
-        vs : numpy.ndarray
-            Three-dimensional array to plot.
-
-        Returns
-        -------
-        zcentergrid : numpy.ndarray
-
-        """
-        verts = self.set_zpts(vs)
-        zcenters = [
-            np.mean(np.array(v).T[1])
-            for i, v in sorted(verts.items())
-            if (i // self.mg.ncpl) % kstep == 0
-        ]
-        return zcenters
-
-    def get_extent(self):
-        """
-        Get the extent of the rotated and offset grid
-
-        Returns
-        -------
-        tuple : (xmin, xmax, ymin, ymax)
-        """
-        xpts = []
-        for _, verts in self.projpts.items():
-            for v in verts:
-                xpts.append(v[0])
-
-        xmin = np.min(xpts)
-        xmax = np.max(xpts)
-
-        ymin = np.min(self.elev)
-        ymax = np.max(self.elev)
-
-        return (xmin, xmax, ymin, ymax)

flopy/utils/__init__.py

@@ -64,7 +64,7 @@ from .zonbud import (
     ZBNetOutput,
 )
 from .mfgrdfile import MfGrdFile
-from .postprocessing import get_transmissivities
+from .postprocessing import get_transmissivities, get_specific_discharge
 from .sfroutputfile import SfrFile
 from .recarray_utils import create_empty_recarray, ra_slice
 from .mtlistfile import MtListBudget

flopy/utils/binaryfile.py

@@ -300,6 +300,9 @@ class BinaryLayerFile(LayerFile):
         header = self._get_header()
         self.nrow = header["nrow"]
         self.ncol = header["ncol"]
+        if header["ilay"] > self.nlay:
+            self.nlay = header["ilay"]
+
         if self.nrow < 0 or self.ncol < 0:
             raise Exception("negative nrow, ncol")
         if self.nrow > 1 and self.nrow * self.ncol > 10000000:

flopy/utils/datafile.py

@@ -209,6 +209,7 @@ class LayerFile(object):
                 delr=np.ones(
                     self.ncol,
                 ),
+                nlay=self.nlay,
                 xoff=0.0,
                 yoff=0.0,
                 angrot=0.0,

flopy/utils/formattedfile.py

@@ -135,6 +135,8 @@ class FormattedLayerFile(LayerFile):
 
         self.nrow = header_info["nrow"]
         self.ncol = header_info["ncol"]
+        if header_info["ilay"] > self.nlay:
+            self.nlay = header_info["ilay"]
 
         ipos = self.file.tell()
         self._store_record(header_info, ipos)

flopy/utils/geometry.py

@@ -855,3 +855,112 @@ def is_clockwise(*geom):
         x = np.append(x, x[-1])
         y = np.append(y, y[-1])
     return np.sum((np.diff(x)) * (y[1:] + y[:-1])) > 0
+
+
+def point_in_polygon(xc, yc, polygon):
+    """
+    Use the ray casting algorithm to determine if a point
+    is within a polygon. Enables very fast
+    intersection calculations!
+
+    Parameters
+    ----------
+    xc : np.ndarray
+        2d array of xpoints
+    yc : np.ndarray
+        2d array of ypoints
+    polygon : iterable (list)
+        polygon vertices [(x0, y0),....(xn, yn)]
+        note: polygon can be open or closed
+
+    Returns
+    -------
+    mask: np.array
+        True value means point is in polygon!
+
+    """
+    x0, y0 = polygon[0]
+    xt, yt = polygon[-1]
+
+    # close polygon if it isn't already
+    if (x0, y0) != (xt, yt):
+        polygon.append((x0, y0))
+
+    ray_count = np.zeros(xc.shape, dtype=int)
+    num = len(polygon)
+    j = num - 1
+    for i in range(num):
+
+        tmp = polygon[i][0] + (polygon[j][0] - polygon[i][0]) * (
+            yc - polygon[i][1]
+        ) / (polygon[j][1] - polygon[i][1])
+
+        comp = np.where(
+            ((polygon[i][1] > yc) ^ (polygon[j][1] > yc)) & (xc < tmp)
+        )
+
+        j = i
+        if len(comp[0]) > 0:
+            ray_count[comp[0], comp[1]] += 1
+
+    mask = np.ones(xc.shape, dtype=bool)
+    mask[ray_count % 2 == 0] = False
+
+    return mask
+
+
+def project_point_onto_xc_line(line, pts, d0=0, direction="x"):
+    """
+    Method to project points onto a cross sectional line
+    that is defined by distance. Used for plotting MODPATH results
+    on to a cross section!
+
+    line : list or np.ndarray
+        numpy array of [(x0, y0), (x1, y1)] that defines the line
+        to project on to
+    pts : list or np.ndarray
+        numpy array of [(x, y),] points to be projected
+    d0 : distance offset along line of min(xl)
+    direction : string
+        projection direction "x" or "y"
+
+    Returns:
+        np.ndarray of projected [(x, y),] points
+    """
+    if isinstance(line, list):
+        line = np.array(line)
+
+    if isinstance(pts, list):
+        pts = np.array(pts)
+
+    x0, x1 = line.T[0, :]
+    y0, y1 = line.T[1, :]
+    dx = np.abs(x0 - x1)
+    dy = np.abs(y0 - y1)
+    m = dy / dx
+    b = y0 - (m * x0)
+    x = pts.T[0]
+    y = pts.T[1]
+
+    if direction == "x":
+        if dy == 0:
+            pass
+        else:
+            y = (x * m) + b
+
+    else:
+        if dx == 0:
+            pass
+        else:
+            x = (y - b) / m
+
+    # now do distance equation on pts from x0, y0
+    asq = (x - x0) ** 2
+    bsq = (y - y0) ** 2
+    dist = np.sqrt(asq + bsq)
+    if direction == "x":
+        x = dist + d0
+    else:
+        y = d0 - dist
+
+    return (x, y)

flopy/utils/postprocessing.py

@@ -177,8 +177,8 @@ def get_saturated_thickness(heads, m, nodata, per_idx=None):
         Heads array.
     m : flopy.modflow.Modflow object
         Must have a flopy.modflow.ModflowDis object attached.
-    nodata : real
-        HDRY value indicating dry cells.
+    nodata : float, list
+        HDRY value indicating dry cells and/or HNOFLO values.
     per_idx : int or sequence of ints
         stress periods to return. If None,
         returns all stress periods (default).
@@ -188,9 +188,17 @@ def get_saturated_thickness(heads, m, nodata, per_idx=None):
     sat_thickness : 3 or 4-D np.ndarray
         Array of saturated thickness
     """
-    # internal calculations done on a masked array
-    heads = np.ma.array(heads, ndmin=4, mask=heads == nodata)
+
+    if not isinstance(nodata, list):
+        nodata = [nodata]
+    heads = np.array(heads, ndmin=4)
+    for mv in nodata:
+        heads[heads == mv] = np.nan
+
+    top = m.dis.top.array
     botm = m.dis.botm.array
+    top.shape = (1,) + botm.shape[1:]
+    top = np.concatenate((top, botm[0:-1]), axis=0)
     thickness = m.dis.thickness.array
     nper, nlay, nrow, ncol = heads.shape
     if per_idx is None:
@@ -199,38 +207,25 @@ def get_saturated_thickness(heads, m, nodata, per_idx=None):
         per_idx = [per_idx]
 
     # get confined or unconfined/convertible info
-    if m.has_package("BCF6") or m.has_package("LPF") or m.has_package("UPW"):
-        if m.has_package("BCF6"):
-            laytyp = m.lpf.laycon.array
-        elif m.has_package("LPF"):
-            laytyp = m.lpf.laytyp.array
-        else:
-            laytyp = m.upw.laytyp.array
-        if len(laytyp) == 1:
-            is_conf = np.full(m.modelgrid.shape, laytyp == 0)
-        else:
-            laytyp = laytyp.reshape(m.modelgrid.nlay, 1, 1)
-            is_conf = np.logical_and(
-                (laytyp == 0), np.full(m.modelgrid.shape, True)
-            )
-    elif m.has_package("NPF"):
-        is_conf = m.npf.icelltype.array == 0
-    else:
-        raise ValueError(
-            "No flow package was found when trying to determine "
-            "the layer type."
+    laytyp = m.laytyp
+    if len(laytyp.shape) == 1:
+        laytyp.shape = (m.nlay, 1, 1)
+        is_conf = np.logical_and(
+            (laytyp == 0), np.full(m.modelgrid.shape, True)
         )
+    else:
+        is_conf = laytyp == 0
 
     # calculate saturated thickness
     sat_thickness = []
     for per in per_idx:
         hds = heads[per]
-        perthickness = hds - botm
-        conf = np.logical_or(perthickness > thickness, is_conf)
-        perthickness[conf] = thickness[conf]
-        # convert to nan-filled array, as is expected(!?)
-        sat_thickness.append(perthickness.filled(np.nan))
-    return np.squeeze(sat_thickness)
+        unconf_thickness = np.where((hds - botm) > top, top, hds - botm)
+        perthickness = np.where(is_conf, thickness, unconf_thickness)
+        sat_thickness.append(perthickness)
+    sat_thickness = np.squeeze(sat_thickness)
+
+    return sat_thickness
 
 
 def get_gradients(heads, m, nodata, per_idx=None):
@@ -344,9 +339,13 @@ def get_extended_budget(
         that the z axis is considered to increase in the upward direction.
     """
     import flopy.utils.binaryfile as bf
+    import flopy.utils.formattedfile as fm
 
     # define useful stuff
-    cbf = bf.CellBudgetFile(cbcfile, precision=precision)
+    if isinstance(cbcfile, bf.CellBudgetFile):
+        cbf = cbcfile
+    else:
+        cbf = bf.CellBudgetFile(cbcfile)
     nlay, nrow, ncol = cbf.nlay, cbf.nrow, cbf.ncol
     rec_names = cbf.get_unique_record_names(decode=True)
     err_msg = " not found in the budget file."
@@ -418,7 +417,13 @@ def get_extended_budget(
             raise ValueError(
                 "hdsfile must be provided when using " "boundary_ifaces"
             )
-        hds = bf.HeadFile(hdsfile, precision=precision)
+        if isinstance(hdsfile, (bf.HeadFile, fm.FormattedHeadFile)):
+            hds = hdsfile
+        else:
+            try:
+                hds = bf.HeadFile(hdsfile)
+            except:
+                hds = fm.FormattedHeadFile(hdsfile, precision=precision)
         head = hds.get_data(idx=idx, kstpkper=kstpkper, totim=totim)
 
         # get hnoflo and hdry values
@@ -599,14 +604,9 @@ def get_extended_budget(
 
 
 def get_specific_discharge(
+    vectors,
     model,
-    cbcfile,
-    precision="single",
-    idx=None,
-    kstpkper=None,
-    totim=None,
-    boundary_ifaces=None,
-    hdsfile=None,
+    head=None,
     position="centers",
 ):
     """
@@ -617,48 +617,16 @@ def get_specific_discharge(
 
     Parameters
     ----------
-    model : flopy.modflow.Modflow object
-        Modflow model instance.
-    cbcfile : str
-        Cell by cell file produced by Modflow.
-    precision : str
-        Binary file precision, default is 'single'.
-    idx : int or list
-            The zero-based record number.
-    kstpkper : tuple of ints
-        A tuple containing the time step and stress period (kstp, kper).
-        The kstp and kper values are zero based.
-    totim : float
-        The simulation time.
-    boundary_ifaces : dictionary {str: int or list}
-        A dictionary defining how to treat stress flows at boundary cells.
-        Only implemented for "classical" MODFLOW versions where the budget is
-        recorded as FLOW RIGHT FACE, FLOW FRONT FACE and FLOW LOWER FACE
-        arrays.
-        The keys are budget terms corresponding to stress packages (same term
-        as in the overall volumetric budget printed in the listing file).
-        The values are either a single iface number to be applied to all cells
-        for the stress package, or a list of lists describing individual
-        boundary cells in the same way as in the package input plus the iface
-        number appended. The iface number indicates the face to which the
-        stress flow is assigned, following the MODPATH convention (see MODPATH
-        user guide).
-        Example:
-        boundary_ifaces = {
-        'RECHARGE': 6,
-        'RIVER LEAKAGE': 6,
-        'WELLS': [[lay, row, col, flux, iface], ...],
-        'HEAD DEP BOUNDS': [[lay, row, col, head, cond, iface], ...]}.
-        Note: stresses that are not informed in boundary_ifaces are implicitly
-        treated as internally-distributed sinks/sources.
-    hdsfile : str
-        Head file produced by MODFLOW. Head is used to calculate saturated
-        thickness and to determine if a cell is inactive or dry. If not
-        provided, all cells are considered fully saturated.
-        hdsfile is also required if the budget term 'HEAD DEP BOUNDS',
-        'RIVER LEAKAGE' or 'DRAINS' is present in boundary_ifaces and that the
-        corresponding value is a list.
-    position : str
+    vectors : tuple, np.recarray
+        either a tuple of (flow right face, flow front face, flow lower face)
+        numpy arrays from a MODFLOW-2005 compatible Cell Budget File
+        or
+        a specific discharge recarray from a MODFLOW 6 Cell Budget File
+    model : object
+        flopy model object
+    head : np.ndarray
+        numpy array of head values for a specific model
+     position : str
         Position at which the specific discharge will be calculated. Possible
         values are "centers" (default), "faces" and "vertices".
 
@@ -672,54 +640,79 @@ def get_specific_discharge(
         in the north direction.
         The sign of qz is such that the z axis is considered to increase
         in the upward direction.
-        Note: if hdsfile is provided, inactive and dry cells are set to NaN.
+        Note: if a head array is provided, inactive and dry cells are
+        set to NaN.
     """
-    import flopy.utils.binaryfile as bf
 
-    # check if budget file has classical budget terms
-    cbf = bf.CellBudgetFile(cbcfile, precision=precision)
-    rec_names = cbf.get_unique_record_names(decode=True)
-    classical_budget_terms = [
-        "FLOW RIGHT FACE",
-        "FLOW FRONT FACE",
-        "FLOW RIGHT FACE",
-    ]
     classical_budget = False
-    for budget_term in classical_budget_terms:
-        matched_name = [s for s in rec_names if budget_term in s]
-        if matched_name:
-            classical_budget = True
-            break
+    spdis, tqx, tqy, tqz = None, None, None, None
+    modelgrid = model.modelgrid
 
-    if hdsfile is not None:
-        hds = bf.HeadFile(hdsfile, precision=precision)
-        head = hds.get_data(idx=idx, kstpkper=kstpkper, totim=totim)
+    if head is not None:
+        head.shape = modelgrid.shape
+
+    if isinstance(vectors, (list, tuple)):
+        classical_budget = True
+        for ix, vector in enumerate(vectors):
+            if vector is None:
+                continue
+            else:
+                tshp = list(modelgrid.shape)[::-1]
+                tshp[ix] += 1
+                ext_shape = tuple(tshp[::-1])
+                break
+
+        if vectors[ix].shape == modelgrid.shape:
+            tqx = np.zeros(
+                (modelgrid.nlay, modelgrid.nrow, modelgrid.ncol + 1),
+                dtype=np.float32,
+            )
+            tqy = np.zeros(
+                (modelgrid.nlay, modelgrid.nrow + 1, modelgrid.ncol),
+                dtype=np.float32,
+            )
+            tqz = np.zeros(
+                (modelgrid.nlay + 1, modelgrid.nrow, modelgrid.ncol),
+                dtype=np.float32,
+            )
+            if vectors[0] is not None:
+                tqx[:, :, 1:] = vectors[0]
+            if vectors[1] is not None:
+                tqy[:, 1:, :] = -vectors[1]
+            if len(vectors) > 2 and vectors[2] is not None:
+                tqz[1:, :, :] = -vectors[2]
+
+        elif vectors[ix].shape == ext_shape:
+            if vectors[0] is not None:
+                tqx = vectors[0]
+            if vectors[1] is not None:
+                tqy = vectors[1]
+            if len(vectors) > 2 and vectors[2] is not None:
+                tqz = vectors[2]
+
+        else:
+            raise IndexError(
+                "Classical budget components must have "
+                "the same shape as the modelgrid"
+            )
+    else:
+        spdis = vectors
 
     if classical_budget:
-        # get extended budget
-        Qx_ext, Qy_ext, Qz_ext = get_extended_budget(
-            cbcfile,
-            precision=precision,
-            idx=idx,
-            kstpkper=kstpkper,
-            totim=totim,
-            boundary_ifaces=boundary_ifaces,
-            hdsfile=hdsfile,
-            model=model,
-        )
-
         # get saturated thickness (head - bottom elev for unconfined layer)
-        if hdsfile is None:
+        if head is None:
             sat_thk = model.dis.thickness.array
         else:
-            sat_thk = get_saturated_thickness(head, model, model.hdry)
-            sat_thk = sat_thk.reshape(model.modelgrid.shape)
+            sat_thk = get_saturated_thickness(
+                head, model, [model.hdry, model.hnoflo]
+            )
+            sat_thk.shape = model.modelgrid.shape
 
         # inform modelgrid of no-flow and dry cells
         modelgrid = model.modelgrid
         if modelgrid._idomain is None:
             modelgrid._idomain = model.dis.ibound
-        if hdsfile is not None:
+        if head is not None:
             noflo_or_dry = np.logical_or(
                 head == model.hnoflo, head == model.hdry
             )
@@ -727,12 +720,10 @@ def get_specific_discharge(
 
         # get cross section areas along x
         delc = np.reshape(modelgrid.delc, (1, modelgrid.nrow, 1))
-        cross_area_x = np.empty(modelgrid.shape, dtype=float)
         cross_area_x = delc * sat_thk
 
         # get cross section areas along y
         delr = np.reshape(modelgrid.delr, (1, 1, modelgrid.ncol))
-        cross_area_y = np.empty(modelgrid.shape, dtype=float)
         cross_area_y = delr * sat_thk
 
         # get cross section areas along z
@@ -740,16 +731,31 @@ def get_specific_discharge(
 
         # calculate qx, qy, qz
         if position == "centers":
-            qx = 0.5 * (Qx_ext[:, :, 1:] + Qx_ext[:, :, :-1]) / cross_area_x
-            qy = 0.5 * (Qy_ext[:, 1:, :] + Qy_ext[:, :-1, :]) / cross_area_y
-            qz = 0.5 * (Qz_ext[1:, :, :] + Qz_ext[:-1, :, :]) / cross_area_z
+            qx = np.zeros(modelgrid.shape, dtype=np.float32)
+            qy = np.zeros(modelgrid.shape, dtype=np.float32)
+            cross_area_x = (
+                delc[:] * 0.5 * (sat_thk[:, :, :-1] + sat_thk[:, :, 1:])
+            )
+            cross_area_y = (
+                delr * 0.5 * (sat_thk[:, 1:, :] + sat_thk[:, :-1, :])
+            )
+            qx[:, :, 1:] = (
+                0.5 * (tqx[:, :, 2:] + tqx[:, :, 1:-1]) / cross_area_x
+            )
+            qx[:, :, 0] = 0.5 * tqx[:, :, 1] / cross_area_x[:, :, 0]
+            qy[:, 1:, :] = (
+                0.5 * (tqy[:, 2:, :] + tqy[:, 1:-1, :]) / cross_area_y
+            )
+            qy[:, 0, :] = 0.5 * tqy[:, 1, :] / cross_area_y[:, 0, :]
+            qz = 0.5 * (tqz[1:, :, :] + tqz[:-1, :, :]) / cross_area_z
+
         elif position == "faces" or position == "vertices":
             cross_area_x = modelgrid.array_at_faces(cross_area_x, "x")
             cross_area_y = modelgrid.array_at_faces(cross_area_y, "y")
             cross_area_z = modelgrid.array_at_faces(cross_area_z, "z")
-            qx = Qx_ext / cross_area_x
-            qy = Qy_ext / cross_area_y
-            qz = Qz_ext / cross_area_z
+            qx = tqx / cross_area_x
+            qy = tqy / cross_area_y
+            qz = tqz / cross_area_z
         else:
             raise ValueError(
                 '"' + position + '" is not a valid value for ' "position"
@@ -760,52 +766,20 @@ def get_specific_discharge(
             qz = modelgrid.array_at_verts(qz)
 
     else:
-        # check valid options
-        if boundary_ifaces is not None:
-            import warnings
-
-            warnings.warn(
-                "the boundary_ifaces option is not implemented "
-                'for "non-classical" MODFLOW versions where the '
-                "budget is not recorded as FLOW RIGHT FACE, "
-                "FLOW FRONT FACE and FLOW LOWER FACE; it will be "
-                "ignored",
-                UserWarning,
-            )
-        if position != "centers":
-            raise NotImplementedError(
-                'position can only be "centers" for '
-                '"non-classical" MODFLOW versions where '
-                "the budget is not recorded as FLOW "
-                "RIGHT FACE, FLOW FRONT FACE and FLOW "
-                "LOWER FACE"
-            )
-
-        is_spdis = [s for s in rec_names if "DATA-SPDIS" in s]
-        if not is_spdis:
-            err_msg = (
-                "Could not find suitable records in the budget file "
-                "to construct the discharge vector."
-            )
-            raise RuntimeError(err_msg)
-        spdis = cbf.get_data(
-            text="DATA-SPDIS", idx=idx, kstpkper=kstpkper, totim=totim
-        )[0]
         nnodes = model.modelgrid.nnodes
-        qx = np.full((nnodes), np.nan)
-        qy = np.full((nnodes), np.nan)
-        qz = np.full((nnodes), np.nan)
+        qx = np.full((nnodes), np.nan, dtype=np.float64)
+        qy = np.full((nnodes), np.nan, dtype=np.float64)
+        qz = np.full((nnodes), np.nan, dtype=np.float64)
         idx = np.array(spdis["node"]) - 1
         qx[idx] = spdis["qx"]
         qy[idx] = spdis["qy"]
         qz[idx] = spdis["qz"]
-        shape = model.modelgrid.shape
-        qx.shape = shape
-        qy.shape = shape
-        qz.shape = shape
+        qx.shape = modelgrid.shape
+        qy.shape = modelgrid.shape
+        qz.shape = modelgrid.shape
 
     # set no-flow and dry cells to NaN
-    if hdsfile is not None and position == "centers":
+    if head is not None and position == "centers":
         noflo_or_dry = np.logical_or(head == model.hnoflo, head == model.hdry)
         qx[noflo_or_dry] = np.nan
         qy[noflo_or_dry] = np.nan

flopy/utils/rasters.py

@@ -1,4 +1,6 @@
 import numpy as np
+import threading
+import queue
 
 try:
     import rasterio
@@ -55,7 +57,7 @@ class Raster(object):
     FLOAT64 = (np.float64,)
     INT8 = (np.int8,)
     INT16 = (np.int16,)
-    INT32 = (int, np.int32, np.int_)
+    INT32 = (int, np.int, np.int32, np.uint32)
     INT64 = (np.int64,)
 
     def __init__(
@@ -84,6 +86,9 @@ class Raster(object):
             )
             raise ImportError(msg)
 
+        from .geometry import point_in_polygon
+
+        self._point_in_polygon = point_in_polygon
         self._array = array
         self._bands = bands
 
@@ -336,7 +341,14 @@ class Raster(object):
 
         return arr_dict[band]
 
-    def resample_to_grid(self, xc, yc, band, method="nearest"):
+    def resample_to_grid(
+        self,
+        modelgrid,
+        band,
+        method="nearest",
+        multithread=False,
+        thread_pool=2,
+    ):
         """
         Method to resample the raster data to a
         user supplied grid of x, y coordinates.
@@ -346,10 +358,8 @@ class Raster(object):
 
         Parameters
         ----------
-        xc : np.ndarray or list
-            an array of x-cell centers
-        yc : np.ndarray or list
-            an array of y-cell centers
+        modelgrid : flopy.Grid object
+            model grid to sample data from
         band : int
             raster band to re-sample
         method : str
@@ -371,26 +381,75 @@ class Raster(object):
         else:
             from scipy.interpolate import griddata
 
-        data_shape = xc.shape
-        xc = xc.flatten()
-        yc = yc.flatten()
-        # step 1: create grid from raster bounds
-        rxc = self.xcenters
-        ryc = self.ycenters
+        method = method.lower()
+        if method in ("linear", "nearest", "cubic"):
+            xc = modelgrid.xcellcenters
+            yc = modelgrid.ycellcenters
 
-        # step 2: flatten grid
-        rxc = rxc.flatten()
-        ryc = ryc.flatten()
+            data_shape = xc.shape
+            xc = xc.flatten()
+            yc = yc.flatten()
+            # step 1: create grid from raster bounds
+            rxc = self.xcenters
+            ryc = self.ycenters
 
-        # step 3: get array
-        if method == "cubic":
-            arr = self.get_array(band, masked=False)
+            # step 2: flatten grid
+            rxc = rxc.flatten()
+            ryc = ryc.flatten()
+
+            # step 3: get array
+            if method == "cubic":
+                arr = self.get_array(band, masked=False)
+            else:
+                arr = self.get_array(band, masked=True)
+            arr = arr.flatten()
+
+            # step 3: use griddata interpolation to snap to grid
+            data = griddata((rxc, ryc), arr, (xc, yc), method=method)
+
+        elif method in ("median", "mean"):
+            # these methods are slow and could use a speed u
+            ncpl = modelgrid.ncpl
+            data_shape = modelgrid.xcellcenters.shape
+            if isinstance(ncpl, (list, np.ndarray)):
+                ncpl = ncpl[0]
+
+            data = np.zeros((ncpl,), dtype=float)
+
+            if multithread:
+                q = queue.Queue()
+                container = threading.BoundedSemaphore(thread_pool)
+                threads = []
+                for node in range(ncpl):
+                    t = threading.Thread(
+                        target=self.__threaded_resampling,
+                        args=(modelgrid, node, band, method, container, q),
+                    )
+                    threads.append(t)
+
+                for thread in threads:
+                    thread.start()
+                for thread in threads:
+                    thread.join()
+
+                for _ in range(len(threads)):
+                    node, val = q.get()
+                    data[node] = val
+            else:
+                for node in range(ncpl):
+                    verts = modelgrid.get_cell_vertices(node)
+                    rstr_data = self.sample_polygon(verts, band)
+                    msk = np.in1d(rstr_data, self.nodatavals)
+                    rstr_data[msk] = np.nan
+
+                    if method == "median":
+                        val = np.nanmedian(rstr_data)
+                    else:
+                        val = np.nanmean(rstr_data)
+
+                    data[node] = val
         else:
-            arr = self.get_array(band, masked=True)
-        arr = arr.flatten()
-
-        # step 3: use griddata interpolation to snap to grid
-        data = griddata((rxc, ryc), arr, (xc, yc), method=method)
+            raise TypeError("{} method not supported".format(method))
 
         # step 4: return grid to user in shape provided
         data.shape = data_shape
@@ -400,6 +459,42 @@ class Raster(object):
 
         return data
 
+    def __threaded_resampling(
+        self, modelgrid, node, band, method, container, q
+    ):
+        """
+        Threaded resampling handler to speed up bottlenecks
+
+        Parameters
+        ----------
+        modelgrid : flopy.discretization.Grid object
+            flopy grid to sample to
+        node : int
+            node number
+        band : int
+            raster band to sample from
+        method : str
+            resampling method
+        container : threading.BoundedSemaphore
+        q : queue.Queue
+
+        Returns
+        -------
+            None
+        """
+        container.acquire()
+        verts = modelgrid.get_cell_vertices(node)
+        rstr_data = self.sample_polygon(verts, band)
+        msk = np.in1d(rstr_data, self.nodatavals)
+        rstr_data[msk] = np.nan
+
+        if method == "median":
+            val = np.nanmedian(rstr_data)
+        else:
+            val = np.nanmean(rstr_data)
+        q.put((node, val))
+        container.release()
+
     def crop(self, polygon, invert=False):
         """
         Method to crop a new raster object
@@ -625,58 +720,6 @@ class Raster(object):
 
         return mask
 
-    @staticmethod
-    def _point_in_polygon(xc, yc, polygon):
-        """
-        Use the ray casting algorithm to determine if a point
-        is within a polygon. Enables very fast
-        intersection calculations!
-
-        Parameters
-        ----------
-        xc : np.ndarray
-            array of xpoints
-        yc : np.ndarray
-            array of ypoints
-        polygon : iterable (list)
-            polygon vertices [(x0, y0),....(xn, yn)]
-            note: polygon can be open or closed
-
-        Returns
-        -------
-        mask: np.array
-            True value means point is in polygon!
-
-        """
-        x0, y0 = polygon[0]
-        xt, yt = polygon[-1]
-
-        # close polygon if it isn't already
-        if (x0, y0) != (xt, yt):
-            polygon.append((x0, y0))
-
-        ray_count = np.zeros(xc.shape, dtype=int)
-        num = len(polygon)
-        j = num - 1
-        for i in range(num):
-
-            tmp = polygon[i][0] + (polygon[j][0] - polygon[i][0]) * (
-                yc - polygon[i][1]
-            ) / (polygon[j][1] - polygon[i][1])
-
-            comp = np.where(
-                ((polygon[i][1] > yc) ^ (polygon[j][1] > yc)) & (xc < tmp)
-            )
-
-            j = i
-            if len(comp[0]) > 0:
-                ray_count[comp[0], comp[1]] += 1
-
-        mask = np.ones(xc.shape, dtype=bool)
-        mask[ray_count % 2 == 0] = False
-
-        return mask
-
     def get_array(self, band, masked=True):
         """
         Method to get a numpy array corresponding to the
@@ -706,6 +749,7 @@ class Raster(object):
 
         if masked:
             for v in self.nodatavals:
+                array = array.astype(float)
                 array[array == v] = np.nan
 
         return array

flopy/utils/reference.py

@@ -138,8 +138,8 @@ class SpatialReference(object):
         length_multiplier=None,
     ):
         warnings.warn(
-            "SpatialReference has been deprecated. Use StructuredGrid"
-            " instead.",
+            "SpatialReference has been deprecated and will be removed in "
+            "version 3.3.5. Use StructuredGrid instead.",
             category=DeprecationWarning,
         )
 
@@ -923,7 +923,7 @@ class SpatialReference(object):
         Get a LineCollection of the grid
 
         """
-        from flopy.plot import ModelMap
+        from ..plot import ModelMap
 
         map = ModelMap(sr=self)
         lc = map.plot_grid(**kwargs)
@@ -1708,8 +1708,8 @@ class SpatialReferenceUnstructured(SpatialReference):
         length_multiplier=1.0,
     ):
         warnings.warn(
-            "SpatialReferenceUnstructured has been deprecated. "
-            "Use VertexGrid instead.",
+            "SpatialReferenceUnstructured has been deprecated and will be "
+            "removed in version 3.3.5. Use VertexGrid instead.",
             category=DeprecationWarning,
         )
         self.xc = xc
@@ -1905,27 +1905,26 @@ class SpatialReferenceUnstructured(SpatialReference):
 
         Returns
         -------
-        quadmesh : matplotlib.collections.QuadMesh
+        pc : matplotlib.collections.PatchCollection
 
         """
-        from ..plot import plotutil
+        from ..plot import ModelMap
 
-        patch_collection = plotutil.plot_cvfd(
-            self.verts, self.iverts, a=a, ax=ax
-        )
-        return patch_collection
+        pmv = ModelMap(sr=self, ax=ax)
+        pc = pmv.plot_array(a)
+
+        return pc
 
     def get_grid_line_collection(self, **kwargs):
         """
         Get a patch collection of the grid
 
         """
-        from ..plot import plotutil
+        from ..plot import ModelMap
 
-        edgecolor = kwargs.pop("colors")
-        pc = plotutil.cvfd_to_patch_collection(self.verts, self.iverts)
-        pc.set(facecolor="none")
-        pc.set(edgecolor=edgecolor)
+        ax = kwargs.pop("ax", None)
+        pmv = ModelMap(sr=self, ax=ax)
+        pc = pmv.plot_grid(**kwargs)
         return pc
 
     def contour_array(self, ax, a, **kwargs):
@@ -1989,7 +1988,9 @@ class epsgRef:
 
     def __init__(self):
         warnings.warn(
-            "epsgRef has been deprecated.", category=DeprecationWarning
+            "epsgRef has been deprecated and will be removed in version "
+            "3.3.5.",
+            category=DeprecationWarning,
         )
         try:
             from appdirs import user_data_dir
@@ -2072,7 +2073,8 @@ class crs(object):
 
     def __init__(self, prj=None, esri_wkt=None, epsg=None):
         warnings.warn(
-            "crs has been deprecated. Use CRS in shapefile_utils instead.",
+            "crs has been deprecated and will be removed in version 3.3.5. "
+            "Use CRS in shapefile_utils instead.",
             category=DeprecationWarning,
         )
         self.wktstr = None
@@ -2293,7 +2295,8 @@ def getprj(epsg, addlocalreference=True, text="esriwkt"):
 
     """
     warnings.warn(
-        "SpatialReference has been deprecated. Use StructuredGrid " "instead.",
+        "SpatialReference has been deprecated and will be removed in version "
+        "3.3.5. Use StructuredGrid instead.",
         category=DeprecationWarning,
     )
     epsgfile = epsgRef()
@@ -2329,7 +2332,8 @@ def get_spatialreference(epsg, text="esriwkt"):
     from flopy.utils.flopy_io import get_url_text
 
     warnings.warn(
-        "SpatialReference has been deprecated. Use StructuredGrid " "instead.",
+        "SpatialReference has been deprecated and will be removed in version "
+        "3.3.5. Use StructuredGrid instead.",
         category=DeprecationWarning,
     )
 
@@ -2373,7 +2377,8 @@ def getproj4(epsg):
 
     """
     warnings.warn(
-        "SpatialReference has been deprecated. Use StructuredGrid " "instead.",
+        "SpatialReference has been deprecated and will be removed in version "
+        "3.3.5. Use StructuredGrid instead.",
         category=DeprecationWarning,
     )
 

flopy/utils/triangle.py

@@ -3,7 +3,6 @@ import numpy as np
 import subprocess
 from ..mbase import which
 from ..utils.cvfdutil import centroid_of_polygon
-from ..plot.plotutil import plot_cvfd
 from ..utils.geospatial_utils import GeoSpatialUtil
 
 
@@ -241,29 +240,31 @@ class Triangle(object):
         None
 
         """
-        try:
-            import matplotlib.pyplot as plt
-        except:
-            err_msg = (
-                "matplotlib must be installed to " + "use triangle.plot()"
-            )
-            raise ImportError(err_msg)
+        from ..plot import PlotMapView
+        from ..discretization import VertexGrid
 
-        if ax is None:
-            ax = plt.gca()
+        cell2d = self.get_cell2d()
+        vertices = self.get_vertices()
+        ncpl = len(cell2d)
 
-        pc = plot_cvfd(
-            self.verts,
-            self.iverts,
-            ax=ax,
-            edgecolor=edgecolor,
-            facecolor=facecolor,
-            cmap=cmap,
-            a=a,
-            masked_values=masked_values,
-            **kwargs
+        modelgrid = VertexGrid(
+            vertices=vertices, cell2d=cell2d, ncpl=ncpl, nlay=1
         )
-        ax.autoscale()
+
+        pmv = PlotMapView(modelgrid=modelgrid, ax=ax, layer=layer)
+        if a is None:
+            pc = pmv.plot_grid(
+                facecolor=facecolor, edgecolor=edgecolor, **kwargs
+            )
+        else:
+            pc = pmv.plot_array(
+                a,
+                masked_values=masked_values,
+                cmap=cmap,
+                edgecolor=edgecolor,
+                **kwargs
+            )
+
         return pc
 
     def get_boundary_marker_array(self):

flopy/utils/voronoi.py

@@ -2,7 +2,6 @@ import numpy as np
 import matplotlib.pyplot as plt
 from scipy.spatial import Voronoi
 from .cvfdutil import get_disv_gridprops
-from ..plot.plotutil import plot_cvfd
 
 
 def get_sorted_vertices(icell_vertices, vertices):
@@ -28,6 +27,7 @@ def get_valid_faces(vor):
     return nvalid_faces
 
 
+# todo: send this to point in polygon method defined in Rasters
 def point_in_cell(point, vertices):
     from shapely.geometry import Point, Polygon
 
@@ -39,6 +39,7 @@ def point_in_cell(point, vertices):
         return False
 
 
+# todo: find out how this is different from get_sorted_vertices()
 def sort_vertices(vlist):
     x, y = zip(*vlist)
     x = np.array(x)
@@ -251,7 +252,12 @@ class VoronoiGrid:
             patch collection of model
 
         """
-        pc = plot_cvfd(self.verts, self.iverts, ax=ax, **kwargs)
+        from ..discretization import VertexGrid
+        from ..plot import PlotMapView
+
+        modelgrid = VertexGrid(**self.get_gridprops_vertexgrid())
+        pmv = PlotMapView(modelgrid=modelgrid, ax=ax)
+        pc = pmv.plot_grid(**kwargs)
         return pc
 
     def plot(self, ax=None, plot_title=True, **kwargs):
@@ -276,8 +282,6 @@ class VoronoiGrid:
         if ax is None:
             ax = plt.subplot(1, 1, 1, aspect="equal")
         pc = self.get_patch_collection(ax, **kwargs)
-        ax.set_xlim(self.verts[:, 0].min(), self.verts[:, 0].max())
-        ax.set_ylim(self.verts[:, 1].min(), self.verts[:, 1].max())
         if plot_title:
             ax.set_title(
                 "ncells: {}; nverts: {}".format(self.ncpl, self.nverts)