From 264bb4fb1d2862f232325886035791b755f70bb5 Mon Sep 17 00:00:00 2001 From: w-bonelli Date: Thu, 23 Mar 2023 10:28:03 -0400 Subject: [PATCH] feat(PRT): add MF6 PRT support utilities * to_coords() method for particle data classes * to_mp7_[pathlines/endpoints] conversion fns * support PRT pathline data in plot routines * add PRT tutorial and example notebook --- .docs/Notebooks/mf6_prt_tutorial01.py | 304 ++++++++ .../modpath7_to_prt_migration_example.py | 420 ++++++++++ .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_model_dot_plot.py | 90 +++ autotest/test_modflow.py | 2 +- autotest/test_particledata.py | 221 +++++- 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 | 80 ++ flopy/modpath/mp7particledata.py | 318 +++++++- flopy/plot/crosssection.py | 133 +++- flopy/plot/map.py | 107 ++- flopy/plot/plotutil.py | 216 +++++ flopy/utils/binaryfile.py | 4 +- 23 files changed, 2807 insertions(+), 847 deletions(-) create mode 100644 .docs/Notebooks/mf6_prt_tutorial01.py create mode 100644 .docs/Notebooks/modpath7_to_prt_migration_example.py 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/.docs/Notebooks/mf6_prt_tutorial01.py b/.docs/Notebooks/mf6_prt_tutorial01.py new file mode 100644 index 0000000000..2b0232a004 --- /dev/null +++ b/.docs/Notebooks/mf6_prt_tutorial01.py @@ -0,0 +1,304 @@ +# --- +# jupyter: +# jupytext: +# notebook_metadata_filter: metadata +# text_representation: +# extension: .py +# format_name: light +# format_version: '1.5' +# jupytext_version: 1.15.0 +# kernelspec: +# display_name: Python 3 +# language: python +# name: python3 +# metadata: +# section: mf6 +# --- + +# # MODFLOW 6 PRT (particle-tracking) tutorial + +# This tutorial runs a GWF model then a PRT model +# in separate simulations via flow model interface. +# +# The grid is a 10x10 square with a single layer, +# the same flow system shown in the FloPy readme. +# +# Particles are released from the top left cell. +# +# ## Initial setup +# +# First, import dependencies, set up a temporary workspace, and define some filenames and model parameters. + + +import os +from pathlib import Path +from tempfile import TemporaryDirectory + +import flopy +import matplotlib.cm as cm +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +from flopy.utils import PathlineFile +from flopy.utils.binaryfile import HeadFile + +# workspace +temp_dir = TemporaryDirectory() +ws = Path(temp_dir.name) + +# model names +name = "prtfmi01" +gwfname = f"{name}_gwf" +prtname = f"{name}_prt" +mp7name = f"{name}_mp7" + +# output file names +gwf_budget_file = f"{gwfname}.bud" +gwf_head_file = f"{gwfname}.hds" +prt_track_file = f"{prtname}.trk" +prt_track_csv_file = f"{prtname}.trk.csv" +mp7_pathline_file = f"{mp7name}.mppth" + +# model info +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 + +# Define release points. We will release three particles from the top left cell of the model. + +# release points +releasepts = [ + # particle index, k, i, j, x, y, z + # (0-based indexing converted to 1-based for mf6 by flopy) + (0, 0, 0, 0, 0.25, 9.25, 0.5), + (1, 0, 0, 0, 0.5, 9.5, 0.5), + (2, 0, 0, 0, 0.75, 9.75, 0.5), +] + +# Particle-tracking models require a flow solution to solve for particle motion — PRT models may either run side-by-side with a GWF model with an exchange, or may consume the output of a previously run GWF model via flow model interface. +# +# In this tutorial we do the latter. First we define a GWF model and simulation, then a separate PRT model/simulation. +# +# Define a function to build a MODFLOW 6 GWF model. + +def build_gwf_sim(ws, mf6): + # create simulation + sim = flopy.mf6.MFSimulation( + sim_name=name, + exe_name=mf6, + version="mf6", + sim_ws=ws, + ) + + # 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 + 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 + + +# Define the PRT simulation. + +def build_prt_sim(ws, mf6): + # create simulation + sim = flopy.mf6.MFSimulation( + sim_name=name, + exe_name=mf6, + version="mf6", + sim_ws=ws, + ) + + # create tdis package + flopy.mf6.modflow.mftdis.ModflowTdis( + sim, + pname="tdis", + time_units="DAYS", + nper=nper, + perioddata=[(perlen, nstp, tsmult)], + ) + + # create prt model + prt = flopy.mf6.ModflowPrt(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 the flow model interface + flopy.mf6.ModflowPrtfmi( + prt, + packagedata=[ + ("GWFHEAD", gwf_head_file), + ("GWFBUDGET", gwf_budget_file), + ], + ) + + # add explicit model solution + ems = flopy.mf6.ModflowEms( + sim, + pname="ems", + filename=f"{prtname}.ems", + ) + sim.register_solution_package(ems, [prt.name]) + + return sim + + +# ## Running models +# +# Next, build and run the models. The flow model must run before the tracking model so the latter can consume the former's cell budget output file (configured above via flow model interface). + +# build mf6 models +gwfsim = build_gwf_sim(ws, "mf6") +prtsim = build_prt_sim(ws, "mf6") + +# run mf6 models +for sim in [gwfsim, prtsim]: + sim.write_simulation() + success, _ = sim.run_simulation() + assert success + +# Extract model and grid objects from the simulations for use with plots. + +# extract model objects +gwf = gwfsim.get_model(gwfname) +prt = prtsim.get_model(prtname) + +# extract model grid +mg = gwf.modelgrid + +# ## Inspecting results +# +# Finally we can load and plot results. First make sure the expected output files were created, then load pathlines, heads, budget, and specific discharge data. + +# check mf6 output files exist +assert (ws / gwf_budget_file).is_file() +assert (ws / gwf_head_file).is_file() +assert (ws / prt_track_file).is_file() +assert (ws / prt_track_csv_file).is_file() + +# load pathlines +pls = pd.read_csv(ws / prt_track_csv_file) + +# extract head, budget, and specific discharge results from GWF model +hds = HeadFile(ws / gwf_head_file).get_data() +bud = gwf.output.budget() +spdis = bud.get_data(text="DATA-SPDIS")[0] +qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(spdis, gwf) + +# Path points/lines can be plotted straightforwardly. + +# plot mf6 pathlines in map view +fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(13, 13)) +ax.set_aspect("equal") +pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax) +pmv.plot_grid() +pmv.plot_array(hds[0], alpha=0.1) +pmv.plot_vector(qx, qy, normalize=True, color="white") +pathlines = pls.groupby(["imdl", "iprp", "irpt", "trelease"]) +for ipl, ((imdl, iprp, irpt, trelease), pl) in enumerate(pathlines): + pl.plot( + title="Pathlines", + kind="line", + x="x", + y="y", + marker="o", + ax=ax, + legend=False, + color=cm.plasma(ipl / len(pathlines)), + ) + +# Alternatively, with FloPy's built-in pathline plotting utilities: + +fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(13, 13)) +ax.set_aspect("equal") +pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax) +pmv.plot_grid() +pmv.plot_array(hds[0], alpha=0.1) +pmv.plot_vector(qx, qy, normalize=True, color="white") +colors = cm.plasma(np.linspace(0, 1, pls.groupby(["imdl", "iprp", "irpt", "trelease"]).ngroups)) +pmv.plot_pathline(pls, layer="all", colors=colors, linewidth=2) diff --git a/.docs/Notebooks/modpath7_to_prt_migration_example.py b/.docs/Notebooks/modpath7_to_prt_migration_example.py new file mode 100644 index 0000000000..6ddd56d8c3 --- /dev/null +++ b/.docs/Notebooks/modpath7_to_prt_migration_example.py @@ -0,0 +1,420 @@ +# --- +# jupyter: +# jupytext: +# notebook_metadata_filter: metadata +# text_representation: +# extension: .py +# format_name: light +# format_version: '1.5' +# jupytext_version: 1.15.0 +# kernelspec: +# display_name: Python 3 +# language: python +# name: python3 +# metadata: +# section: mf6 +# --- + +# # Migrating from MODPATH 7 to MODFLOW 6 PRT + +# This example runs a GWF model, then a PRT model +# in separate simulations via flow model interface, +# then an identical MODPATH 7 model for comparison. +# +# The grid is a 10x10 square with a single layer, +# the same flow system shown in the FloPy readme. +# +# Particles are released from the top left cell. + +# First, import dependencies, set up a temporary workspace, and define some filenames and model parameters. + + +import os +from pathlib import Path +from tempfile import TemporaryDirectory + +import flopy +import matplotlib.cm as cm +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +from flopy.utils import PathlineFile +from flopy.utils.binaryfile import HeadFile + +# workspace +temp_dir = TemporaryDirectory() +ws = Path(temp_dir.name) + +# model names +name = "prtfmi01" +gwfname = f"{name}_gwf" +prtname = f"{name}_prt" +mp7name = f"{name}_mp7" + +# output file names +gwf_budget_file = f"{gwfname}.bud" +gwf_head_file = f"{gwfname}.hds" +prt_track_file = f"{prtname}.trk" +prt_track_csv_file = f"{prtname}.trk.csv" +mp7_pathline_file = f"{mp7name}.mppth" + +# model info +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 + + +# ## MODFLOW 6 GWF setup + +# Both the PRT and MP7 models depend on the same GWF model. + + +def build_gwf_sim(ws, mf6): + # create simulation + sim = flopy.mf6.MFSimulation( + sim_name=name, + exe_name=mf6, + version="mf6", + sim_ws=ws, + ) + + # 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 + 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 + + +# ## Release points + +# FloPy provides utilities to convert MODPATH 7 particle data +# to the appropriate format for MODFLOW 6 PRT PRP package input. + +# We specify release point locations in MODPATH 7 format, then +# convert from local to global coordinates for MODFLOW 6 PRT. + +# release points in mp7 format (using local coordinates) +releasepts_mp7 = [ + # node number, localx, localy, localz + # (0-based indexing converted to 1-based for mp7 by flopy) + (0, float(f"0.{i + 1}"), float(f"0.{i + 1}"), 0.5) + for i in range(9) +] + + +def get_partdata(grid): + return flopy.modpath.ParticleData( + partlocs=[grid.get_lrc(p[0])[0] for p in releasepts_mp7], + structured=True, + localx=[p[1] for p in releasepts_mp7], + localy=[p[2] for p in releasepts_mp7], + localz=[p[3] for p in releasepts_mp7], + timeoffset=0, + drape=0, + ) + + +# ## MODFLOW 6 PRT setup + +# We can now build the PRT simulation, using FloPy's `ParticleData.to_coords()` method to ease +# the conversion from MODPATH 7 to MODFLOW 6 PRT release points. + + +def build_prt_sim(ws, mf6): + # create simulation + sim = flopy.mf6.MFSimulation( + sim_name=name, + exe_name=mf6, + version="mf6", + sim_ws=ws, + ) + + # create tdis package + flopy.mf6.modflow.mftdis.ModflowTdis( + sim, + pname="tdis", + time_units="DAYS", + nper=nper, + perioddata=[(perlen, nstp, tsmult)], + ) + + # create prt model + prt = flopy.mf6.ModflowPrt(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) + + # convert mp7 particledata to prt release points + partdata = get_partdata(prt.modelgrid) + coords = partdata.to_coords(prt.modelgrid) + releasepts = [(i, 0, 0, 0, c[0], c[1], c[2]) for i, c in enumerate(coords)] + + # 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 the flow model interface + flopy.mf6.ModflowPrtfmi( + prt, + packagedata=[ + ("GWFHEAD", gwf_head_file), + ("GWFBUDGET", gwf_budget_file), + ], + ) + + # add explicit model solution + ems = flopy.mf6.ModflowEms( + sim, + pname="ems", + filename=f"{prtname}.ems", + ) + sim.register_solution_package(ems, [prt.name]) + + return sim + + +# ## MODPATH 7 setup + +# Finally we can build the MODPATH 7 model and simulation. + + +def build_mp7_sim(ws, mp7, gwf): + # convert mp7 particledata to prt release points + partdata = get_partdata(gwf.modelgrid) + + # create modpath 7 simulation + pg = flopy.modpath.ParticleGroup( + particlegroupname="G1", + particledata=partdata, + filename=f"{mp7name}.sloc", + ) + mp = flopy.modpath.Modpath7( + modelname=mp7name, + flowmodel=gwf, + exe_name=mp7, + model_ws=ws, + ) + mpbas = flopy.modpath.Modpath7Bas( + mp, + porosity=porosity, + ) + mpsim = flopy.modpath.Modpath7Sim( + mp, + simulationtype="pathline", + trackingdirection="forward", + budgetoutputoption="summary", + stoptimeoption="extend", + particlegroups=[pg], + ) + + return mp + +# ## Running the models + +# Construct and run the models. + +# build mf6 simulations +gwfsim = build_gwf_sim(ws, "mf6") +prtsim = build_prt_sim(ws, "mf6") + +# extract models and grid (useful for plotting etc) +gwf = gwfsim.get_model(gwfname) +prt = prtsim.get_model(prtname) +mg = gwf.modelgrid + +# build mp7 model +mp7sim = build_mp7_sim(ws, "mp7", gwf) + +# run mf6 models +for sim in [gwfsim, prtsim]: + sim.write_simulation() + success, _ = sim.run_simulation() + assert success + +# run mp7 model +mp7sim.write_input() +success, _ = mp7sim.run_model() +assert success + +# ## Inspecting the results + +# Make sure the expected output files exist. + +# check mf6 output files exist +assert (ws / gwf_budget_file).is_file() +assert (ws / gwf_head_file).is_file() +assert (ws / prt_track_file).is_file() +assert (ws / prt_track_csv_file).is_file() + +# check mp7 output files exist +assert (ws / mp7_pathline_file).is_file() + +# Load results from MODPATH 7 and MODFLOW 6 output files. + +# load mf6 pathlines +pls = pd.read_csv(ws / prt_track_csv_file) + +# load mp7 pathlines +plf = PathlineFile(ws / mp7_pathline_file) +mp7_pldata = pd.DataFrame( + plf.get_destination_pathline_data(range(mg.nnodes), to_recarray=True) +) + +# convert mp7 pathline fields from 0- to 1-based indexing +mp7_pldata["particleid"] = mp7_pldata["particleid"] + 1 +mp7_pldata["particlegroup"] = mp7_pldata["particlegroup"] + 1 +mp7_pldata["node"] = mp7_pldata["node"] + 1 +mp7_pldata["k"] = mp7_pldata["k"] + 1 + +# load head, budget, and specific discharge from gwf model +hds = HeadFile(ws / gwf_head_file).get_data() +bud = gwf.output.budget() +spdis = bud.get_data(text="DATA-SPDIS")[0] +qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(spdis, gwf) + +# Finally, plot the MODPATH 7 and MODFLOW 6 pathlines side-by-side in map view. + +# + +# plot mf6 and mp7 pathlines in map view +fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(13, 13)) +for a in ax: + a.set_aspect("equal") + +pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax[0]) +pmv.plot_grid() +pmv.plot_array(hds[0], alpha=0.1) +pmv.plot_vector(qx, qy, normalize=True, color="white") +pathlines = pls.groupby(["imdl", "iprp", "irpt", "trelease"]) +for ipl, ((imdl, iprp, irpt, trelease), pl) in enumerate(pathlines): + pl.plot( + title="MF6 pathlines", + kind="line", + x="x", + y="y", + marker="o", + ax=ax[0], + legend=False, + color=cm.plasma(ipl / len(pathlines)), + ) + +pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax[1]) +pmv.plot_grid() +pmv.plot_array(hds[0], alpha=0.1) +pmv.plot_vector(qx, qy, normalize=True, color="white") +mp7_plines = mp7_pldata.groupby(["particleid"]) +for ipl, (pid, pl) in enumerate(mp7_plines): + pl.plot( + title="MP7 pathlines", + kind="line", + x="x", + y="y", + marker="o", + ax=ax[1], + legend=False, + color=cm.plasma(ipl / len(mp7_plines)), + ) +# - + +# Alternatively, FloPy built-in plotting functions can be used to plot pathlines. + +# + +fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(13, 13)) +for a in ax: + a.set_aspect("equal") + +colors = cm.plasma(np.linspace(0, 1, len(mp7_pldata.particleid.unique()))) + +pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax[0]) +pmv.plot_grid() +pmv.plot_array(hds[0], alpha=0.1) +pmv.plot_vector(qx, qy, normalize=True, color="white") +pmv.plot_pathline(pls, layer="all", colors=colors) +ax[0].set_title("MP7 pathlines") + +pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax[1]) +pmv.plot_grid() +pmv.plot_array(hds[0], alpha=0.1) +pmv.plot_vector(qx, qy, normalize=True, color="white") +pmv.plot_pathline(mp7_pldata, layer="all", colors=colors) +ax[1].set_title("MP7 pathlines") diff --git a/.github/workflows/commit.yml b/.github/workflows/commit.yml index 770eb2460b..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 -c 'import flopy; flopy.mf6.utils.generate_classes(ref="develop", backup=False)' + # 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 -c 'import flopy; flopy.mf6.utils.generate_classes(ref="develop", backup=False)' + # 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 739c408e37..593cbd0d2c 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_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_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..e1a640d647 100644 --- a/autotest/test_particledata.py +++ b/autotest/test_particledata.py @@ -1,8 +1,17 @@ import numpy as np +import pytest +from autotest.test_grid_cases import GridCases -from flopy.modpath import ParticleData +from flopy.discretization import StructuredGrid +from flopy.modpath import ( + CellDataType, + FaceDataType, + LRCParticleData, + NodeParticleData, + ParticleData, +) -structured_plocs = [(1, 1, 1), (1, 1, 2)] +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..03cc1f4294 --- /dev/null +++ b/autotest/test_plotutil.py @@ -0,0 +1,80 @@ +from modflow_devtools.markers import requires_exe + +# test PRT-MP7 pathline conversion functions + + +# @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_to_mp7_pathlines(mf6_prt_sim): +# pass +# +# +# @requires_exe("mf6") +# def test_to_mp7_endpoints(mf6_pr_sim): +# pass diff --git a/flopy/modpath/mp7particledata.py b/flopy/modpath/mp7particledata.py index 9da7e57490..a780a46383 100644 --- a/flopy/modpath/mp7particledata.py +++ b/flopy/modpath/mp7particledata.py @@ -4,6 +4,7 @@ """ +from itertools import product import numpy as np from numpy.lib.recfunctions import unstructured_to_structured @@ -316,19 +317,14 @@ def __init__( self.locationstyle = locationstyle self.particledata = particledata - return - 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"): @@ -358,7 +354,69 @@ def write(self, f=None): for v in d: f.write(fmt.format(*v)) - return + def to_coords(self, grid): + """ + Convert the particle representation to a list of global coordinates. + + Parameters + ---------- + grid : The grid to use for locating particle coordinates + + Returns + ------- + A list of particle tuples ([optional particle ID], x coord, y coord, z coord) + where particle ID is only included if provided to the particle data instance. + """ + + if grid.grid_type == "structured": + if not hasattr(self.particledata, "k"): + raise ValueError( + f"Particle representation is not structured but grid is" + ) + + def to_x(p): + i, j = p[1], p[2] + local = p[3] + verts = grid.get_cell_vertices(i, j) + xs, _ = list(zip(*verts)) + minx, maxx = min(xs), max(xs) + span = maxx - minx + return minx + span * local + + def to_y(p): + i, j = p[1], p[2] + local = p[4] + verts = grid.get_cell_vertices(i, j) + _, ys = list(zip(*verts)) + miny, maxy = min(ys), max(ys) + span = maxy - miny + return miny + span * local + + def to_z(p): + k, i, j = p[0], p[1], p[2] + local = p[5] + minz, maxz = grid.botm[k, i, j], grid.top[i, j] + span = maxz - minz + return minz + span * local + + pdl = self.particledata.tolist() + if hasattr(self.particledata, "id"): + ids = [p[0] for p in pdl] + particles = [r[1:7] for r in pdl] + return [ + (id, to_x(p), to_y(p), to_z(p)) + for id, p in zip(ids, particles) + ] + else: + particles = [r[0:6] for r in pdl] + return [(to_x(p), to_y(p), to_z(p)) for p in particles] + else: + if hasattr(self.particledata, "k"): + raise ValueError( + f"Particle representation is structured but grid is not" + ) + + raise NotImplementedError() def _get_dtype(self, structured, particleid): """ @@ -542,7 +600,6 @@ def __init__( self.columndivisions5 = columndivisions5 self.rowdivisions6 = rowdivisions6 self.columndivisions6 = columndivisions6 - return def write(self, f=None): """ @@ -637,7 +694,6 @@ def __init__( self.columncelldivisions = columncelldivisions self.rowcelldivisions = rowcelldivisions self.layercelldivisions = layercelldivisions - return def write(self, f=None): """ @@ -668,7 +724,176 @@ def write(self, f=None): ) f.write(line) - return + +def get_cell_particles(subdivisiondata, grid, k=None, i=None, j=None, nn=None): + cell_particles = [] + + # 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] + else: + minz, maxz = grid.botm[nn], grid.top[nn] + 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])) + cell_particles = cell_particles + [ + (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])) + cell_particles = cell_particles + [ + (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])) + cell_particles = cell_particles + [ + (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])) + cell_particles = cell_particles + [ + (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])) + cell_particles = cell_particles + [ + (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])) + cell_particles = cell_particles + [ + (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])) + cell_particles = cell_particles + prod + else: + raise ValueError( + f"Unsupported subdivision data type: {type(subdivisiondata)}" + ) + + return cell_particles class LRCParticleData: @@ -686,7 +911,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 +924,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 +965,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 +983,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 +995,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 +1047,37 @@ def write(self, f=None): return + def to_coords(self, grid): + """ + Convert the particle representation to a list of global coordinates. + + Parameters + ---------- + grid : The grid to use for locating particle coordinates + + Returns + ------- + A list of particle tuples (particle ID, x coord, y coord, z coord) + """ + + if grid.grid_type != "structured": + raise ValueError( + f"Particle representation is structured but grid is not" + ) + + particles = [] + 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: + particles = particles + get_cell_particles( + sd, grid, k, i, j + ) + + return particles + class NodeParticleData: """ @@ -986,3 +1240,23 @@ def write(self, f=None): f.write(line) return + + def to_coords(self, grid): + """ + Convert the particle representation to a list of global coordinates. + + Parameters + ---------- + grid : The grid to use for locating particle coordinates + + Returns + ------- + A list of particle tuples (particle ID, x coord, y coord, z coord) + """ + + particles = [] + for nn in self.nodedata: + for sd in self.subdivisiondata: + particles = particles + get_cell_particles(sd, grid, nn=nn) + + return particles 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.