diff --git a/autotest/t004_test_utilarray.py b/autotest/t004_test_utilarray.py index ad47b701d..71f552a9d 100644 --- a/autotest/t004_test_utilarray.py +++ b/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() diff --git a/autotest/t078_lake_connections.py b/autotest/t078_lake_connections.py index b5b768b9b..61f24c414 100644 --- a/autotest/t078_lake_connections.py +++ b/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() diff --git a/examples/Notebooks/flopy3_grid_intersection_demo.ipynb b/examples/Notebooks/flopy3_grid_intersection_demo.ipynb index 3d1240f68..3a988569f 100644 --- a/examples/Notebooks/flopy3_grid_intersection_demo.ipynb +++ b/examples/Notebooks/flopy3_grid_intersection_demo.ipynb @@ -50,7 +50,7 @@ "[Clang 10.0.0 ]\n", "numpy version: 1.19.2\n", "matplotlib version: 3.4.2\n", - "flopy version: 3.3.4\n" + "flopy version: 3.3.5\n" ] } ], @@ -148,7 +148,7 @@ { "data": { "text/plain": [ - "" + "" ] }, "execution_count": 4, @@ -223,7 +223,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "7.48 ms ± 364 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n" + "8.15 ms ± 532 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n" ] } ], @@ -262,11 +262,11 @@ { "data": { "text/plain": [ - "rec.array([((2, 3), (((35.0, 80.0), (40.0, 76.66666666666667), (40.0, 70.0), (30.0, 70.0), (35.0, 80.0)),), 66.66666667, ),\n", - " ((2, 4), (((50.0, 70.0), (40.0, 70.0), (40.0, 76.66666666666667), (50.0, 70.0)),), 33.33333333, ),\n", - " ((3, 2), (((30.0, 70.0), (30.0, 60.0), (25.0, 60.0), (30.0, 70.0)),), 25. , ),\n", - " ((3, 3), (((40.0, 70.0), (40.0, 60.0), (30.0, 60.0), (30.0, 70.0), (40.0, 70.0)),), 100. , ),\n", - " ((3, 4), (((50.0, 70.0), (50.0, 60.0), (40.0, 60.0), (40.0, 70.0), (50.0, 70.0)),), 100. , )],\n", + "rec.array([((2, 3), (((35.0, 80.0), (40.0, 76.66666666666667), (40.0, 70.0), (30.0, 70.0), (35.0, 80.0)),), 66.66666667, ),\n", + " ((2, 4), (((50.0, 70.0), (40.0, 70.0), (40.0, 76.66666666666667), (50.0, 70.0)),), 33.33333333, ),\n", + " ((3, 2), (((30.0, 70.0), (30.0, 60.0), (25.0, 60.0), (30.0, 70.0)),), 25. , ),\n", + " ((3, 3), (((40.0, 70.0), (40.0, 60.0), (30.0, 60.0), (30.0, 70.0), (40.0, 70.0)),), 100. , ),\n", + " ((3, 4), (((50.0, 70.0), (50.0, 60.0), (40.0, 60.0), (40.0, 70.0), (50.0, 70.0)),), 100. , )],\n", " dtype=[('cellids', 'O'), ('vertices', 'O'), ('areas', '),\n", - " ((2, 4), (((50.0, 70.0), (40.0, 70.0), (40.0, 76.66666666666667), (50.0, 70.0)),), 33.33333333, ),\n", - " ((3, 2), (((30.0, 70.0), (30.0, 60.0), (25.0, 60.0), (30.0, 70.0)),), 25. , ),\n", - " ((3, 3), (((40.0, 70.0), (40.0, 60.0), (30.0, 60.0), (30.0, 70.0), (40.0, 70.0)),), 100. , ),\n", - " ((3, 4), (((50.0, 70.0), (50.0, 60.0), (40.0, 60.0), (40.0, 70.0), (50.0, 70.0)),), 100. , )],\n", + "rec.array([((2, 3), (((35.0, 80.0), (40.0, 76.66666666666667), (40.0, 70.0), (30.0, 70.0), (35.0, 80.0)),), 66.66666667, ),\n", + " ((2, 4), (((50.0, 70.0), (40.0, 70.0), (40.0, 76.66666666666667), (50.0, 70.0)),), 33.33333333, ),\n", + " ((3, 2), (((30.0, 70.0), (30.0, 60.0), (25.0, 60.0), (30.0, 70.0)),), 25. , ),\n", + " ((3, 3), (((40.0, 70.0), (40.0, 60.0), (30.0, 60.0), (30.0, 70.0), (40.0, 70.0)),), 100. , ),\n", + " ((3, 4), (((50.0, 70.0), (50.0, 60.0), (40.0, 60.0), (40.0, 70.0), (50.0, 70.0)),), 100. , )],\n", " dtype=[('cellids', 'O'), ('vertices', 'O'), ('areas', '" ] @@ -577,7 +577,7 @@ "source": [ "fig, ax = plt.subplots(1, 1, figsize=(8, 8))\n", "sgr.plot(ax=ax)\n", - "ix.plot_linestring(result, ax=ax)\n", + "ix.plot_linestring(result, ax=ax, cmap=\"viridis\")\n", "\n", "for irow, icol in result.cellids:\n", " h2, = ax.plot(sgr.xcellcenters[0, icol], sgr.ycellcenters[irow, 0], \"kx\", label=\"centroids of intersected gridcells\")\n", @@ -610,7 +610,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "5.41 ms ± 236 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n" + "5.66 ms ± 370 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n" ] } ], @@ -646,7 +646,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "792 µs ± 34.7 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n" + "793 µs ± 24 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n" ] } ], @@ -662,9 +662,9 @@ { "data": { "text/plain": [ - "rec.array([((5, 4), ((45.0, 45.0),), ),\n", - " ((8, 0), ((10.0, 10.0),), ),\n", - " ((9, 4), ((50.0, 0.0),), )],\n", + "rec.array([((5, 4), ((45.0, 45.0),), ),\n", + " ((8, 0), ((10.0, 10.0),), ),\n", + " ((9, 4), ((50.0, 0.0),), )],\n", " dtype=[('cellids', 'O'), ('vertices', 'O'), ('ixshapes', 'O')])" ] }, @@ -732,7 +732,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "325 µs ± 11.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n" + "319 µs ± 11.8 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n" ] } ], @@ -748,9 +748,9 @@ { "data": { "text/plain": [ - "rec.array([((9, 4), ),\n", - " ((5, 4), ),\n", - " ((8, 0), )],\n", + "rec.array([((9, 4), ),\n", + " ((5, 4), ),\n", + " ((8, 0), )],\n", " dtype=[('cellids', 'O'), ('ixshapes', 'O')])" ] }, @@ -804,7 +804,7 @@ { "data": { "text/plain": [ - "" + "" ] }, "execution_count": 31, @@ -956,8 +956,8 @@ { "data": { "text/plain": [ - "rec.array([(1, ((10.0, 10.0), (45.0, 45.0)), ),\n", - " (4, ((50.0, 0.0),), )],\n", + "rec.array([(1, ((10.0, 10.0), (45.0, 45.0)), ),\n", + " (4, ((50.0, 0.0),), )],\n", " dtype=[('cellids', 'O'), ('vertices', 'O'), ('ixshapes', 'O')])" ] }, diff --git a/flopy/mf6/utils/lakpak_utils.py b/flopy/mf6/utils/lakpak_utils.py index 5d4d16c26..332da3c4c 100644 --- a/flopy/mf6/utils/lakpak_utils.py +++ b/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 diff --git a/flopy/utils/gridintersect.py b/flopy/utils/gridintersect.py index f217df436..bfe70ae83 100644 --- a/flopy/utils/gridintersect.py +++ b/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)