From 2b24104e47f7b4d57aa8b863ed7a8750d02bc7f1 Mon Sep 17 00:00:00 2001 From: w-bonelli Date: Thu, 23 Mar 2023 10:28:03 -0400 Subject: [PATCH] feat(PRT): add utilities to configure/plot PRT models * to_coords() and to_prt() methods on particle data classes * to_mp7_[pathlines/endpoints]() dtype conversion functions * support plotting PRT pathlines in map/cross-section plots --- .github/workflows/commit.yml | 48 +- .github/workflows/examples.yml | 8 +- .github/workflows/mf6.yml | 31 +- .github/workflows/rtd.yml | 16 +- autotest/regression/test_mf6.py | 2 +- autotest/test_binaryfile.py | 1 + ...test_cross_section_line_representations.py | 99 +++ autotest/test_export.py | 1 - autotest/test_model_dot_plot.py | 90 +++ autotest/test_model_splitter.py | 60 +- autotest/test_modflow.py | 2 +- autotest/test_particledata.py | 470 ++++++++++- autotest/test_plot.py | 735 ------------------ autotest/test_plot_cross_section.py | 86 ++ autotest/test_plot_map_view.py | 200 +++++ autotest/test_plot_pathlines.py | 379 +++++++++ autotest/test_plot_quasi3d.py | 154 ++++ autotest/test_plotutil.py | 272 +++++++ flopy/modpath/mp7particledata.py | 380 ++++++++- flopy/plot/crosssection.py | 133 +++- flopy/plot/map.py | 107 ++- flopy/plot/plotutil.py | 216 +++++ flopy/utils/binaryfile.py | 4 +- 23 files changed, 2601 insertions(+), 893 deletions(-) create mode 100644 autotest/test_cross_section_line_representations.py create mode 100644 autotest/test_model_dot_plot.py delete mode 100644 autotest/test_plot.py create mode 100644 autotest/test_plot_cross_section.py create mode 100644 autotest/test_plot_map_view.py create mode 100644 autotest/test_plot_pathlines.py create mode 100644 autotest/test_plot_quasi3d.py create mode 100644 autotest/test_plotutil.py diff --git a/.github/workflows/commit.yml b/.github/workflows/commit.yml index 79decdc901..8fec69f0e1 100644 --- a/.github/workflows/commit.yml +++ b/.github/workflows/commit.yml @@ -124,14 +124,28 @@ jobs: run: | pip install --upgrade pip pip install . - pip install ".[test,optional]" - - - name: Install Modflow executables + pip install ".[test, optional]" + + - name: Install executables uses: modflowpy/install-modflow-action@v1 - - - name: Smoke test - working-directory: autotest - run: pytest -v -n=auto --smoke --cov=flopy --cov-report=xml --durations=0 --keep-failed=.failed + + # todo restore below before merging + - name: Install MODFLOW 6 PRT build + uses: modflowpy/install-modflow-action@v1 + with: + # repo: modflow6-nightly-build + owner: aprovost-usgs + repo: modflow6 + + - name: Update FloPy packages + if: runner.os != 'Windows' + # run: python -c 'import flopy; flopy.mf6.utils.generate_classes(ref="develop", backup=False)' + run: python -c 'import flopy; flopy.mf6.utils.generate_classes(owner="aprovost-usgs", ref="PRT", backup=False)' + + - name: Run smoke tests + working-directory: ./autotest + run: | + pytest -v -n=auto --smoke --durations=0 --keep-failed=.failed env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} @@ -150,7 +164,10 @@ jobs: test: name: Test - needs: smoke + needs: + - build + - lint + - smoke runs-on: ${{ matrix.os }} strategy: fail-fast: false @@ -205,22 +222,27 @@ jobs: pip install xmipy pip install . - - name: Install Modflow-related executables + - name: Install executables uses: modflowpy/install-modflow-action@v1 - - name: Install Modflow dev build executables + # todo restore below when prt merges to mf6 develop + - name: Install MODFLOW 6 PRT build uses: modflowpy/install-modflow-action@v1 with: - repo: modflow6-nightly-build + # repo: modflow6-nightly-build + owner: aprovost-usgs + repo: modflow6 - name: Update FloPy packages if: runner.os != 'Windows' - run: python -m flopy.mf6.utils.generate_classes --ref develop --no-backup + # run: python -c 'import flopy; flopy.mf6.utils.generate_classes(ref="develop", backup=False)' + run: python -c 'import flopy; flopy.mf6.utils.generate_classes(owner="aprovost-usgs", ref="PRT", backup=False)' - name: Update FloPy packages if: runner.os == 'Windows' shell: bash -l {0} - run: python -m flopy.mf6.utils.generate_classes --ref develop --no-backup + # run: python -c 'import flopy; flopy.mf6.utils.generate_classes(ref="develop", backup=False)' + run: python -c 'import flopy; flopy.mf6.utils.generate_classes(owner="aprovost-usgs", ref="PRT", backup=False)' - name: Run tests if: runner.os != 'Windows' diff --git a/.github/workflows/examples.yml b/.github/workflows/examples.yml index 62df8cf843..6c2cddf293 100644 --- a/.github/workflows/examples.yml +++ b/.github/workflows/examples.yml @@ -53,7 +53,7 @@ jobs: bash powershell - - name: Install extra Python dependencies + - name: Install Python dependencies if: runner.os == 'Windows' shell: bash -l {0} run: | @@ -96,8 +96,7 @@ jobs: - name: Run example tests if: runner.os != 'Windows' working-directory: ./autotest - run: | - pytest -v -m="example" -n=auto -s --durations=0 --keep-failed=.failed + run: pytest -v -m="example" -n=auto -s --durations=0 --keep-failed=.failed env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} @@ -105,8 +104,7 @@ jobs: if: runner.os == 'Windows' shell: bash -l {0} working-directory: ./autotest - run: | - pytest -v -m="example" -n=auto -s --durations=0 --keep-failed=.failed + run: pytest -v -m="example" -n=auto -s --durations=0 --keep-failed=.failed env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} diff --git a/.github/workflows/mf6.yml b/.github/workflows/mf6.yml index f00b13fe1b..2b74c7c69e 100644 --- a/.github/workflows/mf6.yml +++ b/.github/workflows/mf6.yml @@ -42,6 +42,7 @@ jobs: pip install https://github.com/modflowpy/pymake/zipball/master pip install https://github.com/Deltares/xmipy/zipball/develop pip install https://github.com/MODFLOW-USGS/modflowapi/zipball/develop + pip install . pip install .[test,optional] - name: Install gfortran @@ -50,53 +51,45 @@ jobs: - name: Checkout MODFLOW 6 uses: actions/checkout@v3 with: - repository: MODFLOW-USGS/modflow6 + repository: aprovost-usgs/modflow6 path: modflow6 + ref: PRT - name: Update flopy MODFLOW 6 classes working-directory: modflow6/autotest - run: | - python update_flopy.py + run: python update_flopy.py - name: Install meson - run: | - pip3 install meson ninja + run: pip3 install meson ninja - name: Setup modflow working-directory: modflow6 - run: | - meson setup builddir --buildtype=debugoptimized --prefix=$(pwd) --libdir=bin + run: meson setup builddir --buildtype=debugoptimized --prefix=$(pwd) --libdir=bin - name: Build modflow working-directory: modflow6 - run: | - meson compile -C builddir + run: meson compile -C builddir - name: Install modflow working-directory: modflow6 - run: | - meson install -C builddir + run: meson install -C builddir - name: Get executables working-directory: modflow6/autotest env: GITHUB_TOKEN: ${{ github.token }} - run: | - pytest -v --durations=0 get_exes.py + run: pytest -v --durations=0 get_exes.py - name: Run tests working-directory: modflow6/autotest - run: | - pytest -v -n auto -k "test_gw" --durations=0 --cov=flopy --cov-report=xml + run: pytest -v -n auto -k "test_gw" --durations=0 --cov=flopy --cov-report=xml - name: Print coverage report before upload working-directory: ./modflow6/autotest - run: | - coverage report + run: coverage report - name: Upload coverage to Codecov - if: - github.repository_owner == 'modflowpy' && (github.event_name == 'push' || github.event_name == 'pull_request') + if: github.repository_owner == 'modflowpy' && (github.event_name == 'push' || github.event_name == 'pull_request') uses: codecov/codecov-action@v3 with: files: ./modflow6/autotest/coverage.xml diff --git a/.github/workflows/rtd.yml b/.github/workflows/rtd.yml index 5df817f687..72f977ef87 100644 --- a/.github/workflows/rtd.yml +++ b/.github/workflows/rtd.yml @@ -49,7 +49,9 @@ jobs: run: python -m pip install --upgrade pip - name: Install flopy and dependencies - run: pip install ".[test, doc, optional]" + run: | + pip install . + pip install ".[test, doc, optional]" - name: Workaround OpenGL issue on Linux if: runner.os == 'Linux' @@ -77,6 +79,18 @@ jobs: - name: Install MODFLOW executables uses: modflowpy/install-modflow-action@v1 + # todo restore below before merging + - name: Install MODFLOW 6 PRT build + uses: modflowpy/install-modflow-action@v1 + with: + owner: aprovost-usgs + repo: modflow6 + + - name: Update FloPy + run: | + python -c 'from flopy.mf6.utils import generate_classes; generate_classes(owner="aprovost-usgs", ref="PRT", backup=False)' + pip install . + - name: Run tutorial and example notebooks working-directory: autotest run: pytest -v -n auto test_notebooks.py diff --git a/autotest/regression/test_mf6.py b/autotest/regression/test_mf6.py index 5479014475..5fc98c3fb8 100644 --- a/autotest/regression/test_mf6.py +++ b/autotest/regression/test_mf6.py @@ -3085,7 +3085,7 @@ def test028_create_tests_sfr(function_tmpdir, example_data_path): delc=5000.0, top=top, botm=botm, - #idomain=idomain, + # idomain=idomain, filename=f"{model_name}.dis", ) strt = testutils.read_std_array(os.path.join(pth, "strt.txt"), "float") diff --git a/autotest/test_binaryfile.py b/autotest/test_binaryfile.py index d0d5925205..d58c63e734 100644 --- a/autotest/test_binaryfile.py +++ b/autotest/test_binaryfile.py @@ -1,4 +1,5 @@ import os +from typing import List import numpy as np import pytest diff --git a/autotest/test_cross_section_line_representations.py b/autotest/test_cross_section_line_representations.py new file mode 100644 index 0000000000..15fbf89770 --- /dev/null +++ b/autotest/test_cross_section_line_representations.py @@ -0,0 +1,99 @@ +import numpy as np +import pytest +from modflow_devtools.markers import requires_pkg + +import flopy + + +def structured_square_grid(side: int = 10, thick: int = 10): + """ + Creates a basic 1-layer structured grid with the given thickness and number of cells per side + Parameters + ---------- + side : The number of cells per side + thick : The thickness of the grid's single layer + Returns + ------- + A single-layer StructuredGrid of the given size and thickness + """ + + from flopy.discretization.structuredgrid import StructuredGrid + + delr = np.ones(side) + delc = np.ones(side) + top = np.ones((side, side)) * thick + botm = np.ones((side, side)) * (top - thick).reshape(1, side, side) + return StructuredGrid(delr=delr, delc=delc, top=top, botm=botm) + + +@requires_pkg("shapely") +@pytest.mark.parametrize( + "line", + [(), [], (()), [[]], (0, 0), [0, 0], [[0, 0]]], +) +def test_cross_section_invalid_lines_raise_error(line): + grid = structured_square_grid(side=10) + with pytest.raises(ValueError): + flopy.plot.PlotCrossSection(modelgrid=grid, line={"line": line}) + + +@requires_pkg("shapely") +@pytest.mark.parametrize( + "line", + [ + # diagonal + [(0, 0), (10, 10)], + ([0, 0], [10, 10]), + # horizontal + ([0, 5.5], [10, 5.5]), + [(0, 5.5), (10, 5.5)], + # vertical + [(5.5, 0), (5.5, 10)], + ([5.5, 0], [5.5, 10]), + # multiple segments + [(0, 0), (4, 6), (10, 10)], + ([0, 0], [4, 6], [10, 10]), + ], +) +def test_cross_section_valid_line_representations(line): + from shapely.geometry import LineString as SLS + + from flopy.utils.geometry import LineString as FLS + + grid = structured_square_grid(side=10) + + fls = FLS(line) + sls = SLS(line) + + # use raw, flopy.utils.geometry and shapely.geometry representations + lxc = flopy.plot.PlotCrossSection(modelgrid=grid, line={"line": line}) + fxc = flopy.plot.PlotCrossSection(modelgrid=grid, line={"line": fls}) + sxc = flopy.plot.PlotCrossSection(modelgrid=grid, line={"line": sls}) + + # make sure parsed points are identical for all line representations + assert np.allclose(lxc.pts, fxc.pts) and np.allclose(lxc.pts, sxc.pts) + assert ( + set(lxc.xypts.keys()) == set(fxc.xypts.keys()) == set(sxc.xypts.keys()) + ) + for k in lxc.xypts.keys(): + assert np.allclose(lxc.xypts[k], fxc.xypts[k]) and np.allclose( + lxc.xypts[k], sxc.xypts[k] + ) + + +@pytest.mark.parametrize( + "line", + [ + 0, + [0], + [0, 0], + (0, 0), + [(0, 0)], + ([0, 0]), + ], +) +@requires_pkg("shapely", "geojson") +def test_cross_section_invalid_line_representations_fail(line): + grid = structured_square_grid(side=10) + with pytest.raises(ValueError): + flopy.plot.PlotCrossSection(modelgrid=grid, line={"line": line}) diff --git a/autotest/test_export.py b/autotest/test_export.py index afe9bdf39a..82584537ff 100644 --- a/autotest/test_export.py +++ b/autotest/test_export.py @@ -51,7 +51,6 @@ from flopy.utils.crs import get_authority_crs from flopy.utils.geometry import Polygon - HAS_PYPROJ = has_pkg("pyproj", strict=True) if HAS_PYPROJ: import pyproj diff --git a/autotest/test_model_dot_plot.py b/autotest/test_model_dot_plot.py new file mode 100644 index 0000000000..ab7f261b77 --- /dev/null +++ b/autotest/test_model_dot_plot.py @@ -0,0 +1,90 @@ +import os + +import pytest +from flaky import flaky +from matplotlib import pyplot as plt +from matplotlib import rcParams +from matplotlib.collections import ( + LineCollection, + PatchCollection, + PathCollection, + QuadMesh, +) + +import flopy +from flopy.mf6 import MFSimulation + + +@pytest.mark.mf6 +def test_vertex_model_dot_plot(example_data_path): + rcParams["figure.max_open_warning"] = 36 + + # load up the vertex example problem + sim = MFSimulation.load( + sim_ws=example_data_path / "mf6" / "test003_gwftri_disv" + ) + disv_ml = sim.get_model("gwf_1") + ax = disv_ml.plot() + assert isinstance(ax, list) + assert len(ax) == 36 + + +# occasional _tkinter.TclError: Can't find a usable tk.tcl (or init.tcl) +# similar: https://github.com/microsoft/azure-pipelines-tasks/issues/16426 +@flaky +def test_model_dot_plot(function_tmpdir, example_data_path): + loadpth = example_data_path / "mf2005_test" + ml = flopy.modflow.Modflow.load( + "ibs2k.nam", "mf2k", model_ws=loadpth, check=False + ) + ax = ml.plot() + assert isinstance(ax, list), "ml.plot() ax is is not a list" + assert len(ax) == 18, f"number of axes ({len(ax)}) is not equal to 18" + + +def test_dataset_dot_plot(function_tmpdir, example_data_path): + loadpth = example_data_path / "mf2005_test" + ml = flopy.modflow.Modflow.load( + "ibs2k.nam", "mf2k", model_ws=loadpth, check=False + ) + + # plot specific dataset + ax = ml.bcf6.hy.plot() + assert isinstance(ax, list), "ml.bcf6.hy.plot() ax is is not a list" + assert len(ax) == 2, f"number of hy axes ({len(ax)}) is not equal to 2" + + +def test_dataset_dot_plot_nlay_ne_plottable( + function_tmpdir, example_data_path +): + import matplotlib.pyplot as plt + + loadpth = example_data_path / "mf2005_test" + ml = flopy.modflow.Modflow.load( + "ibs2k.nam", "mf2k", model_ws=loadpth, check=False + ) + # special case where nlay != plottable + ax = ml.bcf6.vcont.plot() + assert isinstance( + ax, plt.Axes + ), "ml.bcf6.vcont.plot() ax is is not of type plt.Axes" + + +def test_model_dot_plot_export(function_tmpdir, example_data_path): + loadpth = example_data_path / "mf2005_test" + ml = flopy.modflow.Modflow.load( + "ibs2k.nam", "mf2k", model_ws=loadpth, check=False + ) + + fh = os.path.join(function_tmpdir, "ibs2k") + ml.plot(mflay=0, filename_base=fh, file_extension="png") + files = [f for f in os.listdir(function_tmpdir) if f.endswith(".png")] + if len(files) < 10: + raise AssertionError( + "ml.plot did not properly export all supported data types" + ) + + for f in files: + t = f.split("_") + if len(t) < 3: + raise AssertionError("Plot filenames not written correctly") diff --git a/autotest/test_model_splitter.py b/autotest/test_model_splitter.py index 921c29dc4a..e743399d91 100644 --- a/autotest/test_model_splitter.py +++ b/autotest/test_model_splitter.py @@ -1,9 +1,9 @@ import numpy as np -import flopy import pytest from autotest.conftest import get_example_data_path from modflow_devtools.markers import requires_exe, requires_pkg +import flopy from flopy.mf6 import MFSimulation from flopy.mf6.utils import Mf6Splitter @@ -260,11 +260,7 @@ def test_control_records(function_tmpdir): tdis = flopy.mf6.ModflowTdis( sim, nper=nper, - perioddata=( - (1., 1, 1.), - (1., 1, 1.), - (1., 1, 1.) - ) + perioddata=((1.0, 1, 1.0), (1.0, 1, 1.0), (1.0, 1, 1.0)), ) gwf = flopy.mf6.ModflowGwf(sim, save_flows=True) @@ -279,20 +275,20 @@ def test_control_records(function_tmpdir): delc=1, top=35, botm=[30, botm2], - idomain=1 + idomain=1, ) ic = flopy.mf6.ModflowGwfic(gwf, strt=32) npf = flopy.mf6.ModflowGwfnpf( gwf, k=[ - 1., + 1.0, { "data": np.ones((10, 10)) * 0.75, "filename": "k.l2.txt", "iprn": 1, - "factor": 1 - } + "factor": 1, + }, ], k33=[ np.ones((nrow, ncol)), @@ -301,41 +297,29 @@ def test_control_records(function_tmpdir): "filename": "k33.l2.bin", "iprn": 1, "factor": 1, - "binary": True - } - ] - + "binary": True, + }, + ], ) - wel_rec = [((0, 4, 5), -10), ] + wel_rec = [ + ((0, 4, 5), -10), + ] spd = { 0: wel_rec, - 1: { - "data": wel_rec, - "filename": "wel.1.txt" - }, - 2: { - "data": wel_rec, - "filename": "wel.2.bin", - "binary": True - } + 1: {"data": wel_rec, "filename": "wel.1.txt"}, + 2: {"data": wel_rec, "filename": "wel.2.bin", "binary": True}, } - wel = flopy.mf6.ModflowGwfwel( - gwf, - stress_period_data=spd - ) + wel = flopy.mf6.ModflowGwfwel(gwf, stress_period_data=spd) chd_rec = [] for cond, j in ((30, 0), (22, 9)): for i in range(10): chd_rec.append(((0, i, j), cond)) - chd = flopy.mf6.ModflowGwfchd( - gwf, - stress_period_data={0: chd_rec} - ) + chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data={0: chd_rec}) arr = np.zeros((10, 10), dtype=int) arr[0:5, :] = 1 @@ -354,14 +338,18 @@ def test_control_records(function_tmpdir): "External ascii files not being preserved for MFArray" ) - k33ls = ml1.npf.k33._data_storage.layer_storage.multi_dim_list + k33ls = ml1.npf.k33._data_storage.layer_storage.multi_dim_list if k33ls[1].data_storage_type.value != 3 or not k33ls[1].binary: raise AssertionError( "Binary file input not being preserved for MFArray" ) - spd_ls1 = ml1.wel.stress_period_data._data_storage[1].layer_storage.multi_dim_list[0] - spd_ls2 = ml1.wel.stress_period_data._data_storage[2].layer_storage.multi_dim_list[0] + spd_ls1 = ml1.wel.stress_period_data._data_storage[ + 1 + ].layer_storage.multi_dim_list[0] + spd_ls2 = ml1.wel.stress_period_data._data_storage[ + 2 + ].layer_storage.multi_dim_list[0] if spd_ls1.data_storage_type.value != 3 or spd_ls1.binary: raise AssertionError( @@ -372,5 +360,3 @@ def test_control_records(function_tmpdir): raise AssertionError( "External binary file input not being preseved for MFList" ) - - diff --git a/autotest/test_modflow.py b/autotest/test_modflow.py index d19a255d37..fd34810f00 100644 --- a/autotest/test_modflow.py +++ b/autotest/test_modflow.py @@ -7,8 +7,8 @@ import numpy as np import pytest from autotest.conftest import get_example_data_path -from modflow_devtools.misc import has_pkg from modflow_devtools.markers import excludes_platform, requires_exe +from modflow_devtools.misc import has_pkg from flopy.discretization import StructuredGrid from flopy.mf6 import MFSimulation diff --git a/autotest/test_particledata.py b/autotest/test_particledata.py index 1f5bdc85fd..968f935077 100644 --- a/autotest/test_particledata.py +++ b/autotest/test_particledata.py @@ -1,8 +1,50 @@ +from itertools import chain + +import matplotlib.pyplot as plt import numpy as np +import pandas as pd +import pytest +from autotest.test_grid_cases import GridCases + +from flopy.discretization import StructuredGrid +from flopy.mf6.modflow.mfgwf import ModflowGwf +from flopy.mf6.modflow.mfgwfdisv import ModflowGwfdisv +from flopy.mf6.modflow.mfgwfic import ModflowGwfic +from flopy.mf6.modflow.mfgwfnpf import ModflowGwfnpf +from flopy.mf6.modflow.mfgwfoc import ModflowGwfoc +from flopy.mf6.modflow.mfims import ModflowIms +from flopy.mf6.modflow.mfsimulation import MFSimulation +from flopy.mf6.modflow.mftdis import ModflowTdis +from flopy.modpath import ( + CellDataType, + FaceDataType, + LRCParticleData, + Modpath7, + Modpath7Bas, + Modpath7Sim, + NodeParticleData, + ParticleData, + ParticleGroupNodeTemplate, +) +from flopy.utils.modpathfile import PathlineFile -from flopy.modpath import ParticleData +# utilities -structured_plocs = [(1, 1, 1), (1, 1, 2)] + +def flatten(a): + return [ + [ + *chain.from_iterable( + xx if isinstance(xx, tuple) else [xx] for xx in x + ) + ] + for x in a + ] + + +# basic structured grid info + +structured_cells = [(0, 1, 1), (0, 1, 2)] structured_dtype = np.dtype( [ ("k", " 0, "Boundary condition was not drawn" - - mapview2 = flopy.plot.PlotMapView(model=ml6c, ax=mapview.ax) - mapview2.plot_bc("MAW") - ax = mapview2.ax - - assert len(ax.collections) > 0, "Boundary condition was not drawn" - - for col in ax.collections: - assert isinstance( - col, (QuadMesh, PathCollection) - ), f"Unexpected collection type: {type(col)}" - - -@pytest.mark.mf6 -@pytest.mark.xfail(reason="sometimes get wrong collection type") -def test_map_view_bc_UZF_3lay(example_data_path): - mpath = example_data_path / "mf6" / "test001e_UZF_3lay" - sim = MFSimulation.load(sim_ws=mpath) - ml6 = sim.get_model("gwf_1") - - mapview = flopy.plot.PlotMapView(model=ml6) - mapview.plot_bc("UZF") - ax = mapview.ax - - if len(ax.collections) == 0: - raise AssertionError("Boundary condition was not drawn") - - for col in ax.collections: - assert isinstance( - col, (QuadMesh, PathCollection) - ), f"Unexpected collection type: {type(col)}" - - -@pytest.mark.mf6 -@pytest.mark.xfail( - reason="sometimes get LineCollections instead of PatchCollections" -) -def test_cross_section_bc_gwfs_disv(example_data_path): - mpath = example_data_path / "mf6" / "test003_gwfs_disv" - sim = MFSimulation.load(sim_ws=mpath) - 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 - - assert len(ax.collections) != 0, "Boundary condition was not drawn" - - for col in ax.collections: - assert isinstance( - col, PatchCollection - ), f"Unexpected collection type: {type(col)}" - - -@pytest.mark.mf6 -@pytest.mark.xfail( - reason="sometimes get LineCollections instead of PatchCollections" -) -def test_cross_section_bc_lake2tr(example_data_path): - mpath = example_data_path / "mf6" / "test045_lake2tr" - sim = MFSimulation.load(sim_ws=mpath) - ml6 = sim.get_model("lakeex2a") - xc = flopy.plot.PlotCrossSection(ml6, line={"row": 10}) - xc.plot_bc("LAK") - xc.plot_bc("SFR") - - ax = xc.ax - assert len(ax.collections) != 0, "Boundary condition was not drawn" - - for col in ax.collections: - assert isinstance( - col, PatchCollection - ), f"Unexpected collection type: {type(col)}" - - -@pytest.mark.mf6 -@pytest.mark.xfail( - reason="sometimes get LineCollections instead of PatchCollections" -) -def test_cross_section_bc_2models_mvr(example_data_path): - mpath = example_data_path / "mf6" / "test006_2models_mvr" - sim = MFSimulation.load(sim_ws=mpath) - ml6 = sim.get_model("parent") - xc = flopy.plot.PlotCrossSection(ml6, line={"column": 1}) - xc.plot_bc("MAW") - - ax = xc.ax - assert len(ax.collections) > 0, "Boundary condition was not drawn" - - for col in ax.collections: - assert isinstance( - col, PatchCollection - ), f"Unexpected collection type: {type(col)}" - - -@pytest.mark.mf6 -@pytest.mark.xfail( - reason="sometimes get LineCollections instead of PatchCollections" -) -def test_cross_section_bc_UZF_3lay(example_data_path): - mpath = example_data_path / "mf6" / "test001e_UZF_3lay" - sim = MFSimulation.load(sim_ws=mpath) - ml6 = sim.get_model("gwf_1") - - xc = flopy.plot.PlotCrossSection(ml6, line={"row": 0}) - xc.plot_bc("UZF") - - ax = xc.ax - assert len(ax.collections) != 0, "Boundary condition was not drawn" - - for col in ax.collections: - assert isinstance( - col, PatchCollection - ), f"Unexpected collection type: {type(col)}" - - -def test_map_view_contour(function_tmpdir): - arr = np.random.rand(10, 10) * 100 - arr[-1, :] = np.nan - delc = np.array([10] * 10, dtype=float) - delr = np.array([8] * 10, dtype=float) - top = np.ones((10, 10), dtype=float) - botm = np.ones((3, 10, 10), dtype=float) - botm[0] = 0.75 - botm[1] = 0.5 - botm[2] = 0.25 - idomain = np.ones((3, 10, 10)) - idomain[0, 0, :] = 0 - vmin = np.nanmin(arr) - vmax = np.nanmax(arr) - levels = np.linspace(vmin, vmax, 7) - - grid = StructuredGrid( - delc=delc, - delr=delr, - top=top, - botm=botm, - idomain=idomain, - lenuni=1, - nlay=3, - nrow=10, - ncol=10, - ) - - pmv = PlotMapView(modelgrid=grid, layer=0) - contours = pmv.contour_array(a=arr) - plt.savefig(function_tmpdir / "map_view_contour.png") - - # if we ever revert from standard contours to tricontours, restore this nan check - # for ix, lev in enumerate(contours.levels): - # if not np.allclose(lev, levels[ix]): - # raise AssertionError("TriContour NaN catch Failed") - - -@pytest.mark.mf6 -def test_vertex_model_dot_plot(example_data_path): - rcParams["figure.max_open_warning"] = 36 - - # load up the vertex example problem - sim = MFSimulation.load( - sim_ws=example_data_path / "mf6" / "test003_gwftri_disv" - ) - disv_ml = sim.get_model("gwf_1") - ax = disv_ml.plot() - assert isinstance(ax, list) - assert len(ax) == 36 - - -# occasional _tkinter.TclError: Can't find a usable tk.tcl (or init.tcl) -# similar: https://github.com/microsoft/azure-pipelines-tasks/issues/16426 -@flaky -def test_model_dot_plot(function_tmpdir, example_data_path): - loadpth = example_data_path / "mf2005_test" - ml = flopy.modflow.Modflow.load( - "ibs2k.nam", "mf2k", model_ws=loadpth, check=False - ) - ax = ml.plot() - assert isinstance(ax, list), "ml.plot() ax is is not a list" - assert len(ax) == 18, f"number of axes ({len(ax)}) is not equal to 18" - - -def test_dataset_dot_plot(function_tmpdir, example_data_path): - loadpth = example_data_path / "mf2005_test" - ml = flopy.modflow.Modflow.load( - "ibs2k.nam", "mf2k", model_ws=loadpth, check=False - ) - - # plot specific dataset - ax = ml.bcf6.hy.plot() - assert isinstance(ax, list), "ml.bcf6.hy.plot() ax is is not a list" - assert len(ax) == 2, f"number of hy axes ({len(ax)}) is not equal to 2" - - -def test_dataset_dot_plot_nlay_ne_plottable( - function_tmpdir, example_data_path -): - import matplotlib.pyplot as plt - - loadpth = example_data_path / "mf2005_test" - ml = flopy.modflow.Modflow.load( - "ibs2k.nam", "mf2k", model_ws=loadpth, check=False - ) - # special case where nlay != plottable - ax = ml.bcf6.vcont.plot() - assert isinstance( - ax, plt.Axes - ), "ml.bcf6.vcont.plot() ax is is not of type plt.Axes" - - -def test_model_dot_plot_export(function_tmpdir, example_data_path): - loadpth = example_data_path / "mf2005_test" - ml = flopy.modflow.Modflow.load( - "ibs2k.nam", "mf2k", model_ws=loadpth, check=False - ) - - fh = os.path.join(function_tmpdir, "ibs2k") - ml.plot(mflay=0, filename_base=fh, file_extension="png") - files = [f for f in os.listdir(function_tmpdir) if f.endswith(".png")] - if len(files) < 10: - raise AssertionError( - "ml.plot did not properly export all supported data types" - ) - - for f in files: - t = f.split("_") - if len(t) < 3: - raise AssertionError("Plot filenames not written correctly") - - -@pytest.fixture -def modpath_model(function_tmpdir, example_data_path): - # test with multi-layer example - load_ws = example_data_path / "mp6" - - ml = Modflow.load("EXAMPLE.nam", model_ws=load_ws, exe_name="mf2005") - ml.change_model_ws(function_tmpdir) - ml.write_input() - ml.run_model() - - mp = Modpath6( - modelname="ex6", - exe_name="mp6", - modflowmodel=ml, - model_ws=function_tmpdir, - ) - - mpb = Modpath6Bas( - mp, hdry=ml.lpf.hdry, laytyp=ml.lpf.laytyp, ibound=1, prsity=0.1 - ) - - sim = mp.create_mpsim( - trackdir="forward", - simtype="pathline", - packages="RCH", - start_time=(2, 0, 1.0), - ) - return ml, mp, sim - - -@requires_exe("mf2005", "mp6") -def test_plot_map_view_mp6_plot_pathline(modpath_model): - ml, mp, sim = modpath_model - mp.write_input() - mp.run_model(silent=False) - - pthobj = PathlineFile(os.path.join(mp.model_ws, "ex6.mppth")) - well_pathlines = pthobj.get_destination_pathline_data( - dest_cells=[(4, 12, 12)] - ) - - def test_plot(pl): - mx = PlotMapView(model=ml) - mx.plot_grid() - mx.plot_bc("WEL", kper=2, color="blue") - pth = mx.plot_pathline(pl, colors="red") - # plt.show() - assert isinstance(pth, LineCollection) - assert len(pth._paths) == 114 - - # support pathlines as list of recarrays - test_plot(well_pathlines) - - # support pathlines as list of dataframes - test_plot([pd.DataFrame(pl) for pl in well_pathlines]) - - # support pathlines as single recarray - test_plot(np.concatenate(well_pathlines)) - - # support pathlines as single dataframe - test_plot(pd.DataFrame(np.concatenate(well_pathlines))) - - -@requires_exe("mf2005", "mp6") -def test_plot_cross_section_mp6_plot_pathline(modpath_model): - ml, mp, sim = modpath_model - mp.write_input() - mp.run_model(silent=False) - - pthobj = PathlineFile(os.path.join(mp.model_ws, "ex6.mppth")) - well_pathlines = pthobj.get_destination_pathline_data( - dest_cells=[(4, 12, 12)] - ) - - def test_plot(pl): - mx = PlotCrossSection(model=ml, line={"row": 4}) - mx.plot_bc("WEL", kper=2, color="blue") - pth = mx.plot_pathline(pl, method="cell", colors="red") - assert isinstance(pth, LineCollection) - assert len(pth._paths) == 6 - - # support pathlines as list of recarrays - test_plot(well_pathlines) - - # support pathlines as list of dataframes - test_plot([pd.DataFrame(pl) for pl in well_pathlines]) - - # support pathlines as single recarray - test_plot(np.concatenate(well_pathlines)) - - # support pathlines as single dataframe - test_plot(pd.DataFrame(np.concatenate(well_pathlines))) - - -@requires_exe("mf2005", "mp6") -def test_plot_map_view_mp6_endpoint(modpath_model): - ml, mp, sim = modpath_model - mp.write_input() - mp.run_model(silent=False) - - pthobj = EndpointFile(os.path.join(mp.model_ws, "ex6.mpend")) - endpts = pthobj.get_alldata() - - # support endpoints as recarray - assert isinstance(endpts, np.recarray) - mv = PlotMapView(model=ml) - mv.plot_bc("WEL", kper=2, color="blue") - ep = mv.plot_endpoint(endpts, direction="ending") - # plt.show() - assert isinstance(ep, PathCollection) - - # support endpoints as dataframe - mv = PlotMapView(model=ml) - mv.plot_bc("WEL", kper=2, color="blue") - ep = mv.plot_endpoint(pd.DataFrame(endpts), direction="ending") - # plt.show() - assert isinstance(ep, PathCollection) - - # test various possibilities for endpoint color configuration. - # first, color kwarg as scalar - mv = PlotMapView(model=ml) - mv.plot_bc("WEL", kper=2, color="blue") - ep = mv.plot_endpoint(endpts, direction="ending", color="red") - # plt.show() - assert isinstance(ep, PathCollection) - - # c kwarg as array - mv = PlotMapView(model=ml) - mv.plot_bc("WEL", kper=2, color="blue") - ep = mv.plot_endpoint( - endpts, - direction="ending", - c=np.random.rand(625) * -1000, - cmap="viridis", - ) - # plt.show() - assert isinstance(ep, PathCollection) - - # colorbar: color by time to termination - mv = PlotMapView(model=ml) - mv.plot_bc("WEL", kper=2, color="blue") - ep = mv.plot_endpoint( - endpts, direction="ending", shrink=0.5, colorbar=True - ) - # plt.show() - assert isinstance(ep, PathCollection) - - # if both color and c are provided, c takes precedence - mv = PlotMapView(model=ml) - mv.plot_bc("WEL", kper=2, color="blue") - ep = mv.plot_endpoint( - endpts, - direction="ending", - color="red", - c=np.random.rand(625) * -1000, - cmap="viridis", - ) - # plt.show() - assert isinstance(ep, PathCollection) - - -@pytest.fixture -def quasi3d_model(function_tmpdir): - mf = Modflow("model_mf", model_ws=function_tmpdir, exe_name="mf2005") - - # Model domain and grid definition - Lx = 1000.0 - Ly = 1000.0 - ztop = 0.0 - zbot = -30.0 - nlay = 3 - nrow = 10 - ncol = 10 - delr = Lx / ncol - delc = Ly / nrow - laycbd = [0] * (nlay) - laycbd[0] = 1 - botm = np.linspace(ztop, zbot, nlay + np.sum(laycbd) + 1)[1:] - - # Create the discretization object - ModflowDis( - mf, - nlay, - nrow, - ncol, - delr=delr, - delc=delc, - top=ztop, - botm=botm, - laycbd=laycbd, - ) - - # Variables for the BAS package - ibound = np.ones((nlay, nrow, ncol), dtype=np.int32) - ibound[:, :, 0] = -1 - ibound[:, :, -1] = -1 - strt = np.ones((nlay, nrow, ncol), dtype=np.float32) - strt[:, :, 0] = 10.0 - strt[:, :, -1] = 0.0 - ModflowBas(mf, ibound=ibound, strt=strt) - - # Add LPF package to the MODFLOW model - ModflowLpf(mf, hk=10.0, vka=10.0, ipakcb=53, vkcb=10) - - # add a well - row = int((nrow - 1) / 2) - col = int((ncol - 1) / 2) - spd = {0: [[1, row, col, -1000]]} - ModflowWel(mf, stress_period_data=spd) - - # Add OC package to the MODFLOW model - spd = {(0, 0): ["save head", "save budget"]} - ModflowOc(mf, stress_period_data=spd, compact=True) - - # Add PCG package to the MODFLOW model - ModflowPcg(mf) - - # Write the MODFLOW model input files - mf.write_input() - - # Run the MODFLOW model - success, buff = mf.run_model() - - assert success, "test_plotting_with_quasi3d_layers() failed" - - return mf - - -@requires_exe("mf2005") -def test_map_plot_with_quasi3d_layers(quasi3d_model): - # read output - hf = HeadFile( - os.path.join(quasi3d_model.model_ws, f"{quasi3d_model.name}.hds") - ) - head = hf.get_data(totim=1.0) - cbb = CellBudgetFile( - os.path.join(quasi3d_model.model_ws, f"{quasi3d_model.name}.cbc") - ) - frf = cbb.get_data(text="FLOW RIGHT FACE", totim=1.0)[0] - fff = cbb.get_data(text="FLOW FRONT FACE", totim=1.0)[0] - flf = cbb.get_data(text="FLOW LOWER FACE", totim=1.0)[0] - - # plot a map - plt.figure() - mv = PlotMapView(model=quasi3d_model, layer=1) - mv.plot_grid() - mv.plot_array(head) - mv.contour_array(head) - mv.plot_ibound() - mv.plot_bc("wel") - mv.plot_vector(frf, fff) - plt.savefig(os.path.join(quasi3d_model.model_ws, "plt01.png")) - - -@requires_exe("mf2005") -def test_cross_section_with_quasi3d_layers(quasi3d_model): - # read output - hf = HeadFile( - os.path.join(quasi3d_model.model_ws, f"{quasi3d_model.name}.hds") - ) - head = hf.get_data(totim=1.0) - cbb = CellBudgetFile( - os.path.join(quasi3d_model.model_ws, f"{quasi3d_model.name}.cbc") - ) - frf = cbb.get_data(text="FLOW RIGHT FACE", totim=1.0)[0] - fff = cbb.get_data(text="FLOW FRONT FACE", totim=1.0)[0] - flf = cbb.get_data(text="FLOW LOWER FACE", totim=1.0)[0] - - # plot a cross-section - plt.figure() - cs = PlotCrossSection( - model=quasi3d_model, - line={"row": int((quasi3d_model.modelgrid.nrow - 1) / 2)}, - ) - cs.plot_grid() - cs.plot_array(head) - cs.contour_array(head) - cs.plot_ibound() - cs.plot_bc("wel") - cs.plot_vector(frf, fff, flf, head=head) - plt.savefig(os.path.join(quasi3d_model.model_ws, "plt02.png")) - plt.close() - - -def structured_square_grid(side: int = 10, thick: int = 10): - """ - Creates a basic 1-layer structured grid with the given thickness and number of cells per side - Parameters - ---------- - side : The number of cells per side - thick : The thickness of the grid's single layer - Returns - ------- - A single-layer StructuredGrid of the given size and thickness - """ - - from flopy.discretization.structuredgrid import StructuredGrid - - delr = np.ones(side) - delc = np.ones(side) - top = np.ones((side, side)) * thick - botm = np.ones((side, side)) * (top - thick).reshape(1, side, side) - return StructuredGrid(delr=delr, delc=delc, top=top, botm=botm) - - -@requires_pkg("shapely") -@pytest.mark.parametrize( - "line", - [(), [], (()), [[]], (0, 0), [0, 0], [[0, 0]]], -) -def test_cross_section_invalid_lines_raise_error(line): - grid = structured_square_grid(side=10) - with pytest.raises(ValueError): - flopy.plot.PlotCrossSection(modelgrid=grid, line={"line": line}) - - -@requires_pkg("shapely") -@pytest.mark.parametrize( - "line", - [ - # diagonal - [(0, 0), (10, 10)], - ([0, 0], [10, 10]), - # horizontal - ([0, 5.5], [10, 5.5]), - [(0, 5.5), (10, 5.5)], - # vertical - [(5.5, 0), (5.5, 10)], - ([5.5, 0], [5.5, 10]), - # multiple segments - [(0, 0), (4, 6), (10, 10)], - ([0, 0], [4, 6], [10, 10]), - ], -) -def test_cross_section_valid_line_representations(line): - from shapely.geometry import LineString as SLS - - from flopy.utils.geometry import LineString as FLS - - grid = structured_square_grid(side=10) - - fls = FLS(line) - sls = SLS(line) - - # use raw, flopy.utils.geometry and shapely.geometry representations - lxc = flopy.plot.PlotCrossSection(modelgrid=grid, line={"line": line}) - fxc = flopy.plot.PlotCrossSection(modelgrid=grid, line={"line": fls}) - sxc = flopy.plot.PlotCrossSection(modelgrid=grid, line={"line": sls}) - - # make sure parsed points are identical for all line representations - assert np.allclose(lxc.pts, fxc.pts) and np.allclose(lxc.pts, sxc.pts) - assert ( - set(lxc.xypts.keys()) == set(fxc.xypts.keys()) == set(sxc.xypts.keys()) - ) - for k in lxc.xypts.keys(): - assert np.allclose(lxc.xypts[k], fxc.xypts[k]) and np.allclose( - lxc.xypts[k], sxc.xypts[k] - ) - - -@pytest.mark.parametrize( - "line", - [ - 0, - [0], - [0, 0], - (0, 0), - [(0, 0)], - ([0, 0]), - ], -) -@requires_pkg("shapely", "geojson") -def test_cross_section_invalid_line_representations_fail(line): - grid = structured_square_grid(side=10) - with pytest.raises(ValueError): - flopy.plot.PlotCrossSection(modelgrid=grid, line={"line": line}) diff --git a/autotest/test_plot_cross_section.py b/autotest/test_plot_cross_section.py new file mode 100644 index 0000000000..d4a507e928 --- /dev/null +++ b/autotest/test_plot_cross_section.py @@ -0,0 +1,86 @@ +import pytest + +import flopy +from flopy.mf6 import MFSimulation + + +@pytest.mark.mf6 +@pytest.mark.xfail( + reason="sometimes get LineCollections instead of PatchCollections" +) +def test_cross_section_bc_gwfs_disv(example_data_path): + mpath = example_data_path / "mf6" / "test003_gwfs_disv" + sim = MFSimulation.load(sim_ws=mpath) + 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 + + assert len(ax.collections) != 0, "Boundary condition was not drawn" + + for col in ax.collections: + assert isinstance( + col, PatchCollection + ), f"Unexpected collection type: {type(col)}" + + +@pytest.mark.mf6 +@pytest.mark.xfail( + reason="sometimes get LineCollections instead of PatchCollections" +) +def test_cross_section_bc_lake2tr(example_data_path): + mpath = example_data_path / "mf6" / "test045_lake2tr" + sim = MFSimulation.load(sim_ws=mpath) + ml6 = sim.get_model("lakeex2a") + xc = flopy.plot.PlotCrossSection(ml6, line={"row": 10}) + xc.plot_bc("LAK") + xc.plot_bc("SFR") + + ax = xc.ax + assert len(ax.collections) != 0, "Boundary condition was not drawn" + + for col in ax.collections: + assert isinstance( + col, PatchCollection + ), f"Unexpected collection type: {type(col)}" + + +@pytest.mark.mf6 +@pytest.mark.xfail( + reason="sometimes get LineCollections instead of PatchCollections" +) +def test_cross_section_bc_2models_mvr(example_data_path): + mpath = example_data_path / "mf6" / "test006_2models_mvr" + sim = MFSimulation.load(sim_ws=mpath) + ml6 = sim.get_model("parent") + xc = flopy.plot.PlotCrossSection(ml6, line={"column": 1}) + xc.plot_bc("MAW") + + ax = xc.ax + assert len(ax.collections) > 0, "Boundary condition was not drawn" + + for col in ax.collections: + assert isinstance( + col, PatchCollection + ), f"Unexpected collection type: {type(col)}" + + +@pytest.mark.mf6 +@pytest.mark.xfail( + reason="sometimes get LineCollections instead of PatchCollections" +) +def test_cross_section_bc_UZF_3lay(example_data_path): + mpath = example_data_path / "mf6" / "test001e_UZF_3lay" + sim = MFSimulation.load(sim_ws=mpath) + ml6 = sim.get_model("gwf_1") + + xc = flopy.plot.PlotCrossSection(ml6, line={"row": 0}) + xc.plot_bc("UZF") + + ax = xc.ax + assert len(ax.collections) != 0, "Boundary condition was not drawn" + + for col in ax.collections: + assert isinstance( + col, PatchCollection + ), f"Unexpected collection type: {type(col)}" diff --git a/autotest/test_plot_map_view.py b/autotest/test_plot_map_view.py new file mode 100644 index 0000000000..7c0aa8470b --- /dev/null +++ b/autotest/test_plot_map_view.py @@ -0,0 +1,200 @@ +import os + +import numpy as np +import pandas as pd +import pytest +from flaky import flaky +from matplotlib import pyplot as plt +from matplotlib import rcParams +from matplotlib.collections import ( + LineCollection, + PatchCollection, + PathCollection, + QuadMesh, +) +from modflow_devtools.markers import requires_exe, requires_pkg + +import flopy +from flopy.discretization import StructuredGrid +from flopy.mf6 import MFSimulation +from flopy.modflow import ( + Modflow, + ModflowBas, + ModflowDis, + ModflowLpf, + ModflowOc, + ModflowPcg, + ModflowWel, +) +from flopy.modpath import Modpath6, Modpath6Bas +from flopy.plot import PlotCrossSection, PlotMapView +from flopy.utils import CellBudgetFile, EndpointFile, HeadFile, PathlineFile + + +@requires_pkg("shapely") +def test_map_view(): + m = flopy.modflow.Modflow(rotation=20.0) + dis = flopy.modflow.ModflowDis( + m, nlay=1, nrow=40, ncol=20, delr=250.0, delc=250.0, top=10, botm=0 + ) + # transformation assigned by arguments + xll, yll, rotation = 500000.0, 2934000.0, 45.0 + + def check_vertices(): + xllp, yllp = pc._paths[0].vertices[0] + assert np.abs(xllp - xll) < 1e-6 + assert np.abs(yllp - yll) < 1e-6 + + m.modelgrid.set_coord_info(xoff=xll, yoff=yll, angrot=rotation) + modelmap = flopy.plot.PlotMapView(model=m) + pc = modelmap.plot_grid() + check_vertices() + + modelmap = flopy.plot.PlotMapView(modelgrid=m.modelgrid) + pc = modelmap.plot_grid() + check_vertices() + + mf = flopy.modflow.Modflow() + + # Model domain and grid definition + dis = flopy.modflow.ModflowDis( + mf, + nlay=1, + nrow=10, + ncol=20, + delr=1.0, + delc=1.0, + ) + xul, yul = 100.0, 210.0 + mg = mf.modelgrid + mf.modelgrid.set_coord_info( + xoff=mg._xul_to_xll(xul, 0.0), yoff=mg._yul_to_yll(yul, 0.0) + ) + verts = [[101.0, 201.0], [119.0, 209.0]] + modelxsect = flopy.plot.PlotCrossSection(model=mf, line={"line": verts}) + patchcollection = modelxsect.plot_grid() + + +@pytest.mark.mf6 +@pytest.mark.xfail(reason="sometimes get wrong collection type") +def test_map_view_bc_gwfs_disv(example_data_path): + mpath = example_data_path / "mf6" / "test003_gwfs_disv" + sim = MFSimulation.load(sim_ws=mpath) + 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: + assert isinstance( + col, (QuadMesh, PathCollection) + ), f"Unexpected collection type: {type(col)}" + + +@pytest.mark.mf6 +@pytest.mark.xfail(reason="sometimes get wrong collection type") +def test_map_view_bc_lake2tr(example_data_path): + mpath = example_data_path / "mf6" / "test045_lake2tr" + sim = MFSimulation.load(sim_ws=mpath) + 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: + assert isinstance( + col, (QuadMesh, PathCollection) + ), f"Unexpected collection type: {type(col)}" + + +@pytest.mark.mf6 +@pytest.mark.xfail(reason="sometimes get wrong collection type") +def test_map_view_bc_2models_mvr(example_data_path): + mpath = example_data_path / "mf6" / "test006_2models_mvr" + sim = MFSimulation.load(sim_ws=mpath) + 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") + ax = mapview.ax + + assert len(ax.collections) > 0, "Boundary condition was not drawn" + + mapview2 = flopy.plot.PlotMapView(model=ml6c, ax=mapview.ax) + mapview2.plot_bc("MAW") + ax = mapview2.ax + + assert len(ax.collections) > 0, "Boundary condition was not drawn" + + for col in ax.collections: + assert isinstance( + col, (QuadMesh, PathCollection) + ), f"Unexpected collection type: {type(col)}" + + +@pytest.mark.mf6 +@pytest.mark.xfail(reason="sometimes get wrong collection type") +def test_map_view_bc_UZF_3lay(example_data_path): + mpath = example_data_path / "mf6" / "test001e_UZF_3lay" + sim = MFSimulation.load(sim_ws=mpath) + ml6 = sim.get_model("gwf_1") + + mapview = flopy.plot.PlotMapView(model=ml6) + mapview.plot_bc("UZF") + ax = mapview.ax + + if len(ax.collections) == 0: + raise AssertionError("Boundary condition was not drawn") + + for col in ax.collections: + assert isinstance( + col, (QuadMesh, PathCollection) + ), f"Unexpected collection type: {type(col)}" + + +def test_map_view_contour(function_tmpdir): + arr = np.random.rand(10, 10) * 100 + arr[-1, :] = np.nan + delc = np.array([10] * 10, dtype=float) + delr = np.array([8] * 10, dtype=float) + top = np.ones((10, 10), dtype=float) + botm = np.ones((3, 10, 10), dtype=float) + botm[0] = 0.75 + botm[1] = 0.5 + botm[2] = 0.25 + idomain = np.ones((3, 10, 10)) + idomain[0, 0, :] = 0 + vmin = np.nanmin(arr) + vmax = np.nanmax(arr) + levels = np.linspace(vmin, vmax, 7) + + grid = StructuredGrid( + delc=delc, + delr=delr, + top=top, + botm=botm, + idomain=idomain, + lenuni=1, + nlay=3, + nrow=10, + ncol=10, + ) + + pmv = PlotMapView(modelgrid=grid, layer=0) + contours = pmv.contour_array(a=arr) + plt.savefig(function_tmpdir / "map_view_contour.png") + + # if we ever revert from standard contours to tricontours, restore this nan check + # for ix, lev in enumerate(contours.levels): + # if not np.allclose(lev, levels[ix]): + # raise AssertionError("TriContour NaN catch Failed") diff --git a/autotest/test_plot_pathlines.py b/autotest/test_plot_pathlines.py new file mode 100644 index 0000000000..6273db68fd --- /dev/null +++ b/autotest/test_plot_pathlines.py @@ -0,0 +1,379 @@ +import os + +import numpy as np +import pandas as pd +import pytest +from flaky import flaky +from matplotlib import pyplot as plt +from matplotlib import rcParams +from matplotlib.collections import ( + LineCollection, + PatchCollection, + PathCollection, + QuadMesh, +) +from modflow_devtools.markers import requires_exe, requires_pkg + +import flopy +from flopy.discretization import StructuredGrid +from flopy.mf6 import MFSimulation +from flopy.modflow import ( + Modflow, + ModflowBas, + ModflowDis, + ModflowLpf, + ModflowOc, + ModflowPcg, + ModflowWel, +) +from flopy.modpath import Modpath6, Modpath6Bas +from flopy.plot import PlotCrossSection, PlotMapView +from flopy.utils import CellBudgetFile, EndpointFile, HeadFile, PathlineFile + +# MP6 + + +@pytest.fixture +def mp6_sim(function_tmpdir, example_data_path): + # test with multi-layer example + load_ws = example_data_path / "mp6" + + ml = Modflow.load("EXAMPLE.nam", model_ws=load_ws, exe_name="mf2005") + ml.change_model_ws(function_tmpdir) + ml.write_input() + ml.run_model() + + mp = Modpath6( + modelname="ex6", + exe_name="mp6", + modflowmodel=ml, + model_ws=function_tmpdir, + ) + + mpb = Modpath6Bas( + mp, hdry=ml.lpf.hdry, laytyp=ml.lpf.laytyp, ibound=1, prsity=0.1 + ) + + sim = mp.create_mpsim( + trackdir="forward", + simtype="pathline", + packages="RCH", + start_time=(2, 0, 1.0), + ) + return ml, mp, sim + + +@requires_exe("mf2005", "mp6") +def test_plot_map_view_mp6_pathline(mp6_sim): + ml, mp, sim = mp6_sim + mp.write_input() + mp.run_model(silent=False) + + pthobj = PathlineFile(os.path.join(mp.model_ws, "ex6.mppth")) + well_pathlines = pthobj.get_destination_pathline_data( + dest_cells=[(4, 12, 12)] + ) + + def test_plot(pl): + mx = PlotMapView(model=ml) + mx.plot_grid() + mx.plot_bc("WEL", kper=2, color="blue") + pth = mx.plot_pathline(pl, colors="red") + # plt.show() + assert isinstance(pth, LineCollection) + assert len(pth._paths) == 114 + + # support pathlines as list of recarrays + test_plot(well_pathlines) + + # support pathlines as list of dataframes + test_plot([pd.DataFrame(pl) for pl in well_pathlines]) + + # support pathlines as single recarray + test_plot(np.concatenate(well_pathlines)) + + # support pathlines as single dataframe + test_plot(pd.DataFrame(np.concatenate(well_pathlines))) + + +@requires_exe("mf2005", "mp6") +def test_plot_cross_section_mp6_pathline(mp6_sim): + ml, mp, sim = mp6_sim + mp.write_input() + mp.run_model(silent=False) + + pthobj = PathlineFile(os.path.join(mp.model_ws, "ex6.mppth")) + well_pathlines = pthobj.get_destination_pathline_data( + dest_cells=[(4, 12, 12)] + ) + + def test_plot(pl): + mx = PlotCrossSection(model=ml, line={"row": 4}) + mx.plot_bc("WEL", kper=2, color="blue") + pth = mx.plot_pathline(pl, method="cell", colors="red") + assert isinstance(pth, LineCollection) + assert len(pth._paths) == 6 + + # support pathlines as list of recarrays + test_plot(well_pathlines) + + # support pathlines as list of dataframes + test_plot([pd.DataFrame(pl) for pl in well_pathlines]) + + # support pathlines as single recarray + test_plot(np.concatenate(well_pathlines)) + + # support pathlines as single dataframe + test_plot(pd.DataFrame(np.concatenate(well_pathlines))) + + +@requires_exe("mf2005", "mp6") +def test_plot_map_view_mp6_endpoint(mp6_sim): + ml, mp, sim = mp6_sim + mp.write_input() + mp.run_model(silent=False) + + pthobj = EndpointFile(os.path.join(mp.model_ws, "ex6.mpend")) + endpts = pthobj.get_alldata() + + # support endpoints as recarray + assert isinstance(endpts, np.recarray) + mv = PlotMapView(model=ml) + mv.plot_bc("WEL", kper=2, color="blue") + ep = mv.plot_endpoint(endpts, direction="ending") + # plt.show() + assert isinstance(ep, PathCollection) + + # support endpoints as dataframe + mv = PlotMapView(model=ml) + mv.plot_bc("WEL", kper=2, color="blue") + ep = mv.plot_endpoint(pd.DataFrame(endpts), direction="ending") + # plt.show() + assert isinstance(ep, PathCollection) + + # test various possibilities for endpoint color configuration. + # first, color kwarg as scalar + mv = PlotMapView(model=ml) + mv.plot_bc("WEL", kper=2, color="blue") + ep = mv.plot_endpoint(endpts, direction="ending", color="red") + # plt.show() + assert isinstance(ep, PathCollection) + + # c kwarg as array + mv = PlotMapView(model=ml) + mv.plot_bc("WEL", kper=2, color="blue") + ep = mv.plot_endpoint( + endpts, + direction="ending", + c=np.random.rand(625) * -1000, + cmap="viridis", + ) + # plt.show() + assert isinstance(ep, PathCollection) + + # colorbar: color by time to termination + mv = PlotMapView(model=ml) + mv.plot_bc("WEL", kper=2, color="blue") + ep = mv.plot_endpoint( + endpts, direction="ending", shrink=0.5, colorbar=True + ) + # plt.show() + assert isinstance(ep, PathCollection) + + # if both color and c are provided, c takes precedence + mv = PlotMapView(model=ml) + mv.plot_bc("WEL", kper=2, color="blue") + ep = mv.plot_endpoint( + endpts, + direction="ending", + color="red", + c=np.random.rand(625) * -1000, + cmap="viridis", + ) + # plt.show() + assert isinstance(ep, PathCollection) + + +# MP7 + + +simname = "test_plot" +nlay = 1 +nrow = 10 +ncol = 10 +top = 1.0 +botm = [0.0] +nper = 1 +perlen = 1.0 +nstp = 1 +tsmult = 1.0 +porosity = 0.1 + + +@pytest.fixture +def mf6_gwf_sim(module_tmpdir): + gwfname = f"{simname}_gwf" + + # create simulation + sim = flopy.mf6.MFSimulation( + sim_name=simname, + exe_name="mf6", + version="mf6", + sim_ws=module_tmpdir, + ) + + # create tdis package + flopy.mf6.modflow.mftdis.ModflowTdis( + sim, + pname="tdis", + time_units="DAYS", + nper=nper, + perioddata=[(perlen, nstp, tsmult)], + ) + + # create gwf model + gwf = flopy.mf6.ModflowGwf(sim, modelname=gwfname, save_flows=True) + + # create gwf discretization + flopy.mf6.modflow.mfgwfdis.ModflowGwfdis( + gwf, + pname="dis", + nlay=nlay, + nrow=nrow, + ncol=ncol, + ) + + # create gwf initial conditions package + flopy.mf6.modflow.mfgwfic.ModflowGwfic(gwf, pname="ic") + + # create gwf node property flow package + flopy.mf6.modflow.mfgwfnpf.ModflowGwfnpf( + gwf, + pname="npf", + save_saturation=True, + save_specific_discharge=True, + ) + + # create gwf chd package + spd = { + 0: [[(0, 0, 0), 1.0, 1.0], [(0, 9, 9), 0.0, 0.0]], + 1: [[(0, 0, 0), 0.0, 0.0], [(0, 9, 9), 1.0, 2.0]], + } + chd = flopy.mf6.ModflowGwfchd( + gwf, + pname="CHD-1", + stress_period_data=spd, + auxiliary=["concentration"], + ) + + # create gwf output control package + gwf_budget_file = f"{gwfname}.bud" + gwf_head_file = f"{gwfname}.hds" + oc = flopy.mf6.ModflowGwfoc( + gwf, + budget_filerecord=gwf_budget_file, + head_filerecord=gwf_head_file, + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + ) + + # create iterative model solution for gwf model + ims = flopy.mf6.ModflowIms(sim) + + return sim + + +@pytest.fixture +def mp7_sim(function_tmpdir, mf6_gwf_sim): + pass + + +@pytest.mark.skip(reason="todo") +@requires_exe("mf6", "mp7") +def test_plot_map_view_mp7_pathline(mp7_sim): + pass + + +@pytest.mark.skip(reason="todo") +@requires_exe("mf6", "mp7") +def test_plot_cross_section_mp7_pathline(mp7_sim): + pass + + +# MF6 PRT + + +# @pytest.fixture +# def mf6_prt_sim(function_tmpdir, mf6_gwf_sim): +# prtname = f"{simname}_prt" +# +# # create prt model +# prt = flopy.mf6.ModflowPrt(mf6_gwf_sim, modelname=prtname) +# +# # create prt discretization +# flopy.mf6.modflow.mfgwfdis.ModflowGwfdis( +# prt, +# pname="dis", +# nlay=nlay, +# nrow=nrow, +# ncol=ncol, +# ) +# +# # create mip package +# flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=porosity) +# +# # create prp package +# flopy.mf6.ModflowPrtprp( +# prt, +# pname="prp1", +# filename=f"{prtname}_1.prp", +# nreleasepts=len(releasepts), +# packagedata=releasepts, +# perioddata={0: ["FIRST"]}, +# ) +# +# # create output control package +# flopy.mf6.ModflowPrtoc( +# prt, +# pname="oc", +# track_filerecord=[prt_track_file], +# trackcsv_filerecord=[prt_track_csv_file], +# ) +# +# # create a flow model interface +# # todo Fienen's report (crash when FMI created but not needed) +# # flopy.mf6.ModflowPrtfmi( +# # prt, +# # packagedata=[ +# # ("GWFHEAD", gwf_head_file), +# # ("GWFBUDGET", gwf_budget_file), +# # ], +# # ) +# +# # create exchange +# flopy.mf6.ModflowGwfprt( +# sim, +# exgtype="GWF6-PRT6", +# exgmnamea=gwfname, +# exgmnameb=prtname, +# filename=f"{gwfname}.gwfprt", +# ) +# +# # add explicit model solution +# ems = flopy.mf6.ModflowEms( +# sim, +# pname="ems", +# filename=f"{prtname}.ems", +# ) +# sim.register_solution_package(ems, [prt.name]) +# +# return sim +# +# +# @requires_exe("mf6") +# def test_plot_map_view_prt_pathline(mf6_prt_sim): +# pass +# +# +# @requires_exe("mf6") +# def test_plot_cross_section_prt_pathline(mf6_prt_sim): +# pass diff --git a/autotest/test_plot_quasi3d.py b/autotest/test_plot_quasi3d.py new file mode 100644 index 0000000000..f26eb03dac --- /dev/null +++ b/autotest/test_plot_quasi3d.py @@ -0,0 +1,154 @@ +import os + +import numpy as np +import pandas as pd +import pytest +from flaky import flaky +from matplotlib import pyplot as plt +from matplotlib import rcParams +from matplotlib.collections import ( + LineCollection, + PatchCollection, + PathCollection, + QuadMesh, +) +from modflow_devtools.markers import requires_exe, requires_pkg + +import flopy +from flopy.discretization import StructuredGrid +from flopy.mf6 import MFSimulation +from flopy.modflow import ( + Modflow, + ModflowBas, + ModflowDis, + ModflowLpf, + ModflowOc, + ModflowPcg, + ModflowWel, +) +from flopy.modpath import Modpath6, Modpath6Bas +from flopy.plot import PlotCrossSection, PlotMapView +from flopy.utils import CellBudgetFile, EndpointFile, HeadFile, PathlineFile + + +@pytest.fixture +def quasi3d_model(function_tmpdir): + mf = Modflow("model_mf", model_ws=function_tmpdir, exe_name="mf2005") + + # Model domain and grid definition + Lx = 1000.0 + Ly = 1000.0 + ztop = 0.0 + zbot = -30.0 + nlay = 3 + nrow = 10 + ncol = 10 + delr = Lx / ncol + delc = Ly / nrow + laycbd = [0] * (nlay) + laycbd[0] = 1 + botm = np.linspace(ztop, zbot, nlay + np.sum(laycbd) + 1)[1:] + + # Create the discretization object + ModflowDis( + mf, + nlay, + nrow, + ncol, + delr=delr, + delc=delc, + top=ztop, + botm=botm, + laycbd=laycbd, + ) + + # Variables for the BAS package + ibound = np.ones((nlay, nrow, ncol), dtype=np.int32) + ibound[:, :, 0] = -1 + ibound[:, :, -1] = -1 + strt = np.ones((nlay, nrow, ncol), dtype=np.float32) + strt[:, :, 0] = 10.0 + strt[:, :, -1] = 0.0 + ModflowBas(mf, ibound=ibound, strt=strt) + + # Add LPF package to the MODFLOW model + ModflowLpf(mf, hk=10.0, vka=10.0, ipakcb=53, vkcb=10) + + # add a well + row = int((nrow - 1) / 2) + col = int((ncol - 1) / 2) + spd = {0: [[1, row, col, -1000]]} + ModflowWel(mf, stress_period_data=spd) + + # Add OC package to the MODFLOW model + spd = {(0, 0): ["save head", "save budget"]} + ModflowOc(mf, stress_period_data=spd, compact=True) + + # Add PCG package to the MODFLOW model + ModflowPcg(mf) + + # Write the MODFLOW model input files + mf.write_input() + + # Run the MODFLOW model + success, buff = mf.run_model() + + assert success, "test_plotting_with_quasi3d_layers() failed" + + return mf + + +@requires_exe("mf2005") +def test_map_plot_with_quasi3d_layers(quasi3d_model): + # read output + hf = HeadFile( + os.path.join(quasi3d_model.model_ws, f"{quasi3d_model.name}.hds") + ) + head = hf.get_data(totim=1.0) + cbb = CellBudgetFile( + os.path.join(quasi3d_model.model_ws, f"{quasi3d_model.name}.cbc") + ) + frf = cbb.get_data(text="FLOW RIGHT FACE", totim=1.0)[0] + fff = cbb.get_data(text="FLOW FRONT FACE", totim=1.0)[0] + flf = cbb.get_data(text="FLOW LOWER FACE", totim=1.0)[0] + + # plot a map + plt.figure() + mv = PlotMapView(model=quasi3d_model, layer=1) + mv.plot_grid() + mv.plot_array(head) + mv.contour_array(head) + mv.plot_ibound() + mv.plot_bc("wel") + mv.plot_vector(frf, fff) + plt.savefig(os.path.join(quasi3d_model.model_ws, "plt01.png")) + + +@requires_exe("mf2005") +def test_cross_section_with_quasi3d_layers(quasi3d_model): + # read output + hf = HeadFile( + os.path.join(quasi3d_model.model_ws, f"{quasi3d_model.name}.hds") + ) + head = hf.get_data(totim=1.0) + cbb = CellBudgetFile( + os.path.join(quasi3d_model.model_ws, f"{quasi3d_model.name}.cbc") + ) + frf = cbb.get_data(text="FLOW RIGHT FACE", totim=1.0)[0] + fff = cbb.get_data(text="FLOW FRONT FACE", totim=1.0)[0] + flf = cbb.get_data(text="FLOW LOWER FACE", totim=1.0)[0] + + # plot a cross-section + plt.figure() + cs = PlotCrossSection( + model=quasi3d_model, + line={"row": int((quasi3d_model.modelgrid.nrow - 1) / 2)}, + ) + cs.plot_grid() + cs.plot_array(head) + cs.contour_array(head) + cs.plot_ibound() + cs.plot_bc("wel") + cs.plot_vector(frf, fff, flf, head=head) + plt.savefig(os.path.join(quasi3d_model.model_ws, "plt02.png")) + plt.close() diff --git a/autotest/test_plotutil.py b/autotest/test_plotutil.py new file mode 100644 index 0000000000..997a4f865e --- /dev/null +++ b/autotest/test_plotutil.py @@ -0,0 +1,272 @@ +import numpy as np +import pandas as pd +from modflow_devtools.markers import requires_exe + +from flopy.plot.plotutil import to_mp7_endpoints, to_mp7_pathlines + +# test PRT-MP7 pathline conversion functions + + +mp7_pl_cols = [ + "particleid", + "particlegroup", + "sequencenumber", + "particleidloc", + "time", + "xloc", + "yloc", + "zloc", + "x", + "y", + "z", + "node", + "k", + "stressperiod", + "timestep", +] + + +mp7_ep_cols = [ + "particleid", + "particlegroup", + "particleidloc", + "time", + "time0", + "xloc", + "xloc0", + "yloc", + "yloc0", + "zloc", + "zloc0", + "x0", + "y0", + "z0", + "x", + "y", + "z", + "node", + "node0", + "k", + "k0", + "zone", + "zone0", + "initialcellface", + "cellface", + "status", +] + + +pls = pd.DataFrame.from_records( + [ + [ + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 0, + 0, + 0, + 0.0, + 0.000000, + 0.100000, + 9.1, + 0.5, + "PRP000000001", + ], + [ + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 0, + 1, + 1, + 0.0, + 0.063460, + 0.111111, + 9.0, + 0.5, + "PRP000000001", + ], + [ + 1, + 1, + 1, + 1, + 1, + 1, + 11, + 0, + 1, + 1, + 0.0, + 0.830431, + 0.184020, + 8.0, + 0.5, + "PRP000000001", + ], + [ + 1, + 1, + 1, + 1, + 1, + 1, + 21, + 0, + 1, + 1, + 0.0, + 2.026390, + 0.267596, + 7.0, + 0.5, + "PRP000000001", + ], + [ + 1, + 1, + 1, + 1, + 1, + 1, + 31, + 0, + 1, + 1, + 0.0, + 3.704265, + 0.360604, + 6.0, + 0.5, + "PRP000000001", + ], + [ + 1, + 1, + 1, + 1, + 1, + 1, + 60, + 0, + 1, + 1, + 0.0, + 39.087992, + 9.639587, + 4.0, + 0.5, + "PRP000000001", + ], + [ + 1, + 1, + 1, + 1, + 1, + 1, + 70, + 0, + 1, + 1, + 0.0, + 40.765791, + 9.732597, + 3.0, + 0.5, + "PRP000000001", + ], + [ + 1, + 1, + 1, + 1, + 1, + 1, + 80, + 0, + 1, + 1, + 0.0, + 41.961755, + 9.816110, + 2.0, + 0.5, + "PRP000000001", + ], + [ + 1, + 1, + 1, + 1, + 1, + 1, + 90, + 0, + 1, + 1, + 0.0, + 42.728752, + 9.888968, + 1.0, + 0.5, + "PRP000000001", + ], + [ + 1, + 1, + 1, + 1, + 1, + 1, + 100, + 0, + 5, + 3, + 0.0, + 42.728752, + 9.888968, + 1.0, + 0.5, + "PRP000000001", + ], + ], + columns=[ + "kper", + "kstp", + "imdl", + "iprp", + "irpt", + "ilay", + "icell", + "izone", + "istatus", + "ireason", + "trelease", + "t", + "x", + "y", + "z", + "name", + ], +) + + +def test_to_mp7_pathlines(): + mp7_pls = to_mp7_pathlines(pls) + assert len(mp7_pls) == 10 + assert set(dict(mp7_pls.dtypes).keys()) == set(mp7_pl_cols) + + +def test_to_mp7_endpoints(): + mp7_eps = to_mp7_endpoints(pls) + assert len(mp7_eps) == 1 + assert set(dict(mp7_eps.dtypes).keys()) == set(mp7_ep_cols) diff --git a/flopy/modpath/mp7particledata.py b/flopy/modpath/mp7particledata.py index 9f8b8b393d..7217c1fc8e 100644 --- a/flopy/modpath/mp7particledata.py +++ b/flopy/modpath/mp7particledata.py @@ -4,8 +4,11 @@ """ +from itertools import product +from typing import Generator, Iterator, Tuple import numpy as np +import pandas as pd from numpy.lib.recfunctions import unstructured_to_structured from ..utils.recarray_utils import create_empty_recarray @@ -314,21 +317,17 @@ def __init__( self.particlecount = particledata.shape[0] self.particleidoption = particleidoption self.locationstyle = locationstyle - self.particledata = particledata - - return + self.dtype = particledata.dtype + self.particledata = pd.DataFrame.from_records(particledata) def write(self, f=None): """ + Write the particle data to disk. Parameters ---------- f : fileobject Fileobject that is open with write access - - Returns - ------- - """ # validate that a valid file object was passed if not hasattr(f, "write"): @@ -338,7 +337,7 @@ def write(self, f=None): ) # particle data item 4 and 5 - d = np.recarray.copy(self.particledata) + d = np.recarray.copy(self.particledata.to_records(index=False)) lnames = [name.lower() for name in d.dtype.names] # Add one to the kij and node indices for idx in ( @@ -358,7 +357,101 @@ def write(self, f=None): for v in d: f.write(fmt.format(*v)) - return + def to_coords(self, grid): + """ + Convert particle data to global coordinates. + + Parameters + ---------- + grid : The grid to use for locating particle coordinates + + Returns + ------- + A list of coordinate tuples (x, y, z) + """ + + def cvt_xy(p, vs): + mn, mx = min(vs), max(vs) + span = mx - mn + return mn + span * p + + if grid.grid_type == "structured": + if not hasattr(self.particledata, "k"): + raise ValueError( + f"Particle representation is not structured but grid is" + ) + + def cvt_z(p, k, i, j): + mn, mx = grid.botm[k, i, j], grid.top[i, j] + span = mx - mn + return mn + span * p + + def convert(row) -> Tuple[float, float, float]: + verts = grid.get_cell_vertices(row.i, row.j) + xs, ys = list(zip(*verts)) + return ( + cvt_xy(row.localx, xs), + cvt_xy(row.localy, ys), + cvt_z(row.localz, row.k, row.i, row.j), + ) + + else: + if hasattr(self.particledata, "k"): + raise ValueError( + f"Particle representation is structured but grid is not" + ) + + def cvt_z(p, nn): + k, j = grid.get_lni([nn])[0] + mn, mx = grid.botm[k, j], grid.top[j] + span = mx - mn + return mn + span * p + + def convert(row) -> Tuple[float, float, float]: + verts = grid.get_cell_vertices(row.node) + xs, ys = list(zip(*verts)) + return ( + cvt_xy(row.localx, xs), + cvt_xy(row.localy, ys), + cvt_z(row.localz, row.node), + ) + + for row in self.particledata.itertuples(): + yield convert(row) + + def to_prp(self, grid): + """ + Convert particle data to PRT particle data format. + + Parameters + ---------- + grid : The grid to use for locating particle coordinates + + Returns + ------- + A list of PRT particle release locations, as defined + for particle release point (PRP) packagedata tuples: + release point id, cell id (k, i, j) if structured or + (k, j) if unstructured, x coordinate, y coordinate, + z coordinate. + """ + + for i, (p, c) in enumerate( + zip( + self.particledata.itertuples(index=False), + self.to_coords(grid), + ) + ): + yield ( + i, # release point index + # cell id for vertex grid + grid.get_lni([p.node])[0] if "node" in self.particledata + # cell id for structured grid + else (p.k, p.i, p.j), + c[0], # x + c[1], # y + c[2], # z + ) def _get_dtype(self, structured, particleid): """ @@ -421,7 +514,7 @@ def _fmt_string(self): """ fmts = [] - for field in self.particledata.dtype.descr: + for field in self.dtype.descr: vtype = field[1][1].lower() if vtype == "i" or vtype == "b": fmts.append("{:9d}") @@ -542,7 +635,6 @@ def __init__( self.columndivisions5 = columndivisions5 self.rowdivisions6 = rowdivisions6 self.columndivisions6 = columndivisions6 - return def write(self, f=None): """ @@ -637,7 +729,6 @@ def __init__( self.columncelldivisions = columncelldivisions self.rowcelldivisions = rowcelldivisions self.layercelldivisions = layercelldivisions - return def write(self, f=None): """ @@ -668,7 +759,177 @@ def write(self, f=None): ) f.write(line) - return + +def get_release_points(subdivisiondata, grid, k=None, i=None, j=None, nn=None): + if nn is None and (k is None or i is None or j is None): + raise ValueError( + f"A cell (node) must be specified by indices (for structured grids) or node number (for vertex/unstructured)" + ) + + rpts = [] + template = [k, i, j] if nn is None else [nn] + + # get cell coords and span in each dimension + if not (k is None or i is None or j is None): + verts = grid.get_cell_vertices(i, j) + minz, maxz = grid.botm[k, i, j], grid.top[i, j] + elif nn is not None: + verts = grid.get_cell_vertices(nn) + if grid.grid_type == "structured": + k, i, j = grid.get_lrc([nn])[0] + minz, maxz = ( + grid.botm[k, i, j], + grid.top[i, j] if k == 0 else grid.botm[k - 1, i, j], + ) + else: + k, j = grid.get_lni([nn])[0] + minz, maxz = ( + grid.botm[k, j], + grid.top[j] if k == 0 else grid.botm[k - 1, j], + ) + else: + raise ValueError( + f"A cell (node) must be specified by indices (for structured grids) or node number (for vertex/unstructured)" + ) + xs, ys = list(zip(*verts)) + minx, maxx = min(xs), max(xs) + miny, maxy = min(ys), max(ys) + xspan = maxx - minx + yspan = maxy - miny + zspan = maxz - minz + + if isinstance(subdivisiondata, FaceDataType): + # face perpendicular to x (outer) + if ( + subdivisiondata.verticaldivisions1 > 0 + and subdivisiondata.horizontaldivisions1 > 0 + ): + yincr = yspan / subdivisiondata.horizontaldivisions1 + ylocs = [ + (miny + (yincr * 0.5) + (yincr * d)) + for d in range(subdivisiondata.horizontaldivisions1) + ] + zincr = zspan / subdivisiondata.verticaldivisions1 + zlocs = [ + (minz + (zincr * 0.5) + (zincr * d)) + for d in range(subdivisiondata.verticaldivisions1) + ] + prod = list(product(*[ylocs, zlocs])) + rpts = rpts + [template + [minx, p[0], p[1]] for p in prod] + + # face perpendicular to x (outer) + if ( + subdivisiondata.verticaldivisions2 > 0 + and subdivisiondata.horizontaldivisions2 > 0 + ): + yincr = yspan / subdivisiondata.horizontaldivisions2 + ylocs = [ + (miny + (yincr * 0.5) + (yincr * d)) + for d in range(subdivisiondata.horizontaldivisions2) + ] + zincr = zspan / subdivisiondata.verticaldivisions2 + zlocs = [ + (minz + (zincr * 0.5) + (zincr * d)) + for d in range(subdivisiondata.verticaldivisions2) + ] + prod = list(product(*[ylocs, zlocs])) + rpts = rpts + [template + [maxx, p[0], p[1]] for p in prod] + + # face perpendicular to y (inner) + if ( + subdivisiondata.verticaldivisions3 > 0 + and subdivisiondata.horizontaldivisions3 > 0 + ): + xincr = xspan / subdivisiondata.horizontaldivisions3 + xlocs = [ + (minx + (xincr * 0.5) + (xincr * rd)) + for rd in range(subdivisiondata.horizontaldivisions3) + ] + zincr = zspan / subdivisiondata.verticaldivisions3 + zlocs = [ + (minz + (zincr * 0.5) + (zincr * d)) + for d in range(subdivisiondata.verticaldivisions3) + ] + prod = list(product(*[xlocs, zlocs])) + rpts = rpts + [template + [p[0], miny, p[1]] for p in prod] + + # face perpendicular to y (outer) + if ( + subdivisiondata.verticaldivisions4 > 0 + and subdivisiondata.horizontaldivisions4 > 0 + ): + xincr = xspan / subdivisiondata.horizontaldivisions4 + xlocs = [ + (minx + (xincr * 0.5) + (xincr * rd)) + for rd in range(subdivisiondata.horizontaldivisions4) + ] + zincr = zspan / subdivisiondata.verticaldivisions4 + zlocs = [ + (minz + (zincr * 0.5) + (zincr * d)) + for d in range(subdivisiondata.verticaldivisions4) + ] + prod = list(product(*[xlocs, zlocs])) + rpts = rpts + [template + [p[0], maxy, p[1]] for p in prod] + + # bottom face + if ( + subdivisiondata.rowdivisions5 > 0 + and subdivisiondata.columndivisions5 > 0 + ): + xincr = xspan / subdivisiondata.columndivisions5 + xlocs = [ + (minx + (xincr * 0.5) + (xincr * rd)) + for rd in range(subdivisiondata.columndivisions5) + ] + yincr = yspan / subdivisiondata.rowdivisions5 + ylocs = [ + (miny + (yincr * 0.5) + (yincr * rd)) + for rd in range(subdivisiondata.rowdivisions5) + ] + prod = list(product(*[xlocs, ylocs])) + rpts = rpts + [template + [p[0], p[1], minz] for p in prod] + + # top face + if ( + subdivisiondata.rowdivisions6 > 0 + and subdivisiondata.columndivisions6 > 0 + ): + xincr = xspan / subdivisiondata.columndivisions6 + xlocs = [ + (minx + (xincr * 0.5) + (xincr * rd)) + for rd in range(subdivisiondata.columndivisions6) + ] + yincr = yspan / subdivisiondata.rowdivisions5 + ylocs = [ + (miny + (yincr * 0.5) + (yincr * rd)) + for rd in range(subdivisiondata.rowdivisions6) + ] + prod = list(product(*[xlocs, ylocs])) + rpts = rpts + [template + [p[0], p[1], maxz] for p in prod] + elif isinstance(subdivisiondata, CellDataType): + xincr = xspan / subdivisiondata.columncelldivisions + xlocs = [ + (minx + (xincr * 0.5) + (xincr * rd)) + for rd in range(subdivisiondata.columncelldivisions) + ] + yincr = yspan / subdivisiondata.rowcelldivisions + ylocs = [ + (miny + (yincr * 0.5) + (yincr * d)) + for d in range(subdivisiondata.rowcelldivisions) + ] + zincr = zspan / subdivisiondata.layercelldivisions + zlocs = [ + (minz + (zincr * 0.5) + (zincr * d)) + for d in range(subdivisiondata.layercelldivisions) + ] + prod = list(product(*[xlocs, ylocs, zlocs])) + rpts = rpts + [template + [p[0], p[1], p[2]] for p in prod] + else: + raise ValueError( + f"Unsupported subdivision data type: {type(subdivisiondata)}" + ) + + return rpts class LRCParticleData: @@ -686,7 +947,7 @@ class LRCParticleData: CellDataTypes that are used to create one or more particle templates in a particle group. If subdivisiondata is None, a default CellDataType with 27 particles per cell will be created (default is None). - lrcregions : list of lists tuples or np.ndarrays + lrcregions : list of lists of tuples or np.ndarrays Layer, row, column (zero-based) regions with particles created using the specified template parameters. A region is defined as a list/tuple of minlayer, minrow, mincolumn, maxlayer, maxrow, maxcolumn values. @@ -699,7 +960,7 @@ class LRCParticleData: -------- >>> import flopy - >>> pg = flopy.modpath.LRCParticleData(lrcregions=[0, 0, 0, 3, 10, 10]) + >>> pg = flopy.modpath.LRCParticleData(lrcregions=[[0, 0, 0, 3, 10, 10]]) """ @@ -740,8 +1001,8 @@ def __init__(self, subdivisiondata=None, lrcregions=None): "a list with lists, tuples, or arrays".format(self.name) ) t = [] - for lrcregion in lrcregions: - t.append(np.array(lrcregion, dtype=np.int32)) + for region in lrcregions: + t.append(np.array(region, dtype=np.int32)) lrcregions = t else: raise TypeError( @@ -758,11 +1019,11 @@ def __init__(self, subdivisiondata=None, lrcregions=None): ) # validate that there are 6 columns in each lrcregions entry - for idx, lrcregion in enumerate(lrcregions): - shapel = lrcregion.shape + for idx, region in enumerate(lrcregions): + shapel = np.array(region).shape if len(shapel) == 1: - lrcregions[idx] = lrcregion.reshape(1, shapel) - shapel = lrcregion[idx].shape + region = region.reshape(1, shapel[0]) + shapel = region.shape if shapel[1] != 6: raise ValueError( "{}: Each lrcregions entry must " @@ -770,17 +1031,15 @@ def __init__(self, subdivisiondata=None, lrcregions=None): "{} columns".format(self.name, shapel[1]) ) - # totalcellregioncount = 0 - for lrcregion in lrcregions: - totalcellregioncount += lrcregion.shape[0] + for region in lrcregions: + totalcellregioncount += region.shape[0] # assign attributes self.particletemplatecount = shape self.totalcellregioncount = totalcellregioncount self.subdivisiondata = subdivisiondata self.lrcregions = lrcregions - return def write(self, f=None): """ @@ -824,6 +1083,47 @@ def write(self, f=None): return + def to_prp(self, grid) -> Iterator[tuple]: + """ + Convert particle data to PRT particle release point (PRP) + package data entries. + + Parameters + ---------- + grid : The grid to use for locating particle coordinates + + Returns + ------- + A list of particle release points as defined for PRT + particle release point (PRP) package data recarrays: + release point id, cell id (k, i, j), x coordinate, y + coordinate, z coordinate. + """ + + if grid.grid_type != "structured": + raise ValueError( + f"Particle representation is structured but grid is not" + ) + + for region in self.lrcregions: + mink, mini, minj, maxk, maxi, maxj = region + for k in range(mink, maxk + 1): + for i in range(mini, maxi + 1): + for j in range(minj, maxj + 1): + for sd in self.subdivisiondata: + rpts = get_release_points(sd, grid, k, i, j) + for irpt, rpt in enumerate(rpts): + assert rpt[0] == k + assert rpt[1] == i + assert rpt[2] == j + yield ( + irpt, + (k, i, j), + rpt[3], + rpt[4], + rpt[5], + ) + class NodeParticleData: """ @@ -939,7 +1239,6 @@ def __init__(self, subdivisiondata=None, nodes=None): self.totalcellcount = totalcellcount self.subdivisiondata = subdivisiondata self.nodedata = nodes - return def write(self, f=None): """ @@ -985,4 +1284,31 @@ def write(self, f=None): line += "\n" f.write(line) - return + def to_prp(self, grid) -> Iterator[tuple]: + """ + Convert particle data to PRT particle release point (PRP) + package data entries. + + Parameters + ---------- + grid : The grid to use for locating particle coordinates + + Returns + ------- + A list of particle release points as defined for PRT + particle release point (PRP) package data recarrays: + release point id, cell id (k, j), x coordinate, y + coordinate, z coordinate. + """ + + for sd in self.subdivisiondata: + for nd in self.nodedata: + rpts = get_release_points(sd, grid, nn=int(nd[0])) + for i, rpt in enumerate(rpts): + yield ( + i, # release point index + grid.get_lni([rpt[0]])[0], # get (k, j) from nn + rpt[1], # x + rpt[2], # y + rpt[3], # z + ) diff --git a/flopy/plot/crosssection.py b/flopy/plot/crosssection.py index 1e50ee69e8..8e870aa748 100644 --- a/flopy/plot/crosssection.py +++ b/flopy/plot/crosssection.py @@ -6,10 +6,12 @@ import numpy as np import pandas as pd from matplotlib.patches import Polygon +from numpy.lib.recfunctions import stack_arrays from ..utils import geometry, import_optional_dependency from ..utils.geospatial_utils import GeoSpatialUtil from . import plotutil +from .plotutil import to_mp7_endpoints, to_mp7_pathlines warnings.simplefilter("always", PendingDeprecationWarning) @@ -1054,15 +1056,27 @@ def plot_pathline( self, pl, travel_time=None, method="cell", head=None, **kwargs ): """ - Plot the MODPATH pathlines + Plot particle pathlines. Compatible with MODFLOW 6 PRT particle track + data format, or MODPATH 6 or 7 pathline data format. Parameters ---------- - pl : list of rec arrays or a single rec array - rec array or list of rec arrays is data returned from - modpathfile PathlineFile get_data() or get_alldata() - methods. Data in rec array is 'x', 'y', 'z', 'time', - 'k', and 'particleid'. + pl : list of recarrays or dataframes, or a single recarray or dataframe + Particle pathline data. If a list of recarrays or dataframes, + each must contain the path of only a single particle. If just + one recarray or dataframe, it should contain the paths of all + particles. The flopy.utils.modpathfile.PathlineFile.get_data() + or get_alldata() return value may be passed directly as this + argument. + + For MODPATH 6 or 7 pathlines, columns must include 'x', 'y', 'z', + 'time', 'k', and 'particleid'. Additional columns are ignored. + + For MODFLOW 6 PRT pathlines, columns must include 'x', 'y', 'z', + 't', 'trelease', 'imdl', 'iprp', 'irpt', and 'ilay'. Additional + columns are ignored. Note that MODFLOW 6 PRT does not assign to + particles a unique ID, but infers particle identity from 'imdl', + 'iprp', 'irpt', and 'trelease' combos (i.e. via composite key). travel_time : float or str travel_time is a travel time selection for the displayed pathlines. If a float is passed then pathlines with times @@ -1086,32 +1100,67 @@ def plot_pathline( Returns ------- lc : matplotlib.collections.LineCollection - + The pathlines added to the plot. """ from matplotlib.collections import LineCollection - # make sure pathlines is a list + # make sure pl is a list if not isinstance(pl, list): - pids = np.unique(pl["particleid"]) - if len(pids) > 1: - pl = [pl[pl["particleid"] == pid] for pid in pids] - else: - pl = [pl] + if not isinstance(pl, (np.ndarray, pd.DataFrame)): + raise TypeError( + "Pathline data must be a list of recarrays or dataframes, " + f"or a single recarray or dataframe, got {type(pl)}" + ) + pl = [pl] - # make sure each element in pl is a recarray + # convert dataframes to recarrays if needed pl = [ p.to_records(index=False) if isinstance(p, pd.DataFrame) else p for p in pl ] + def convert(p): + mp7_names = ["x", "y", "z", "time", "k", "particleid"] + prt_names = [ + "x", + "y", + "z", + "t", + "trelease", + "imdl", + "iprp", + "irpt", + "ilay", + ] + + # check format + if not ( + all(n in p.dtype.names for n in mp7_names) + or all(n in p.dtype.names for n in prt_names) + ): + raise ValueError( + "Pathline data must contain all of the following columns: " + f"{mp7_names} for MODPATH 6-7, or {prt_names} for MODFLOW 6 PRT" + ) + + return to_mp7_pathlines(p) if "t" in p.dtype.names else p + + # convert prt to mp7 format if needed + pl = [convert(p) for p in pl] + + # merge pathlines then split on particleid + pls = stack_arrays(pl, asrecarray=True, usemask=False) + pids = np.unique(pls["particleid"]) + pl = [pls[pls["particleid"] == pid] for pid in pids] + + # configure plot settings 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" @@ -1176,7 +1225,7 @@ def plot_timeseries( self, ts, travel_time=None, method="cell", head=None, **kwargs ): """ - Plot the MODPATH timeseries. + Plot the MODPATH timeseries. Not compatible with MODFLOW 6 PRT. Parameters ---------- @@ -1221,22 +1270,64 @@ def plot_endpoint( **kwargs, ): """ + Plot particle endpoints. Compatible with MODFLOW 6 PRT particle + track data format, or MODPATH 6 or 7 endpoint data format. Parameters ---------- - + ep : recarray or dataframe + A numpy recarray with the endpoint particle data from the + MODPATH endpoint file. + + For MODFLOW 6 PRT pathlines, columns must include 'x', 'y', 'z', + 't', 'trelease', 'imdl', 'iprp', 'irpt', and 'ilay'. Additional + columns are ignored. Note that MODFLOW 6 PRT does not assign to + particles a unique ID, but infers particle identity from 'imdl', + 'iprp', 'irpt', and 'trelease' combos (i.e. via composite key). + direction : str + String defining if starting or ending particle locations should be + considered. (default is 'ending') + selection : tuple + tuple that defines the zero-base layer, row, column location + (l, r, c) to use to make a selection of particle endpoints. + The selection could be a well location to determine capture zone + for the well. If selection is None, all particle endpoints for + the user-sepcified direction will be plotted. (default is None) + selection_direction : str + String defining is a selection should be made on starting or + ending particle locations. If selection is not None and + selection_direction is None, the selection direction will be set + to the opposite of direction. (default is None) + + kwargs : ax, c, s or size, colorbar, colorbar_label, shrink. The + remaining kwargs are passed into the matplotlib scatter + method. If colorbar is True a colorbar will be added to the plot. + If colorbar_label is passed in and colorbar is True then + colorbar_label will be passed to the colorbar set_label() + method. If shrink is passed in and colorbar is True then + the colorbar size will be set using shrink. Returns ------- + sp : matplotlib.collections.PathCollection + The PathCollection added to the plot. """ - ax = kwargs.pop("ax", self.ax) - # colorbar kwargs + # convert ep to recarray if needed + if isinstance(ep, pd.DataFrame): + ep = ep.to_records(index=False) + + # convert ep from prt to mp7 format if needed + if "t" in ep.dtype.names: + from .plotutil import to_mp7_endpoints + + ep = to_mp7_endpoints(ep) + + # configure plot settings + ax = kwargs.pop("ax", self.ax) createcb = kwargs.pop("colorbar", False) colorbar_label = kwargs.pop("colorbar_label", "Endpoint Time") shrink = float(kwargs.pop("shrink", 1.0)) - - # marker kwargs s = kwargs.pop("s", np.sqrt(50)) s = float(kwargs.pop("size", s)) ** 2.0 diff --git a/flopy/plot/map.py b/flopy/plot/map.py index 98ecdd44d8..2bef039da7 100644 --- a/flopy/plot/map.py +++ b/flopy/plot/map.py @@ -6,9 +6,11 @@ import pandas as pd from matplotlib.collections import LineCollection, PathCollection from matplotlib.path import Path +from numpy.lib.recfunctions import stack_arrays from ..utils import geometry from . import plotutil +from .plotutil import to_mp7_endpoints, to_mp7_pathlines warnings.simplefilter("always", PendingDeprecationWarning) @@ -696,7 +698,8 @@ def plot_vector( def plot_pathline(self, pl, travel_time=None, **kwargs): """ - Plot MODPATH pathlines. + Plot particle pathlines. Compatible with MODFLOW 6 PRT particle track + data format, or MODPATH 6 or 7 pathline data format. Parameters ---------- @@ -704,11 +707,18 @@ def plot_pathline(self, pl, travel_time=None, **kwargs): Particle pathline data. If a list of recarrays or dataframes, each must contain the path of only a single particle. If just one recarray or dataframe, it should contain the paths of all - particles. Pathline data returned from PathlineFile.get_data() - or get_alldata() can be passed directly as this argument. Data - columns should be 'x', 'y', 'z', 'time', 'k', and 'particleid' - at minimum. Additional columns are ignored. The 'particleid' - column must be unique to each particle path. + particles. The flopy.utils.modpathfile.PathlineFile.get_data() + or get_alldata() return value may be passed directly as this + argument. + + For MODPATH 6 or 7 pathlines, columns must include 'x', 'y', 'z', + 'time', 'k', and 'particleid'. Additional columns are ignored. + + For MODFLOW 6 PRT pathlines, columns must include 'x', 'y', 'z', + 't', 'trelease', 'imdl', 'iprp', 'irpt', and 'ilay'. Additional + columns are ignored. Note that MODFLOW 6 PRT does not assign to + particles a unique ID, but infers particle identity from 'imdl', + 'iprp', 'irpt', and 'trelease' combos (i.e. via composite key). travel_time : float or str Travel time selection. If a float, then pathlines with total time less than or equal to the given value are plotted. If a @@ -729,19 +739,56 @@ def plot_pathline(self, pl, travel_time=None, **kwargs): from matplotlib.collections import LineCollection - # make sure pathlines is a list + # make sure pl is a list if not isinstance(pl, list): - pids = np.unique(pl["particleid"]) - if len(pids) > 1: - pl = [pl[pl["particleid"] == pid] for pid in pids] - else: - pl = [pl] + if not isinstance(pl, (np.ndarray, pd.DataFrame)): + raise TypeError( + "Pathline data must be a list of recarrays or dataframes, " + f"or a single recarray or dataframe, got {type(pl)}" + ) + pl = [pl] + # convert dataframes to recarrays if needed pl = [ p.to_records(index=False) if isinstance(p, pd.DataFrame) else p for p in pl ] + def convert(p): + mp7_names = ["x", "y", "z", "time", "k", "particleid"] + prt_names = [ + "x", + "y", + "z", + "t", + "trelease", + "imdl", + "iprp", + "irpt", + "ilay", + ] + + # check format + if not ( + all(n in p.dtype.names for n in mp7_names) + or all(n in p.dtype.names for n in prt_names) + ): + raise ValueError( + "Pathline data must contain all of the following columns: " + f"{mp7_names} for MODPATH 6-7, or {prt_names} for MODFLOW 6 PRT" + ) + + return to_mp7_pathlines(p) if "t" in p.dtype.names else p + + # convert prt to mp7 format if needed + pl = [convert(p) for p in pl] + + # merge pathlines then split on particleid + pls = stack_arrays(pl, asrecarray=True, usemask=False) + pids = np.unique(pls["particleid"]) + pl = [pls[pls["particleid"] == pid] for pid in pids] + + # configure layer if "layer" in kwargs: kon = kwargs.pop("layer") if isinstance(kon, bytes): @@ -754,19 +801,21 @@ def plot_pathline(self, pl, travel_time=None, **kwargs): else: kon = self.layer + # configure plot settings 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" + # compose pathlines linecol = [] markers = [] for p in pl: + # filter by travel time tp = plotutil.filter_modpath_by_travel_time(p, travel_time) # transform data! @@ -777,8 +826,10 @@ def plot_pathline(self, pl, travel_time=None, **kwargs): self.mg.yoffset, self.mg.angrot_radians, ) + # build polyline array arr = np.vstack((x0r, y0r)).T + # select based on layer if kon >= 0: kk = p["k"].copy().reshape(p.shape[0], 1) @@ -786,7 +837,8 @@ def plot_pathline(self, pl, travel_time=None, **kwargs): arr = np.ma.masked_where((kk != kon), arr) else: arr = np.ma.asarray(arr) - # append line to linecol if there is some unmasked segment + + # append pathline if there are any unmasked segments if not arr.mask.all(): linecol.append(arr) if not arr.mask.all(): @@ -795,6 +847,7 @@ def plot_pathline(self, pl, travel_time=None, **kwargs): for xy in arr[::markerevery]: if not np.all(xy.mask): markers.append(xy) + # create line collection lc = None if len(linecol) > 0: @@ -811,13 +864,15 @@ def plot_pathline(self, pl, travel_time=None, **kwargs): ms=markersize, ) + # set axis limits ax.set_xlim(self.extent[0], self.extent[1]) ax.set_ylim(self.extent[2], self.extent[3]) + return lc def plot_timeseries(self, ts, travel_time=None, **kwargs): """ - Plot MODPATH timeseries. + Plot MODPATH 6 or 7 timeseries. Incompatible with MODFLOW 6 PRT. Parameters ---------- @@ -861,13 +916,20 @@ def plot_endpoint( **kwargs, ): """ - Plot MODPATH endpoints. + Plot particle endpoints. Compatible with MODFLOW 6 PRT particle + track data format, or MODPATH 6 or 7 endpoint data format. Parameters ---------- ep : recarray or dataframe A numpy recarray with the endpoint particle data from the - MODPATH endpoint file + MODPATH endpoint file. + + For MODFLOW 6 PRT pathlines, columns must include 'x', 'y', 'z', + 't', 'trelease', 'imdl', 'iprp', 'irpt', and 'ilay'. Additional + columns are ignored. Note that MODFLOW 6 PRT does not assign to + particles a unique ID, but infers particle identity from 'imdl', + 'iprp', 'irpt', and 'trelease' combos (i.e. via composite key). direction : str String defining if starting or ending particle locations should be considered. (default is 'ending') @@ -898,6 +960,17 @@ def plot_endpoint( """ + # convert ep to recarray if needed + if isinstance(ep, pd.DataFrame): + ep = ep.to_records(index=False) + + # convert ep from prt to mp7 format if needed + if "t" in ep.dtype.names: + from .plotutil import to_mp7_endpoints + + ep = to_mp7_endpoints(ep) + + # parse selection options ax = kwargs.pop("ax", self.ax) tep, _, xp, yp = plotutil.parse_modpath_selection_options( ep, direction, selection, selection_direction diff --git a/flopy/plot/plotutil.py b/flopy/plot/plotutil.py index 52001cbd87..192edcda73 100644 --- a/flopy/plot/plotutil.py +++ b/flopy/plot/plotutil.py @@ -10,6 +10,7 @@ import matplotlib.pyplot as plt import numpy as np +import pandas as pd from ..datbase import DataInterface, DataType from ..utils import Util3d, import_optional_dependency @@ -2591,3 +2592,218 @@ def parse_modpath_selection_options( tep = ep.copy() return tep, istart, xp, yp + + +def to_mp7_pathlines(data) -> pd.DataFrame: + """ + Convert MODFLOW 6 PRT pathline data to MODPATH 7 pathline format. + + Parameters + ---------- + data : np.recarray or pd.DataFrame + MODFLOW 6 PRT pathline data + + Returns + ------- + pd.DataFrame + """ + + # convert to dataframe if needed + if isinstance(data, np.recarray): + data = pd.DataFrame(data) + + # assign a unique particle index column incrementing an integer + # for each unique combination of irpt, iprp, imdl, and trelease + data = data.sort_values(["imdl", "iprp", "irpt", "trelease"]) + particles = data.groupby(["imdl", "iprp", "irpt", "trelease"]) + seqn_key = "sequencenumber" + data[seqn_key] = particles.ngroup() + + # convert to recarray + data = data.to_records(index=False) + + # define mp7 dtype + mp7_dtypes = np.dtype( + [ + ("particleid", np.int32), # same as sequencenumber + ("particlegroup", np.int32), + ( + seqn_key, + np.int32, + ), # mp7 sequencenumber (globally unique auto-generated ID) + ( + "particleidloc", + np.int32, + ), # mp7 particle ID (unique within a group, user-assigned or autogenerated) + ("time", np.float32), + ("x", np.float32), + ("y", np.float32), + ("z", np.float32), + ("k", np.int32), + ("node", np.int32), + ("xloc", np.float32), + ("yloc", np.float32), + ("zloc", np.float32), + ("stressperiod", np.int32), + ("timestep", np.int32), + ] + ) + + # build mp7 format dataframe + return pd.DataFrame( + np.core.records.fromarrays( + [ + data[seqn_key], + data["iprp"], + data[seqn_key], + data["irpt"], + data["t"], + data["x"], + data["y"], + data["z"], + data["ilay"], + data["icell"], + # todo local coords (xloc, yloc, zloc) + np.zeros(data.shape[0]), + np.zeros(data.shape[0]), + np.zeros(data.shape[0]), + data["kper"], + data["kstp"], + ], + dtype=mp7_dtypes, + ) + ) + + +def to_mp7_endpoints(data) -> pd.DataFrame: + """ + Convert MODFLOW 6 PRT pathline data to MODPATH 7 endpoint format. + + Parameters + ---------- + data : np.recarray or pd.DataFrame + MODFLOW 6 PRT pathline data + + Returns + ------- + pd.DataFrame + """ + + # convert to dataframe if needed + if isinstance(data, np.recarray): + data = pd.DataFrame(data) + + # assign a unique particle index column incrementing an integer + # for each unique combination of irpt, iprp, imdl, and trelease + data = data.sort_values(["imdl", "iprp", "irpt", "trelease"]) + particles = data.groupby(["imdl", "iprp", "irpt", "trelease"]) + seqn_key = "sequencenumber" + data[seqn_key] = particles.ngroup() + + # select startpoints and endpoints, sort by sequencenumber + startpts = ( + data.sort_values("t").groupby(seqn_key).head(1).sort_values(seqn_key) + ) + endpts = ( + data.sort_values("t").groupby(seqn_key).tail(1).sort_values(seqn_key) + ) + + # add new columns for: + # - initial coordinates + # - initial zone + # - initial node number + # - initial layer + starttimes = startpts["t"].to_numpy() + startxs = startpts["x"].to_numpy() + startys = startpts["y"].to_numpy() + startzs = startpts["z"].to_numpy() + startzones = startpts["izone"].to_numpy() + startnodes = startpts["icell"].to_numpy() + startlayers = startpts["ilay"].to_numpy() + conditions = [ + startpts[seqn_key].eq(row[seqn_key]) for i, row in startpts.iterrows() + ] + endpts["x0"] = np.select(conditions, startxs) + endpts["y0"] = np.select(conditions, startys) + endpts["z0"] = np.select(conditions, startzs) + endpts["zone0"] = np.select(conditions, startzones) + endpts["node0"] = np.select(conditions, startnodes) + endpts["k0"] = np.select(conditions, startlayers) + + # convert to recarray + endpts = endpts.to_records(index=False) + + # define mp7 dtype + mp7_dtype = np.dtype( + [ + ( + "particleid", + np.int32, + ), # mp7 sequencenumber (globally unique auto-generated ID) + ("particlegroup", np.int32), + ( + "particleidloc", + np.int32, + ), # mp7 particle ID (unique within a group, user-assigned or autogenerated) + ("status", np.int32), + ("time0", np.float32), + ("time", np.float32), + ("node0", np.int32), + ("k0", np.int32), + ("xloc0", np.float32), + ("yloc0", np.float32), + ("zloc0", np.float32), + ("x0", np.float32), + ("y0", np.float32), + ("z0", np.float32), + ("zone0", np.int32), + ("initialcellface", np.int32), + ("node", np.int32), + ("k", np.int32), + ("xloc", np.float32), + ("yloc", np.float32), + ("zloc", np.float32), + ("x", np.float32), + ("y", np.float32), + ("z", np.float32), + ("zone", np.int32), + ("cellface", np.int32), + ] + ) + + # build mp7 format dataframe + return pd.DataFrame( + np.core.records.fromarrays( + [ + endpts["sequencenumber"], + endpts["iprp"], + endpts["irpt"], + endpts["istatus"], + endpts["trelease"], + endpts["t"], + endpts["node0"], + endpts["k0"], + # todo initial local coords (xloc0, yloc0, zloc0) + np.zeros(endpts.shape[0]), + np.zeros(endpts.shape[0]), + np.zeros(endpts.shape[0]), + endpts["x0"], + endpts["y0"], + endpts["z0"], + endpts["zone0"], + np.zeros(endpts.shape[0]), # todo initial cell face? + endpts["icell"], + endpts["ilay"], + # todo local coords (xloc, yloc, zloc) + np.zeros(endpts.shape[0]), + np.zeros(endpts.shape[0]), + np.zeros(endpts.shape[0]), + endpts["x"], + endpts["y"], + endpts["z"], + endpts["izone"], + np.zeros(endpts.shape[0]), # todo cell face? + ], + dtype=mp7_dtype, + ) + ) diff --git a/flopy/utils/binaryfile.py b/flopy/utils/binaryfile.py index 9bd90b90d0..e78b290b45 100644 --- a/flopy/utils/binaryfile.py +++ b/flopy/utils/binaryfile.py @@ -11,7 +11,7 @@ import os import warnings from pathlib import Path -from typing import Optional, Union +from typing import List, Optional, Union import numpy as np @@ -1576,7 +1576,7 @@ def get_data( text=None, paknam=None, full3D=False, - ): + ) -> Union[List, np.ndarray]: """ Get data from the binary budget file.