feat(lak6): support none lake bedleak values (#1189)

Also modified gridintersect.plot_linestring to accept a matplotlib
colormap string.

autotest/t004_test_utilarray.py

@@ -896,22 +896,33 @@ def test_mflist_fromfile():
     import pandas as pd
     import flopy
 
-    wel_data = pd.DataFrame([(0,1,2, -50.0),(0,5,5, -50.0)],
-                 columns=['k', 'i', 'j', 'flux'])
-    wel_data.to_csv('wel_000.dat', index=False, sep=' ',header=False)
-
-    nwt_model = flopy.modflow.Modflow('nwt_testmodel',verbose=True)
-    dis = flopy.modflow.ModflowDis(nwt_model, nlay=1, nrow=10, ncol=10, delr=500.0,
-                        delc=500.0, top=100.0, botm=50.0)
-    wel = flopy.modflow.ModflowWel(nwt_model, stress_period_data={0: 'wel_000.dat'})
+    wel_data = pd.DataFrame(
+        [(0, 1, 2, -50.0), (0, 5, 5, -50.0)], columns=["k", "i", "j", "flux"]
+    )
+    wel_data.to_csv("wel_000.dat", index=False, sep=" ", header=False)
+
+    nwt_model = flopy.modflow.Modflow("nwt_testmodel", verbose=True)
+    dis = flopy.modflow.ModflowDis(
+        nwt_model,
+        nlay=1,
+        nrow=10,
+        ncol=10,
+        delr=500.0,
+        delc=500.0,
+        top=100.0,
+        botm=50.0,
+    )
+    wel = flopy.modflow.ModflowWel(
+        nwt_model, stress_period_data={0: "wel_000.dat"}
+    )
     flx_array = wel.stress_period_data.array["flux"][0]
-    for k,i,j,flx in zip(wel_data.k,wel_data.i,wel_data.j,wel_data.flux):
-        assert flx_array[k,i,j] == flx
+    for k, i, j, flx in zip(wel_data.k, wel_data.i, wel_data.j, wel_data.flux):
+        assert flx_array[k, i, j] == flx
 
 
 if __name__ == "__main__":
     # test_util3d_reset()
-    #test_mflist()
+    # test_mflist()
     test_mflist_fromfile()
     # test_new_get_file_entry()
     # test_arrayformat()

autotest/t078_lake_connections.py

@@ -612,7 +612,82 @@ def test_embedded_lak_prudic():
     return
 
 
+def test_embedded_lak_prudic_mixed():
+    lakebed_leakance = 1.0  # Lakebed leakance ($ft^{-1}$)
+    nlay = 8  # Number of layers
+    nrow = 36  # Number of rows
+    ncol = 23  # Number of columns
+    delr = float(405.665)  # Column width ($ft$)
+    delc = float(403.717)  # Row width ($ft$)
+    delv = 15.0  # Layer thickness ($ft$)
+    top = 100.0  # Top of the model ($ft$)
+
+    shape2d = (nrow, ncol)
+    shape3d = (nlay, nrow, ncol)
+
+    # load data from text files
+    data_ws = os.path.join("..", "examples", "data", "mf6_test")
+    fname = os.path.join(data_ws, "prudic2004t2_bot1.dat")
+    bot0 = np.loadtxt(fname)
+    botm = np.array(
+        [bot0]
+        + [
+            np.ones(shape2d, dtype=float) * (bot0 - (delv * k))
+            for k in range(1, nlay)
+        ]
+    )
+    fname = os.path.join(data_ws, "prudic2004t2_idomain1.dat")
+    idomain0 = np.loadtxt(fname, dtype=np.int32)
+    idomain = np.array(nlay * [idomain0], dtype=np.int32)
+    fname = os.path.join(data_ws, "prudic2004t2_lakibd.dat")
+    lakibd = np.loadtxt(fname, dtype=int)
+    lake_map = np.ones(shape3d, dtype=np.int32) * -1
+    lake_map[0, :, :] = lakibd[:, :] - 1
+
+    lakebed_leakance = np.zeros(shape2d, dtype=object)
+    idx = np.where(lake_map[0, :, :] == 0)
+    lakebed_leakance[idx] = "none"
+    idx = np.where(lake_map[0, :, :] == 1)
+    lakebed_leakance[idx] = 1.0
+    lakebed_leakance = lakebed_leakance.tolist()
+
+    # build StructuredGrid
+    model_grid = flopy.discretization.StructuredGrid(
+        nlay=nlay,
+        nrow=nrow,
+        ncol=ncol,
+        delr=np.ones(ncol, dtype=float) * delr,
+        delc=np.ones(nrow, dtype=float) * delc,
+        top=np.ones(shape2d, dtype=float) * top,
+        botm=botm,
+        idomain=idomain,
+    )
+
+    # test mixed lakebed leakance list
+    (_, _, connectiondata,) = flopy.mf6.utils.get_lak_connections(
+        model_grid,
+        lake_map,
+        idomain=idomain,
+        bedleak=lakebed_leakance,
+    )
+
+    # test the connections
+    for data in connectiondata:
+        lakeno, bedleak = data[0], data[4]
+        if lakeno == 0:
+            assert (
+                bedleak == "none"
+            ), "bedleak for lake 0 " "is not 'none' ({})".format(bedleak)
+        else:
+            assert (
+                bedleak == 1.0
+            ), "bedleak for lake 1 " "is not 1.0 ({})".format(bedleak)
+
+    return
+
+
 if __name__ == "__main__":
+    test_embedded_lak_prudic_mixed()
     test_base_run()
     test_lake()
     test_embedded_lak_ex01()

flopy/mf6/utils/lakpak_utils.py

@@ -1,7 +1,7 @@
 import numpy as np
 
 
-def get_lak_connections(modelgrid, lake_map, idomain=None, bedleak=0.1):
+def get_lak_connections(modelgrid, lake_map, idomain=None, bedleak=None):
     """
     Function to create lake package connection data from a zero-based
     integer array of lake numbers. If the shape of lake number array is
@@ -32,7 +32,10 @@ def get_lak_connections(modelgrid, lake_map, idomain=None, bedleak=0.1):
         bed leakance for lakes in the model domain. If bedleak is a float the
         same bed leakance is applied to each lake connection in the model.
         If bedleak is of size (nrow, ncol) or (ncpl) then all lake
-        connections for the cellid are given the same bed leakance value.
+        connections for the cellid are given the same bed leakance value. If
+        bedleak is None, lake conductance is only a function of aquifer
+        properties for all lakes. Can also pass mixed values as list or
+        ndarray of size (nrow, ncol) or (ncpl) with floats and 'none' strings.
 
     Returns
     -------
@@ -90,10 +93,20 @@ def get_lak_connections(modelgrid, lake_map, idomain=None, bedleak=0.1):
         )
 
     # convert bedleak to numpy array if necessary
-    if isinstance(bedleak, (float, int)):
+    if bedleak is None:
+        bedleak = np.chararray(shape2d, itemsize=4, unicode=True)
+        bedleak[:] = "none"
+    elif isinstance(bedleak, (float, int)):
         bedleak = np.ones(shape2d, dtype=float) * float(bedleak)
     elif isinstance(bedleak, (list, tuple)):
-        bedleak = np.array(bedleak, dtype=float)
+        bedleak = np.array(bedleak, dtype=object)
+
+    # check the dimensions of the bedleak array
+    if bedleak.shape != shape2d:
+        raise ValueError(
+            "shape of bedleak "
+            "({}) not equal to {}".format(bedleak.shape, shape2d)
+        )
 
     # get the model grid elevations and reset lake_map using idomain
     # if lake is embedded and in an inactive cell

flopy/utils/gridintersect.py

@@ -1481,7 +1481,7 @@ def plot_polygon(rec, ax=None, **kwargs):
         return ax
 
     @staticmethod
-    def plot_linestring(rec, ax=None, **kwargs):
+    def plot_linestring(rec, ax=None, cmap=None, **kwargs):
         """method to plot the linestring intersection results from the
         resulting numpy.recarray.
 
@@ -1494,6 +1494,8 @@ def plot_linestring(rec, ax=None, **kwargs):
             (the resulting shapes)
         ax : matplotlib.pyplot.axes, optional
             axes to plot onto, if not provided, creates a new figure
+        cmap : str
+            matplotlib colormap
         **kwargs:
             passed to the plot function
 
@@ -1509,13 +1511,24 @@ def plot_linestring(rec, ax=None, **kwargs):
         if ax is None:
             _, ax = plt.subplots()
 
+        specified_color = True
+        if "c" in kwargs:
+            c = kwargs.pop("c")
+        elif "color" in kwargs:
+            c = kwargs.pop("color")
+        else:
+            specified_color = False
+
+        if cmap is not None:
+            colormap = plt.get_cmap(cmap)
+            colors = colormap(np.linspace(0, 1, rec.shape[0]))
+
         for i, ishp in enumerate(rec.ixshapes):
-            if "c" in kwargs:
-                c = kwargs.pop("c")
-            elif "color" in kwargs:
-                c = kwargs.pop("color")
-            else:
-                c = "C{}".format(i % 10)
+            if not specified_color:
+                if cmap is None:
+                    c = "C{}".format(i % 10)
+                else:
+                    c = colors[i]
             if ishp.type == "MultiLineString":
                 for part in ishp:
                     ax.plot(part.xy[0], part.xy[1], ls="-", c=c, **kwargs)