feat(plot_bc): updated plot_bc for MAW, UZF, SFR, and multiple bc pac…

…kages (#808)

* fixed bug with mf2k loading #802
* added shapefile support for output_helper #807
* updated test cases

autotest/t007_test.py

@@ -163,6 +163,19 @@ def export_shapefile_modelgrid_override(namfile):
         raise Exception(msg)
 
 
+def test_output_helper_shapefile_export():
+    ws = os.path.join('..', 'examples', 'data', 'freyberg_multilayer_transient')
+    name = 'freyberg.nam'
+
+    ml = flopy.modflow.Modflow.load(name, model_ws=ws)
+
+    head = flopy.utils.HeadFile(os.path.join(ws, 'freyberg.hds'))
+    cbc = flopy.utils.CellBudgetFile(os.path.join(ws, "freyberg.cbc"))
+    flopy.export.utils.output_helper(os.path.join('temp', 'test.shp'), ml,
+                                     {'HDS': head, 'cbc': cbc},
+                                     mflay=1, kper=10)
+
+
 def test_freyberg_export():
     from flopy.discretization import StructuredGrid
     namfile = 'freyberg.nam'
@@ -1119,6 +1132,174 @@ def check_vertices():
     plt.close()
 
 
+def test_mapview_plot_bc():
+    from matplotlib.collections import QuadMesh, PatchCollection
+    import matplotlib.pyplot as plt
+
+
+    sim_name = 'mfsim.nam'
+    sim_path = os.path.join("..", "examples", "data", "mf6",
+                            "test003_gwfs_disv")
+    sim = flopy.mf6.MFSimulation.load(sim_name=sim_name,
+                                      sim_ws=sim_path)
+    ml6 = sim.get_model("gwf_1")
+    ml6.modelgrid.set_coord_info(angrot=-14)
+    mapview = flopy.plot.PlotMapView(model=ml6)
+    mapview.plot_bc('CHD')
+    ax = mapview.ax
+
+    if len(ax.collections) == 0:
+        raise AssertionError("Boundary condition was not drawn")
+
+    for col in ax.collections:
+        if not isinstance(col, PatchCollection):
+            raise AssertionError("Unexpected collection type")
+    plt.close()
+
+    sim_name = 'mfsim.nam'
+    sim_path = os.path.join('..', 'examples', 'data', 'mf6', 'test045_lake2tr')
+    sim = flopy.mf6.MFSimulation.load(sim_name=sim_name,
+                                      sim_ws=sim_path)
+
+    ml6 = sim.get_model("lakeex2a")
+    mapview = flopy.plot.PlotMapView(model=ml6)
+    mapview.plot_bc('LAK')
+    mapview.plot_bc("SFR")
+
+    ax = mapview.ax
+    if len(ax.collections) == 0:
+        raise AssertionError("Boundary condition was not drawn")
+
+    for col in ax.collections:
+        if not isinstance(col, QuadMesh):
+            raise AssertionError("Unexpected collection type")
+    plt.close()
+
+    sim_name = 'mfsim.nam'
+    sim_path = os.path.join('..', 'examples', 'data', 'mf6',
+                            'test006_2models_mvr')
+    sim = flopy.mf6.MFSimulation.load(sim_name=sim_name,
+                                      sim_ws=sim_path)
+
+    ml6 = sim.get_model("parent")
+    ml6c = sim.get_model('child')
+    ml6c.modelgrid.set_coord_info(xoff=700, yoff=0, angrot=0)
+
+    mapview = flopy.plot.PlotMapView(model=ml6)
+    mapview.plot_bc("MAW")
+
+    mapview2 = flopy.plot.PlotMapView(model=ml6c, ax=mapview.ax)
+    mapview2.plot_bc("MAW")
+    ax = mapview2.ax
+
+    if len(ax.collections) == 0:
+        raise AssertionError("Boundary condition was not drawn")
+
+    for col in ax.collections:
+        if not isinstance(col, QuadMesh):
+            raise AssertionError("Unexpected collection type")
+    plt.close()
+
+    sim_name = 'mfsim.nam'
+    sim_path = os.path.join('..', 'examples', 'data', 'mf6',
+                            'test001e_UZF_3lay')
+    sim = flopy.mf6.MFSimulation.load(sim_name=sim_name,
+                                      sim_ws=sim_path)
+    ml6 = sim.get_model("gwf_1")
+
+    mapview = flopy.plot.PlotMapView(model=ml6)
+    mapview.plot_bc("UZF")
+
+    if len(ax.collections) == 0:
+        raise AssertionError("Boundary condition was not drawn")
+
+    for col in ax.collections:
+        if not isinstance(col, QuadMesh):
+            raise AssertionError("Unexpected collection type")
+    plt.close()
+
+
+def test_crosssection_plot_bc():
+    from matplotlib.collections import PatchCollection
+    import matplotlib.pyplot as plt
+
+    sim_name = 'mfsim.nam'
+    sim_path = os.path.join("..", "examples", "data", "mf6",
+                            "test003_gwfs_disv")
+    sim = flopy.mf6.MFSimulation.load(sim_name=sim_name,
+                                      sim_ws=sim_path)
+    ml6 = sim.get_model("gwf_1")
+    xc = flopy.plot.PlotCrossSection(ml6, line={'line': ([0, 5.5],
+                                                         [10, 5.5])})
+    xc.plot_bc('CHD')
+    ax = xc.ax
+
+    if len(ax.collections) == 0:
+        raise AssertionError("Boundary condition was not drawn")
+
+    for col in ax.collections:
+        if not isinstance(col, PatchCollection):
+            raise AssertionError("Unexpected collection type")
+    plt.close()
+
+    sim_name = 'mfsim.nam'
+    sim_path = os.path.join('..', 'examples', 'data', 'mf6', 'test045_lake2tr')
+    sim = flopy.mf6.MFSimulation.load(sim_name=sim_name,
+                                      sim_ws=sim_path)
+
+    ml6 = sim.get_model("lakeex2a")
+    xc = flopy.plot.PlotCrossSection(ml6, line={'row': 10})
+    xc.plot_bc('LAK')
+    xc.plot_bc("SFR")
+
+    ax = xc.ax
+    if len(ax.collections) == 0:
+        raise AssertionError("Boundary condition was not drawn")
+
+    for col in ax.collections:
+        if not isinstance(col, PatchCollection):
+            raise AssertionError("Unexpected collection type")
+    plt.close()
+
+    sim_name = 'mfsim.nam'
+    sim_path = os.path.join('..', 'examples', 'data', 'mf6',
+                            'test006_2models_mvr')
+    sim = flopy.mf6.MFSimulation.load(sim_name=sim_name,
+                                      sim_ws=sim_path)
+
+    ml6 = sim.get_model("parent")
+    xc = flopy.plot.PlotCrossSection(ml6, line={'column': 1})
+    xc.plot_bc("MAW")
+
+    ax = xc.ax
+    if len(ax.collections) == 0:
+        raise AssertionError("Boundary condition was not drawn")
+
+    for col in ax.collections:
+        if not isinstance(col, PatchCollection):
+            raise AssertionError("Unexpected collection type")
+    plt.close()
+
+    sim_name = 'mfsim.nam'
+    sim_path = os.path.join('..', 'examples', 'data', 'mf6',
+                            'test001e_UZF_3lay')
+    sim = flopy.mf6.MFSimulation.load(sim_name=sim_name,
+                                      sim_ws=sim_path)
+    ml6 = sim.get_model("gwf_1")
+
+    xc = flopy.plot.PlotCrossSection(ml6, line={"row": 0})
+    xc.plot_bc("UZF")
+
+    ax = xc.ax
+    if len(ax.collections) == 0:
+        raise AssertionError("Boundary condition was not drawn")
+
+    for col in ax.collections:
+        if not isinstance(col, PatchCollection):
+            raise AssertionError("Unexpected collection type")
+    plt.close()
+
+
 def test_tricontour_NaN():
     from flopy.plot import PlotMapView
     import numpy as np
@@ -1258,43 +1439,6 @@ def test_netcdf_classmethods():
     assert len(diff) == 0, str(diff)
 
 
-# def test_netcdf_overloads():
-#     import os
-#     import flopy
-#     nam_file = "freyberg.nam"
-#     model_ws = os.path.join('..', 'examples', 'data', 'freyberg_multilayer_transient')
-#     ml = flopy.modflow.Modflow.load(nam_file,model_ws=model_ws,check=False,
-#                                     verbose=False,load_only=[])
-#
-#     f = ml.export(os.path.join("temp","freyberg.nc"))
-#     fzero = flopy.export.NetCdf.zeros_like(f)
-#     assert fzero.nc.variables["model_top"][:].sum() == 0
-#     print(f.nc.variables["model_top"][0,:])
-#     fplus1 = f + 1
-#     assert fplus1.nc.variables["model_top"][0,0] == f.nc.variables["model_top"][0,0] + 1
-#     assert (f + fplus1).nc.variables["model_top"][0,0] ==\
-#            f.nc.variables["model_top"][0,0] + \
-#            fplus1.nc.variables["model_top"][0,0]
-#
-#     fminus1 = f - 1
-#     assert fminus1.nc.variables["model_top"][0,0] == f.nc.variables["model_top"][0,0] - 1
-#     assert (f - fminus1).nc.variables["model_top"][0,0]==\
-#            f.nc.variables["model_top"][0,0] - \
-#            fminus1.nc.variables["model_top"][0,0]
-#
-#     ftimes2 = f * 2
-#     assert ftimes2.nc.variables["model_top"][0,0] == f.nc.variables["model_top"][0,0] * 2
-#     assert (f * ftimes2).nc.variables["model_top"][0,0] ==\
-#             f.nc.variables["model_top"][0,0] * \
-#            ftimes2.nc.variables["model_top"][0,0]
-#
-#     fdiv2 = f / 2
-#     assert fdiv2.nc.variables["model_top"][0,0] == f.nc.variables["model_top"][0,0] / 2
-#     assert (f / fdiv2).nc.variables["model_top"][0,0] == \
-#          f.nc.variables["model_top"][0,0] / \
-#            fdiv2.nc.variables["model_top"][0,0]
-#
-#     assert f.nc.variables["ibound"][0,0,0] == 1
 def test_wkt_parse():
     """Test parsing of Coordinate Reference System parameters
     from well-known-text in .prj files."""
@@ -1458,15 +1602,14 @@ def test_export_contourf():
 def main():
     # test_shapefile()
     # test_shapefile_ibound()
+    # test_netcdf_classmethods()
 
-    test_netcdf_classmethods()
-
-    for namfile in namfiles:
-        export_mf2005_netcdf(namfile)
-        export_shapefile(namfile)
+    # for namfile in namfiles:
+    #    export_mf2005_netcdf(namfile)
+    #    export_shapefile(namfile)
 
-    for namfile in namfiles[0:2]:
-        export_shapefile_modelgrid_override(namfile)
+    # for namfile in namfiles[0:2]:
+    #     export_shapefile_modelgrid_override(namfile)
 
     # test_netcdf_overloads()
     # test_netcdf_classmethods()
@@ -1480,7 +1623,7 @@ def main():
     # test_vertex_model_dot_plot()
     # test_sr_with_Map()
     # test_modelgrid_with_PlotMapView()
-    test_epsgs()
+    # test_epsgs()
     # test_sr_scaling()
     # test_read_usgs_model_reference()
     # test_dynamic_xll_yll()
@@ -1500,7 +1643,9 @@ def main():
     # test_export_contourf()
     # test_sr()
     # test_shapefile_polygon_closed()
-
+    test_mapview_plot_bc()
+    test_crosssection_plot_bc()
+    test_output_helper_shapefile_export()
 
 if __name__ == '__main__':
 

autotest/t012_test.py

@@ -61,7 +61,7 @@ def test_mf2000_p07():
     pth = os.path.join(pth2000, 'P07')
     namfile = 'p7mf2k.nam'
     mf = flopy.modflow.Modflow.load(namfile, model_ws=pth,
-                                    version='mf2k', verbose=True,
+                                    verbose=True,
                                     exe_name=mf2k_exe)
 
     cpth = os.path.join(newpth, 'P07_2K')

flopy/export/utils.py

@@ -278,6 +278,10 @@ def output_helper(f, ml, oudic, **kwargs):
         modelgrid : flopy.discretizaiton.Grid
             user supplied model grid instance that will be used for export
             in lieu of the models model grid instance
+        mflay : int
+            zero based model layer which can be used in shapefile exporting
+        kper : int
+            zero based stress period which can be used for shapefile exporting
 
     Returns
     -------
@@ -294,6 +298,8 @@ def output_helper(f, ml, oudic, **kwargs):
     forgive = kwargs.pop("forgive", False)
     kwargs.pop("suffix", None)
     mask_vals = []
+    mflay = kwargs.pop('mflay', None)
+    kper = kwargs.pop('kper', None)
     if "masked_vals" in kwargs:
         mask_vals = kwargs.pop("masked_vals")
     if len(kwargs) > 0 and logger is not None:
@@ -427,11 +433,65 @@ def output_helper(f, ml, oudic, **kwargs):
                                                    zonebud.flux,
                                                    logger)
 
