diff --git a/autotest/test.py b/autotest/test.py new file mode 100644 index 000000000..f20d464e6 --- /dev/null +++ b/autotest/test.py @@ -0,0 +1,87 @@ +import os +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np +import pytest +from modflow_devtools.markers import requires_exe, requires_pkg +from test_mp7_cases import Mp7Cases + +from flopy.mf6 import ( + MFSimulation, + ModflowGwf, + ModflowGwfdis, + ModflowGwfic, + ModflowGwfnpf, + ModflowGwfoc, + ModflowGwfrcha, + ModflowGwfriv, + ModflowGwfwel, + ModflowIms, + ModflowTdis, +) +from flopy.modpath import ( + CellDataType, + FaceDataType, + LRCParticleData, + Modpath7, + Modpath7Bas, + Modpath7Sim, + NodeParticleData, + ParticleData, + ParticleGroup, + ParticleGroupLRCTemplate, + ParticleGroupNodeTemplate, +) +from flopy.plot import PlotMapView +from flopy.utils import EndpointFile, PathlineFile + +function_tmpdir = Path("temp/test") + +mf6sim = Mp7Cases.mf6(function_tmpdir) +mf6sim.write_simulation() +mf6sim.run_simulation() + +# create mp7 model +mp = Modpath7( + modelname=f"{mf6sim.name}_mp", + flowmodel=mf6sim.get_model(mf6sim.name), + exe_name="mp7", + model_ws=mf6sim.sim_path, + verbose=True, +) +defaultiface6 = {"RCH": 6, "EVT": 6} +mpbas = Modpath7Bas(mp, porosity=0.1, defaultiface=defaultiface6) +mpsim = Modpath7Sim( + mp, + simulationtype="combined", + trackingdirection="forward", + weaksinkoption="pass_through", + weaksourceoption="pass_through", + budgetoutputoption="summary", + budgetcellnumbers=[1049, 1259], + traceparticledata=[1, 1000], + referencetime=[0, 0, 0.0], + stoptimeoption="extend", + timepointdata=[500, 1000.0], + zonedataoption="on", + zones=Mp7Cases.zones, + particlegroups=Mp7Cases.particlegroups, +) +# add a duplicate mp7sim package +mpsim = Modpath7Sim( + mp, + simulationtype="combined", + trackingdirection="forward", + weaksinkoption="pass_through", + weaksourceoption="pass_through", + budgetoutputoption="summary", + budgetcellnumbers=[1049, 1259], + traceparticledata=[1, 1000], + referencetime=[0, 0, 0.0], + stoptimeoption="extend", + timepointdata=[500, 1000.0], + zonedataoption="on", + zones=Mp7Cases.zones, + particlegroups=Mp7Cases.particlegroups, +) diff --git a/autotest/test_binaryfile.py b/autotest/test_binaryfile.py index bb6795b70..d0d592520 100644 --- a/autotest/test_binaryfile.py +++ b/autotest/test_binaryfile.py @@ -6,12 +6,14 @@ from matplotlib import pyplot as plt from matplotlib.axes import Axes +from flopy.mf6.modflow import MFSimulation from flopy.modflow import Modflow from flopy.utils import ( BinaryHeader, CellBudgetFile, HeadFile, HeadUFile, + UcnFile, Util2d, ) from flopy.utils.binaryfile import ( @@ -197,6 +199,9 @@ def test_headufile_get_ts(example_data_path): heads.get_ts([i, i + 1, i + 2]) +# precision + + def test_get_headfile_precision(example_data_path): precision = get_headfile_precision( example_data_path / "freyberg" / "freyberg.githds" @@ -249,6 +254,9 @@ def test_budgetfile_detect_precision_double(path): assert file.realtype == np.float64 +# write + + def test_write_head(function_tmpdir): file_path = function_tmpdir / "headfile" head_data = np.random.random((10, 10)) @@ -286,6 +294,9 @@ def test_write_budget(function_tmpdir): assert np.array_equal(content1, content2) +# read + + def test_binaryfile_read(function_tmpdir, freyberg_model_path): h = HeadFile(freyberg_model_path / "freyberg.githds") assert isinstance(h, HeadFile) @@ -336,212 +347,59 @@ def test_binaryfile_read_context(freyberg_model_path): assert str(e.value) == "seek of closed file", str(e.value) -def test_cellbudgetfile_read_context(example_data_path): - mf2005_model_path = example_data_path / "mf2005_test" - cbc_path = mf2005_model_path / "mnw1.gitcbc" - with CellBudgetFile(cbc_path) as v: - data = v.get_data(text="DRAINS")[0] - assert data.min() < 0, data.min() - assert not v.file.closed - assert v.file.closed - - with pytest.raises(ValueError) as e: - v.get_data(text="DRAINS") - assert str(e.value) == "seek of closed file", str(e.value) - - -def test_cellbudgetfile_read(example_data_path): - mf2005_model_path = example_data_path / "mf2005_test" - v = CellBudgetFile(mf2005_model_path / "mnw1.gitcbc") - assert isinstance(v, CellBudgetFile) - - kstpkper = v.get_kstpkper() - assert len(kstpkper) == 5, "length of kstpkper != 5" - - records = v.get_unique_record_names() - idx = 0 - for t in kstpkper: - for record in records: - t0 = v.get_data(kstpkper=t, text=record, full3D=True)[0] - t1 = v.get_data(idx=idx, text=record, full3D=True)[0] - assert np.array_equal(t0, t1), ( - f"binary budget item {record} read using kstpkper != binary " - f"budget item {record} read using idx" - ) - idx += 1 - v.close() - - -def test_cellbudgetfile_position(function_tmpdir, zonbud_model_path): - fpth = zonbud_model_path / "freyberg.gitcbc" - v = CellBudgetFile(fpth) - assert isinstance(v, CellBudgetFile) - - # starting position of data - idx = 8767 - ipos = v.get_position(idx) - ival = 50235424 - assert ipos == ival, f"position of index 8767 != {ival}" - - ipos = v.get_position(idx, header=True) - ival = 50235372 - assert ipos == ival, f"position of index 8767 header != {ival}" - - cbcd = [] - for i in range(idx, v.get_nrecords()): - cbcd.append(v.get_data(i)[0]) - v.close() - - # write the last entry as a new binary file - fin = open(fpth, "rb") - fin.seek(ipos) - length = os.path.getsize(fpth) - ipos - - buffsize = 32 - opth = str(function_tmpdir / "end.cbc") - with open(opth, "wb") as fout: - while length: - chunk = min(buffsize, length) - data = fin.read(chunk) - fout.write(data) - length -= chunk - fin.close() - - v2 = CellBudgetFile(opth, verbose=True) - - try: - v2.list_records() - except: - assert False, f"could not list records on {opth}" - - names = v2.get_unique_record_names(decode=True) - - cbcd2 = [] - for i in range(0, v2.get_nrecords()): - cbcd2.append(v2.get_data(i)[0]) - v2.close() - - for i, (d1, d2) in enumerate(zip(cbcd, cbcd2)): - msg = f"{names[i].rstrip()} data from slice is not identical" - assert np.array_equal(d1, d2), msg - - # Check error when reading empty file - fname = function_tmpdir / "empty.gitcbc" - with open(fname, "w"): - pass - with pytest.raises(ValueError): - CellBudgetFile(fname) - - -def test_cellbudgetfile_readrecord(example_data_path): - mf2005_model_path = example_data_path / "mf2005_test" - cbc_fname = mf2005_model_path / "test1tr.gitcbc" - v = CellBudgetFile(cbc_fname) - assert isinstance(v, CellBudgetFile) - - kstpkper = v.get_kstpkper() - assert len(kstpkper) == 30, "length of kstpkper != 30" - - with pytest.raises(TypeError) as e: - v.get_data() - assert str(e.value).startswith( - "get_data() missing 1 required argument" - ), str(e.exception) - - t = v.get_data(text="STREAM LEAKAGE") - assert len(t) == 30, "length of stream leakage data != 30" - assert ( - t[0].shape[0] == 36 - ), "sfr budget data does not have 36 reach entries" - - t = v.get_data(text="STREAM LEAKAGE", full3D=True) - assert t[0].shape == (1, 15, 10), ( - "3D sfr budget data does not have correct shape (1, 15,10) - " - "returned shape {}".format(t[0].shape) - ) +# reverse - for kk in kstpkper: - t = v.get_data(kstpkper=kk, text="STREAM LEAKAGE", full3D=True)[0] - assert t.shape == (1, 15, 10), ( - "3D sfr budget data for kstpkper {} " - "does not have correct shape (1, 15,10) - " - "returned shape {}".format(kk, t[0].shape) - ) - - idx = v.get_indices() - assert idx is None, "get_indices() without record did not return None" - - records = v.get_unique_record_names() - for record in records: - indices = v.get_indices(text=record.strip()) - for idx, kk in enumerate(kstpkper): - t0 = v.get_data(kstpkper=kk, text=record.strip())[0] - t1 = v.get_data(idx=indices[idx], text=record)[0] - assert np.array_equal( - t0, t1 - ), "binary budget item {0} read using kstpkper != binary budget item {0} read using idx".format( - record - ) - - # idx can be either an int or a list of ints - s9 = v.get_data(idx=9) - assert len(s9) == 1 - s09 = v.get_data(idx=[0, 9]) - assert len(s09) == 2 - assert (s09[1] == s9).all() - - v.close() - - -def test_cellbudgetfile_readrecord_waux(example_data_path): - mf2005_model_path = example_data_path / "mf2005_test" - cbc_fname = mf2005_model_path / "test1tr.gitcbc" - v = CellBudgetFile(cbc_fname) - assert isinstance(v, CellBudgetFile) - - kstpkper = v.get_kstpkper() - assert len(kstpkper) == 30, "length of kstpkper != 30" - - t = v.get_data(text="WELLS") - assert len(t) == 30, "length of well data != 30" - assert t[0].shape[0] == 10, "wel budget data does not have 10 well entries" - assert t[0].dtype.names == ("node", "q", "IFACE") - np.testing.assert_array_equal( - t[0]["node"], - [54, 55, 64, 65, 74, 75, 84, 85, 94, 95], - ) - np.testing.assert_array_equal(t[0]["q"], np.repeat(np.float32(-10.0), 10)) - np.testing.assert_array_equal( - t[0]["IFACE"], - np.array([1, 2, 3, 4, 5, 6, 0, 0, 0, 0], np.float32), - ) - t = v.get_data(text="WELLS", full3D=True) - assert t[0].shape == (1, 15, 10), ( - "3D wel budget data does not have correct shape (1, 15,10) - " - "returned shape {}".format(t[0].shape) +def test_headfile_reverse_mf6(example_data_path, function_tmpdir): + # load simulation and extract tdis + sim_name = "test006_gwf3" + sim = MFSimulation.load( + sim_name=sim_name, sim_ws=example_data_path / "mf6" / sim_name ) - - for kk in kstpkper: - t = v.get_data(kstpkper=kk, text="wells", full3D=True)[0] - assert t.shape == (1, 15, 10), ( - "3D wel budget data for kstpkper {} " - "does not have correct shape (1, 15,10) - " - "returned shape {}".format(kk, t[0].shape) - ) - - idx = v.get_indices() - assert idx is None, "get_indices() without record did not return None" - - records = v.get_unique_record_names() - for record in records: - indices = v.get_indices(text=record.strip()) - for idx, kk in enumerate(kstpkper): - t0 = v.get_data(kstpkper=kk, text=record.strip())[0] - t1 = v.get_data(idx=indices[idx], text=record)[0] - assert np.array_equal( - t0, t1 - ), "binary budget item {0} read using kstpkper != binary budget item {0} read using idx".format( - record - ) - v.close() + tdis = sim.get_package("tdis") + + # load cell budget file, providing tdis as kwarg + model_path = example_data_path / "mf6" / sim_name + file_stem = "flow_adj" + file_path = model_path / "expected_output" / f"{file_stem}.hds" + f = HeadFile(file_path, tdis=tdis) + assert isinstance(f, HeadFile) + + # reverse the file + rf_name = f"{file_stem}_rev.hds" + f.reverse(filename=function_tmpdir / rf_name) + rf = HeadFile(function_tmpdir / rf_name) + assert isinstance(rf, HeadFile) + + # check that data from both files have the same shape + f_data = f.get_alldata() + f_shape = f_data.shape + rf_data = rf.get_alldata() + rf_shape = rf_data.shape + assert f_shape == rf_shape + + # check number of records + nrecords = f.get_nrecords() + assert nrecords == rf.get_nrecords() + + # check that the data are reversed + for idx in range(nrecords - 1, -1, -1): + # check headers + f_header = list(f.recordarray[nrecords - idx - 1]) + rf_header = list(rf.recordarray[idx]) + f_totim = f_header.pop(9) # todo check totim + rf_totim = rf_header.pop(9) + assert f_header == rf_header + assert f_header == rf_header + + # check data + f_data = f.get_data(idx=idx)[0] + rf_data = rf.get_data(idx=nrecords - idx - 1)[0] + assert f_data.shape == rf_data.shape + if f_data.ndim == 1: + for row in range(len(f_data)): + f_datum = f_data[row] + rf_datum = rf_data[row] + assert f_datum == rf_datum + else: + assert np.array_equal(f_data[0][0], rf_data[0][0]) diff --git a/autotest/test_cellbudgetfile.py b/autotest/test_cellbudgetfile.py new file mode 100644 index 000000000..175df6037 --- /dev/null +++ b/autotest/test_cellbudgetfile.py @@ -0,0 +1,308 @@ +import os + +import numpy as np +import pytest + +from flopy.mf6.modflow.mfsimulation import MFSimulation +from flopy.utils.binaryfile import CellBudgetFile + + +@pytest.fixture +def zonbud_model_path(example_data_path): + return example_data_path / "zonbud_examples" + + +def test_cellbudgetfile_position(function_tmpdir, zonbud_model_path): + fpth = zonbud_model_path / "freyberg.gitcbc" + v = CellBudgetFile(fpth) + assert isinstance(v, CellBudgetFile) + + # starting position of data + idx = 8767 + ipos = v.get_position(idx) + ival = 50235424 + assert ipos == ival, f"position of index 8767 != {ival}" + + ipos = v.get_position(idx, header=True) + ival = 50235372 + assert ipos == ival, f"position of index 8767 header != {ival}" + + cbcd = [] + for i in range(idx, v.get_nrecords()): + cbcd.append(v.get_data(i)[0]) + v.close() + + # write the last entry as a new binary file + fin = open(fpth, "rb") + fin.seek(ipos) + length = os.path.getsize(fpth) - ipos + + buffsize = 32 + opth = str(function_tmpdir / "end.cbc") + with open(opth, "wb") as fout: + while length: + chunk = min(buffsize, length) + data = fin.read(chunk) + fout.write(data) + length -= chunk + fin.close() + + v2 = CellBudgetFile(opth, verbose=True) + + try: + v2.list_records() + except: + assert False, f"could not list records on {opth}" + + names = v2.get_unique_record_names(decode=True) + + cbcd2 = [] + for i in range(0, v2.get_nrecords()): + cbcd2.append(v2.get_data(i)[0]) + v2.close() + + for i, (d1, d2) in enumerate(zip(cbcd, cbcd2)): + msg = f"{names[i].rstrip()} data from slice is not identical" + assert np.array_equal(d1, d2), msg + + # Check error when reading empty file + fname = function_tmpdir / "empty.gitcbc" + with open(fname, "w"): + pass + with pytest.raises(ValueError): + CellBudgetFile(fname) + + +# read context + + +def test_cellbudgetfile_read_context(example_data_path): + mf2005_model_path = example_data_path / "mf2005_test" + cbc_path = mf2005_model_path / "mnw1.gitcbc" + with CellBudgetFile(cbc_path) as v: + data = v.get_data(text="DRAINS")[0] + assert data.min() < 0, data.min() + assert not v.file.closed + assert v.file.closed + + with pytest.raises(ValueError) as e: + v.get_data(text="DRAINS") + assert str(e.value) == "seek of closed file", str(e.value) + + +# read + + +def test_cellbudgetfile_read(example_data_path): + mf2005_model_path = example_data_path / "mf2005_test" + v = CellBudgetFile(mf2005_model_path / "mnw1.gitcbc") + assert isinstance(v, CellBudgetFile) + + kstpkper = v.get_kstpkper() + assert len(kstpkper) == 5, "length of kstpkper != 5" + + records = v.get_unique_record_names() + idx = 0 + for t in kstpkper: + for record in records: + t0 = v.get_data(kstpkper=t, text=record, full3D=True)[0] + t1 = v.get_data(idx=idx, text=record, full3D=True)[0] + assert np.array_equal(t0, t1), ( + f"binary budget item {record} read using kstpkper != binary " + f"budget item {record} read using idx" + ) + idx += 1 + v.close() + + +# readrecord + + +def test_cellbudgetfile_readrecord(example_data_path): + mf2005_model_path = example_data_path / "mf2005_test" + cbc_fname = mf2005_model_path / "test1tr.gitcbc" + v = CellBudgetFile(cbc_fname) + assert isinstance(v, CellBudgetFile) + + kstpkper = v.get_kstpkper() + assert len(kstpkper) == 30, "length of kstpkper != 30" + + with pytest.raises(TypeError) as e: + v.get_data() + assert str(e.value).startswith( + "get_data() missing 1 required argument" + ), str(e.exception) + + t = v.get_data(text="STREAM LEAKAGE") + assert len(t) == 30, "length of stream leakage data != 30" + assert ( + t[0].shape[0] == 36 + ), "sfr budget data does not have 36 reach entries" + + t = v.get_data(text="STREAM LEAKAGE", full3D=True) + assert t[0].shape == (1, 15, 10), ( + "3D sfr budget data does not have correct shape (1, 15,10) - " + "returned shape {}".format(t[0].shape) + ) + + for kk in kstpkper: + t = v.get_data(kstpkper=kk, text="STREAM LEAKAGE", full3D=True)[0] + assert t.shape == (1, 15, 10), ( + "3D sfr budget data for kstpkper {} " + "does not have correct shape (1, 15,10) - " + "returned shape {}".format(kk, t[0].shape) + ) + + idx = v.get_indices() + assert idx is None, "get_indices() without record did not return None" + + records = v.get_unique_record_names() + for record in records: + indices = v.get_indices(text=record.strip()) + for idx, kk in enumerate(kstpkper): + t0 = v.get_data(kstpkper=kk, text=record.strip())[0] + t1 = v.get_data(idx=indices[idx], text=record)[0] + assert np.array_equal( + t0, t1 + ), "binary budget item {0} read using kstpkper != binary budget item {0} read using idx".format( + record + ) + + # idx can be either an int or a list of ints + s9 = v.get_data(idx=9) + assert len(s9) == 1 + s09 = v.get_data(idx=[0, 9]) + assert len(s09) == 2 + assert (s09[1] == s9).all() + + v.close() + + +def test_cellbudgetfile_readrecord_waux(example_data_path): + mf2005_model_path = example_data_path / "mf2005_test" + cbc_fname = mf2005_model_path / "test1tr.gitcbc" + v = CellBudgetFile(cbc_fname) + assert isinstance(v, CellBudgetFile) + + kstpkper = v.get_kstpkper() + assert len(kstpkper) == 30, "length of kstpkper != 30" + + t = v.get_data(text="WELLS") + assert len(t) == 30, "length of well data != 30" + assert t[0].shape[0] == 10, "wel budget data does not have 10 well entries" + assert t[0].dtype.names == ("node", "q", "IFACE") + np.testing.assert_array_equal( + t[0]["node"], + [54, 55, 64, 65, 74, 75, 84, 85, 94, 95], + ) + np.testing.assert_array_equal(t[0]["q"], np.repeat(np.float32(-10.0), 10)) + np.testing.assert_array_equal( + t[0]["IFACE"], + np.array([1, 2, 3, 4, 5, 6, 0, 0, 0, 0], np.float32), + ) + + t = v.get_data(text="WELLS", full3D=True) + assert t[0].shape == (1, 15, 10), ( + "3D wel budget data does not have correct shape (1, 15,10) - " + "returned shape {}".format(t[0].shape) + ) + + for kk in kstpkper: + t = v.get_data(kstpkper=kk, text="wells", full3D=True)[0] + assert t.shape == (1, 15, 10), ( + "3D wel budget data for kstpkper {} " + "does not have correct shape (1, 15,10) - " + "returned shape {}".format(kk, t[0].shape) + ) + + idx = v.get_indices() + assert idx is None, "get_indices() without record did not return None" + + records = v.get_unique_record_names() + for record in records: + indices = v.get_indices(text=record.strip()) + for idx, kk in enumerate(kstpkper): + t0 = v.get_data(kstpkper=kk, text=record.strip())[0] + t1 = v.get_data(idx=indices[idx], text=record)[0] + assert np.array_equal( + t0, t1 + ), "binary budget item {0} read using kstpkper != binary budget item {0} read using idx".format( + record + ) + v.close() + + +# reverse + + +@pytest.mark.skip( + reason="failing, need to modify CellBudgetFile.reverse to support mf2005?" +) +def test_cellbudgetfile_reverse_mf2005(example_data_path, function_tmpdir): + sim_name = "test1tr" + + # load simulation and extract tdis + sim = MFSimulation.load( + sim_name=sim_name, sim_ws=example_data_path / "mf2005_test" + ) + tdis = sim.get_package("tdis") + + mf2005_model_path = example_data_path / sim_name + cbc_fname = mf2005_model_path / f"{sim_name}.gitcbc" + f = CellBudgetFile(cbc_fname, tdis=tdis) + assert isinstance(f, CellBudgetFile) + + rf_name = "test1tr_rev.gitcbc" + f.reverse(function_tmpdir / rf_name) + rf = CellBudgetFile(function_tmpdir / rf_name) + assert isinstance(rf, CellBudgetFile) + + +def test_cellbudgetfile_reverse_mf6(example_data_path, function_tmpdir): + # load simulation and extract tdis + sim_name = "test006_gwf3" + sim = MFSimulation.load( + sim_name=sim_name, sim_ws=example_data_path / "mf6" / sim_name + ) + tdis = sim.get_package("tdis") + + # load cell budget file, providing tdis as kwarg + model_path = example_data_path / "mf6" / sim_name + file_stem = "flow_adj" + file_path = model_path / "expected_output" / f"{file_stem}.cbc" + f = CellBudgetFile(file_path, tdis=tdis) + assert isinstance(f, CellBudgetFile) + + # reverse the file + rf_name = f"{file_stem}_rev.cbc" + f.reverse(filename=function_tmpdir / rf_name) + rf = CellBudgetFile(function_tmpdir / rf_name) + assert isinstance(rf, CellBudgetFile) + + # check that both files have the same number of records + nrecords = f.get_nrecords() + assert nrecords == rf.get_nrecords() + + # check data were reversed + for idx in range(nrecords - 1, -1, -1): + # check headers + f_header = list(f.recordarray[nrecords - idx - 1]) + rf_header = list(rf.recordarray[idx]) + f_totim = f_header.pop(9) # todo check totim + rf_totim = rf_header.pop(9) + assert f_header == rf_header + + # check data + f_data = f.get_data(idx=idx)[0] + rf_data = rf.get_data(idx=nrecords - idx - 1)[0] + assert f_data.shape == rf_data.shape + if f_data.ndim == 1: + for row in range(len(f_data)): + f_datum = f_data[row] + rf_datum = rf_data[row] + # flows should be negated + rf_datum[2] = -rf_datum[2] + assert f_datum == rf_datum + else: + # flows should be negated + assert np.array_equal(f_data[0][0], -rf_data[0][0]) diff --git a/flopy/utils/binaryfile.py b/flopy/utils/binaryfile.py index db7ffaf9f..07445cbcb 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 Union +from typing import Optional, Union import numpy as np @@ -165,16 +165,14 @@ def write_budget( class BinaryHeader(Header): """ - The binary_header class is a class to create headers for MODFLOW - binary files. + Represents data headers for binary output files. Parameters ---------- bintype : str - is the type of file being opened (head and ucn file currently - supported) + Type of file being opened. Accepted values are 'head' and 'ucn'. precision : str - is the precision of the floating point data in the file + Precision of floating point data in the file. """ @@ -423,9 +421,16 @@ def get_headfile_precision(filename: Union[str, os.PathLike]): class BinaryLayerFile(LayerFile): """ - The BinaryLayerFile class is the super class from which specific derived - classes are formed. This class should not be instantiated directly + The BinaryLayerFile class is a parent class from which concrete + classes inherit. This class should not be instantiated directly. + Notes + ----- + + The BinaryLayerFile class is built on a record array consisting of + headers, which are record arrays of the modflow header information + (kstp, kper, pertim, totim, text, nrow, ncol, ilay), and long ints + pointing to the 1st byte of data for the corresponding data arrays. """ def __init__( @@ -578,39 +583,25 @@ def get_ts(self, idx): class HeadFile(BinaryLayerFile): """ - HeadFile Class. + The HeadFile class provides simple ways to retrieve and manipulate + 2D or 3D head arrays, or time series arrays for one or more cells, + from a binary head output file. A utility method is also provided + to reverse the order of head data, for use with particle tracking + simulations in which particles are tracked backwards in time from + terminating to release locations (e.g., to compute capture zones). Parameters ---------- filename : str or PathLike - Path of the concentration file + Path of the head file. text : string - Name of the text string in the head file. Default is 'head' + Name of the text string in the head file. Default is 'head'. precision : string - 'auto', 'single' or 'double'. Default is 'auto'. + Precision of floating point head data in the value. Accepted + values are 'auto', 'single' or 'double'. Default is 'auto', + which enables automatic detection of precision. verbose : bool - Write information to the screen. Default is False. - - Attributes - ---------- - - Methods - ------- - - See Also - -------- - - Notes - ----- - The HeadFile class provides simple ways to retrieve 2d and 3d - head arrays from a MODFLOW binary head file and time series - arrays for one or more cells. - - The BinaryLayerFile class is built on a record array consisting of - headers, which are record arrays of the modflow header information - (kstp, kper, pertim, totim, text, nrow, ncol, ilay) - and long integers, which are pointers to first bytes of data for - the corresponding data array. + Toggle logging output. Default is False. Examples -------- @@ -624,7 +615,6 @@ class HeadFile(BinaryLayerFile): >>> ddnobj.list_records() >>> rec = ddnobj.get_data(totim=100.) - """ def __init__( @@ -647,6 +637,88 @@ def __init__( ) super().__init__(filename, precision, verbose, kwargs) + def reverse(self, filename: Optional[os.PathLike] = None): + """ + Write a new binary head file with the records in reverse order. + If a new filename is not provided, or if the filename is the same + as the existing filename, the file will be overwritten and data + reloaded from the rewritten/reversed file. + + Parameters + ---------- + + filename : str or PathLike + Path of the new reversed binary file to create. + """ + + filename = ( + Path(filename).expanduser().absolute() + if filename + else self.filename + ) + + # header array formats + dt = np.dtype( + [ + ("kstp", np.int32), + ("kper", np.int32), + ("pertim", np.float64), + ("totim", np.float64), + ("text", "S16"), + ("ncol", np.int32), + ("nrow", np.int32), + ("ilay", np.int32), + ] + ) + + # make sure we have tdis + if self.tdis is None or not any(self.tdis.perioddata.get_data()): + raise ValueError("tdis mu/st be known to reverse head file") + + # extract period data + pd = self.tdis.perioddata.get_data() + + # get maximum period number and total simulation time + kpermx = len(pd) - 1 + tsimtotal = 0.0 + for tpd in pd: + tsimtotal += tpd[0] + + # get total number of records + nrecords = self.recordarray.shape[0] + + # open backward file + with open(filename, "wb") as fbin: + # loop over head file records in reverse order + for idx in range(nrecords - 1, -1, -1): + # load header array + header = self.recordarray[idx].copy() + + # reverse kstp and kper in the header array + (kstp, kper) = (header["kstp"] - 1, header["kper"] - 1) + kstpmx = pd[kper][1] - 1 + kstpb = kstpmx - kstp + kperb = kpermx - kper + (header["kstp"], header["kper"]) = (kstpb + 1, kperb + 1) + + # reverse totim and pertim in the header array + header["totim"] = tsimtotal - header["totim"] + perlen = pd[kper][0] + header["pertim"] = perlen - header["pertim"] + + # write header information + h = np.array(header, dtype=dt) + h.tofile(fbin) + + # load and write data + data = self.get_data(idx=idx)[0][0] + data = np.array(data, dtype=np.float64) + data.tofile(fbin) + + # if we rewrote the original file, reinitialize + if filename == self.filename: + super().__init__(self.filename, self.precision, self.verbose, {}) + class UcnFile(BinaryLayerFile): """ @@ -716,34 +788,190 @@ def __init__( return -class BudgetIndexError(Exception): - pass - - -class CellBudgetFile: +class HeadUFile(BinaryLayerFile): """ - CellBudgetFile Class. + The HeadUFile class provides simple ways to retrieve a list of + head arrays from a MODFLOW-USG binary head file and time series + arrays for one or more cells. Parameters ---------- filename : str or PathLike - Path of the cell budget file + Path of the head file + text : string + Name of the text string in the head file. Default is 'headu'. precision : string - 'single' or 'double'. Default is 'single'. + Precision of the floating point head data in the file. Accepted + values are 'auto', 'single' or 'double'. Default is 'auto', which + enables precision to be automatically detected. verbose : bool - Write information to the screen. Default is False. + Toggle logging output. Default is False. - Attributes - ---------- + Notes + ----- - Methods - ------- + The BinaryLayerFile class is built on a record array consisting of + headers, which are record arrays of the modflow header information + (kstp, kper, pertim, totim, text, nrow, ncol, ilay), and long ints + pointing to the 1st byte of data for the corresponding data arrays. + This class overrides methods in the parent class so that the proper + sized arrays are created: for unstructured grids, nrow and ncol are + the starting and ending node numbers for layer, ilay. - See Also + When the get_data method is called for this class, a list of + one-dimensional arrays will be returned, where each array is the head + array for a layer. If the heads for a layer were not saved, then + None will be returned for that layer. + + Examples -------- - Notes - ----- + >>> import flopy.utils.binaryfile as bf + >>> hdobj = bf.HeadUFile('model.hds') + >>> hdobj.list_records() + >>> usgheads = hdobj.get_data(kstpkper=(1, 50)) + + """ + + def __init__( + self, + filename: Union[str, os.PathLike], + text="headu", + precision="auto", + verbose=False, + **kwargs, + ): + """ + Class constructor + """ + self.text = text.encode() + if precision == "auto": + precision = get_headfile_precision(filename) + if precision == "unknown": + s = f"Error. Precision could not be determined for {filename}" + print(s) + raise Exception() + self.header_dtype = BinaryHeader.set_dtype( + bintype="Head", precision=precision + ) + super().__init__(filename, precision, verbose, kwargs) + + def _get_data_array(self, totim=0.0): + """ + Get a list of 1D arrays for the + specified kstp and kper value or totim value. + + """ + + if totim >= 0.0: + keyindices = np.where(self.recordarray["totim"] == totim)[0] + if len(keyindices) == 0: + msg = f"totim value ({totim}) not found in file..." + raise Exception(msg) + else: + raise Exception("Data not found...") + + # fill a list of 1d arrays with heads from binary file + data = self.nlay * [None] + for idx in keyindices: + ipos = self.iposarray[idx] + ilay = self.recordarray["ilay"][idx] + nstrt = self.recordarray["ncol"][idx] + nend = self.recordarray["nrow"][idx] + npl = nend - nstrt + 1 + if self.verbose: + print(f"Byte position in file: {ipos} for layer {ilay}") + self.file.seek(ipos, 0) + data[ilay - 1] = binaryread(self.file, self.realtype, shape=(npl,)) + return data + + def get_databytes(self, header): + """ + + Parameters + ---------- + header : datafile.Header + header object + + Returns + ------- + databytes : int + size of the data array, in bytes, following the header + + """ + # unstructured head files contain node starting and ending indices + # for each layer + nstrt = np.int64(header["ncol"]) + nend = np.int64(header["nrow"]) + npl = nend - nstrt + 1 + return npl * np.int64(self.realtype(1).nbytes) + + def get_ts(self, idx): + """ + Get a time series from the binary HeadUFile + + Parameters + ---------- + idx : int or list of ints + idx can be nodenumber or it can be a list in the form + [nodenumber, nodenumber, ...]. The nodenumber, + values must be zero based. + + Returns + ---------- + out : numpy array + Array has size (ntimes, ncells + 1). The first column in the + data array will contain time (totim). + + """ + times = self.get_times() + data = self.get_data(totim=times[0]) + layers = len(data) + ncpl = [len(data[l]) for l in range(layers)] + result = [] + + if isinstance(idx, int): + layer, nn = get_lni(ncpl, [idx])[0] + for i, time in enumerate(times): + data = self.get_data(totim=time) + value = data[layer][nn] + result.append([time, value]) + elif isinstance(idx, list) and all(isinstance(x, int) for x in idx): + for i, time in enumerate(times): + data = self.get_data(totim=time) + row = [time] + lni = get_lni(ncpl, idx) + for layer, nn in lni: + value = data[layer][nn] + row += [value] + result.append(row) + else: + raise ValueError("idx must be an integer or a list of integers") + + return np.array(result) + + +class BudgetIndexError(Exception): + pass + + +class CellBudgetFile: + """ + The CellBudgetFile class provides convenient ways to retrieve and + manipulate budget data from a binary cell budget file. A utility + method is also provided to reverse the budget records for particle + tracking simulations in which particles are tracked backwards from + terminating to release locations (e.g., to compute capture zones). + + Parameters + ---------- + filename : str or PathLike + Path of the cell budget file. + precision : string + Precision of floating point budget data in the file. Accepted + values are 'single' or 'double'. Default is 'single'. + verbose : bool + Toggle logging output. Default is False. Examples -------- @@ -796,6 +1024,8 @@ def __init__( if "dis" in kwargs.keys(): self.dis = kwargs.pop("dis") self.modelgrid = self.dis.parent.modelgrid + if "tdis" in kwargs.keys(): + self.tdis = kwargs.pop("tdis") if "sr" in kwargs.keys(): from ..discretization import StructuredGrid, UnstructuredGrid @@ -1969,177 +2199,158 @@ def close(self): Close the file handle """ self.file.close() - return - -class HeadUFile(BinaryLayerFile): - """ - Unstructured MODFLOW-USG HeadUFile Class. - - Parameters - ---------- - filename : string - Name of the concentration file - text : string - Name of the text string in the head file. Default is 'headu' - precision : string - 'auto', 'single' or 'double'. Default is 'auto'. - verbose : bool - Write information to the screen. Default is False. - - Attributes - ---------- - - Methods - ------- - - See Also - -------- - - Notes - ----- - The HeadUFile class provides simple ways to retrieve a list of - head arrays from a MODFLOW-USG binary head file and time series - arrays for one or more cells. - - The BinaryLayerFile class is built on a record array consisting of - headers, which are record arrays of the modflow header information - (kstp, kper, pertim, totim, text, nrow, ncol, ilay) - and long integers, which are pointers to first bytes of data for - the corresponding data array. For unstructured grids, nrow and ncol - are the starting and ending node numbers for layer, ilay. This class - overrides methods in the parent class so that the proper sized arrays - are created. - - When the get_data method is called for this class, a list of - one-dimensional arrays will be returned, where each array is the head - array for a layer. If the heads for a layer were not saved, then - None will be returned for that layer. - - Examples - -------- - - >>> import flopy.utils.binaryfile as bf - >>> hdobj = bf.HeadUFile('model.hds') - >>> hdobj.list_records() - >>> usgheads = hdobj.get_data(kstpkper=(1, 50)) - - - """ - - def __init__( - self, - filename: Union[str, os.PathLike], - text="headu", - precision="auto", - verbose=False, - **kwargs, - ): - """ - Class constructor - """ - self.text = text.encode() - if precision == "auto": - precision = get_headfile_precision(filename) - if precision == "unknown": - s = f"Error. Precision could not be determined for {filename}" - print(s) - raise Exception() - self.header_dtype = BinaryHeader.set_dtype( - bintype="Head", precision=precision - ) - super().__init__(filename, precision, verbose, kwargs) - - def _get_data_array(self, totim=0.0): - """ - Get a list of 1D arrays for the - specified kstp and kper value or totim value. - - """ - - if totim >= 0.0: - keyindices = np.where(self.recordarray["totim"] == totim)[0] - if len(keyindices) == 0: - msg = f"totim value ({totim}) not found in file..." - raise Exception(msg) - else: - raise Exception("Data not found...") - - # fill a list of 1d arrays with heads from binary file - data = self.nlay * [None] - for idx in keyindices: - ipos = self.iposarray[idx] - ilay = self.recordarray["ilay"][idx] - nstrt = self.recordarray["ncol"][idx] - nend = self.recordarray["nrow"][idx] - npl = nend - nstrt + 1 - if self.verbose: - print(f"Byte position in file: {ipos} for layer {ilay}") - self.file.seek(ipos, 0) - data[ilay - 1] = binaryread(self.file, self.realtype, shape=(npl,)) - return data - - def get_databytes(self, header): + def reverse(self, filename: Optional[os.PathLike] = None): """ + Write a binary cell budget file with the records in reverse order. + If a new filename is not provided, or if the filename is the same + as the existing filename, the file will be overwritten and data + reloaded from the rewritten/reversed file. Parameters ---------- - header : datafile.Header - header object - - Returns - ------- - databytes : int - size of the data array, in bytes, following the header - - """ - # unstructured head files contain node starting and ending indices - # for each layer - nstrt = np.int64(header["ncol"]) - nend = np.int64(header["nrow"]) - npl = nend - nstrt + 1 - return npl * np.int64(self.realtype(1).nbytes) - def get_ts(self, idx): + filename : str or PathLike, optional + Path of the new reversed binary cell budget file to create. """ - Get a time series from the binary HeadUFile - Parameters - ---------- - idx : int or list of ints - idx can be nodenumber or it can be a list in the form - [nodenumber, nodenumber, ...]. The nodenumber, - values must be zero based. + filename = ( + Path(filename).expanduser().absolute() + if filename + else self.filename + ) - Returns - ---------- - out : numpy array - Array has size (ntimes, ncells + 1). The first column in the - data array will contain time (totim). + # header array formats + dt1 = np.dtype( + [ + ("kstp", np.int32), + ("kper", np.int32), + ("text", "S16"), + ("ndim1", np.int32), + ("ndim2", np.int32), + ("ndim3", np.int32), + ("imeth", np.int32), + ("delt", np.float64), + ("pertim", np.float64), + ("totim", np.float64), + ] + ) + dt2 = np.dtype( + [ + ("text1id1", "S16"), + ("text1id2", "S16"), + ("text2id1", "S16"), + ("text2id2", "S16"), + ] + ) - """ - times = self.get_times() - data = self.get_data(totim=times[0]) - layers = len(data) - ncpl = [len(data[l]) for l in range(layers)] - result = [] + # make sure we have tdis + if self.tdis is None or not any(self.tdis.perioddata.get_data()): + raise ValueError( + "tdis must be known to reverse a cell budget file" + ) - if isinstance(idx, int): - layer, nn = get_lni(ncpl, [idx])[0] - for i, time in enumerate(times): - data = self.get_data(totim=time) - value = data[layer][nn] - result.append([time, value]) - elif isinstance(idx, list) and all(isinstance(x, int) for x in idx): - for i, time in enumerate(times): - data = self.get_data(totim=time) - row = [time] - lni = get_lni(ncpl, idx) - for layer, nn in lni: - value = data[layer][nn] - row += [value] - result.append(row) - else: - raise ValueError("idx must be an integer or a list of integers") + # extract perioddata + pd = self.tdis.perioddata.get_data() + + # get maximum period number and total simulation time + nper = len(pd) + kpermx = nper - 1 + tsimtotal = 0.0 + for tpd in pd: + tsimtotal += tpd[0] + + # get number of records + nrecords = self.get_nrecords() + + # open backward budget file + with open(filename, "wb") as fbin: + # loop over budget file records in reverse order + for idx in range(nrecords - 1, -1, -1): + # load header array + header = self.recordarray[idx] + + # reverse kstp and kper in the header array + (kstp, kper) = (header["kstp"] - 1, header["kper"] - 1) + kstpmx = pd[kper][1] - 1 + kstpb = kstpmx - kstp + kperb = kpermx - kper + (header["kstp"], header["kper"]) = (kstpb + 1, kperb + 1) + + # reverse totim and pertim in the header array + header["totim"] = tsimtotal - header["totim"] + perlen = pd[kper][0] + header["pertim"] = perlen - header["pertim"] + + # Write main header information to backward budget file + h = header[ + [ + "kstp", + "kper", + "text", + "ncol", + "nrow", + "nlay", + "imeth", + "delt", + "pertim", + "totim", + ] + ] + # Note: much of the code below is based on binary_file_writer.py + h = np.array(h, dtype=dt1) + h.tofile(fbin) + if header["imeth"] == 6: + # Write additional header information to the backward budget file + h = header[ + [ + "modelnam", + "paknam", + "modelnam2", + "paknam2", + ] + ] + h = np.array(h, dtype=dt2) + h.tofile(fbin) + # Load data + data = self.get_data(idx)[0] + data = np.array(data) + # Negate flows + data["q"] = -data["q"] + # Write ndat (number of floating point columns) + colnames = data.dtype.names + ndat = len(colnames) - 2 + dt = np.dtype([("ndat", np.int32)]) + h = np.array([(ndat,)], dtype=dt) + h.tofile(fbin) + # Write auxiliary column names + naux = ndat - 1 + if naux > 0: + auxtxt = [ + "{:16}".format(colname) for colname in colnames[3:] + ] + auxtxt = tuple(auxtxt) + dt = np.dtype( + [(colname, "S16") for colname in colnames[3:]] + ) + h = np.array(auxtxt, dtype=dt) + h.tofile(fbin) + # Write nlist + nlist = data.shape[0] + dt = np.dtype([("nlist", np.int32)]) + h = np.array([(nlist,)], dtype=dt) + h.tofile(fbin) + elif header["imeth"] == 1: + # Load data + data = self.get_data(idx)[0][0][0] + data = np.array(data, dtype=np.float64) + # Negate flows + data = -data + else: + raise ValueError("not expecting imeth " + header["imeth"]) + # Write data + data.tofile(fbin) - return np.array(result) + # if we rewrote the original file, reinitialize + if filename == self.filename: + self.__init__(self.filename, self.precision, self.verbose) diff --git a/flopy/utils/datafile.py b/flopy/utils/datafile.py index 87fb4d9d9..81b84901c 100644 --- a/flopy/utils/datafile.py +++ b/flopy/utils/datafile.py @@ -196,6 +196,8 @@ def __init__( if "dis" in kwargs.keys(): self.dis = kwargs.pop("dis") self.mg = self.dis.parent.modelgrid + if "tdis" in kwargs.keys(): + self.tdis = kwargs.pop("tdis") if "modelgrid" in kwargs.keys(): self.mg = kwargs.pop("modelgrid") if len(kwargs.keys()) > 0: @@ -422,6 +424,11 @@ def list_records(self): print(header) return + def get_nrecords(self): + if isinstance(self.recordarray, np.recarray): + return self.recordarray.shape[0] + return 0 + def _get_data_array(self, totim=0): """ Get the three dimensional data array for the