-            # todo: write the zone array to standard output
+            # write the zone array to standard output
             _add_output_nc_variable(f, times, shape3d, zonebud,
                                     "budget_zones", logger=logger,
                                     mask_vals=mask_vals,
                                     mask_array3d=mask_array3d)
+
+    elif isinstance(f, str) and f.endswith('.shp'):
+        attrib_dict = {}
+        for _, out_obj in oudic.items():
+
+            if isinstance(out_obj, HeadFile) or \
+                    isinstance(out_obj, FormattedHeadFile) or \
+                    isinstance(out_obj, UcnFile):
+                if isinstance(out_obj, UcnFile):
+                    attrib_name = 'conc'
+                else:
+                    attrib_name = 'head'
+                plotarray = np.atleast_3d(out_obj.get_alldata()
+                                          .transpose()).transpose()
+
+
+                for per in range(plotarray.shape[0]):
+                    for k in range(plotarray.shape[1]):
+                        if kper is not None and per != kper:
+                            continue
+                        if mflay is not None and k != mflay:
+                            continue
+                        name = attrib_name + '{}_{}'.format(per, k)
+                        attrib_dict[name] = plotarray[per][k]
+
+            elif isinstance(out_obj, CellBudgetFile):
+                names = out_obj.get_unique_record_names(decode=True)
+
+                for attrib_name in names:
+                    plotarray = np.atleast_3d(out_obj.get_data(
+                                              text=attrib_name,
+                                              full3D=True))
+
+                    attrib_name = attrib_name.strip()
+                    if attrib_name == "FLOW RIGHT FACE":
+                        attrib_name = 'FRF'
+                    elif attrib_name == "FLOW FRONT FACE":
+                        attrib_name = "FFF"
+                    elif attrib_name == "FLOW LOWER FACE":
+                        attrib_name = "FLF"
+                    else:
+                        pass
+                    for per in range(plotarray.shape[0]):
+                        for k in range(plotarray.shape[1]):
+                            if kper is not None and per != kper:
+                                continue
+                            if mflay is not None and k != mflay:
+                                continue
+                            name = attrib_name + '{}_{}'.format(per, k)
+                            attrib_dict[name] = plotarray[per][k]
+
+        if attrib_dict:
+            shapefile_utils.write_grid_shapefile(f, ml.modelgrid, attrib_dict)
+
     else:
         if logger:
             logger.lraise("unrecognized export argument:{0}".format(f))

flopy/modflow/mf.py

@@ -682,6 +682,8 @@ def load(f, version='mf2005', exe_name='mf2005.exe', verbose=False,
         if 'NWT' in ext_pkg_d or 'UPW' in ext_pkg_d:
             version = 'mfnwt'
         if 'GLOBAL' in ext_pkg_d:
+            if version != "mf2k":
+                ml.glo = ModflowGlobal(ml)
             version = 'mf2k'
         if 'SMS' in ext_pkg_d:
             version = 'mfusg'

flopy/pakbase.py

@@ -502,8 +502,8 @@ def parent(self, parent):
 
     @property
     def package_type(self):
-        if len(self.extension) > 0:
-            return self.extension[0]
+        if len(self.name) > 0:
+            return self.name[0].lower()
 
     @property
     def plotable(self):