From 4eac63176f1328a22a55e1cf131db55b8ca929a8 Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Tue, 26 Nov 2024 21:55:38 -0500 Subject: [PATCH] feat(binaryfile): add head/budget file reversal script (#2383) Close #2382. Also remove the need for tdis in the reversal methods by inferring period length etc from data headers. --- autotest/test_binaryfile.py | 14 +-- .../{binaryfile.py => binaryfile/__init__.py} | 96 +++++++++---------- flopy/utils/binaryfile/reverse.py | 31 ++++++ 3 files changed, 86 insertions(+), 55 deletions(-) rename flopy/utils/{binaryfile.py => binaryfile/__init__.py} (97%) create mode 100644 flopy/utils/binaryfile/reverse.py diff --git a/autotest/test_binaryfile.py b/autotest/test_binaryfile.py index ce993439c..dd5b8161b 100644 --- a/autotest/test_binaryfile.py +++ b/autotest/test_binaryfile.py @@ -502,7 +502,7 @@ def test_binaryfile_reverse_mf6_dis(function_tmpdir): assert success, pformat(buff) # reverse head file in place and check reversal - head_file = flopy.utils.HeadFile(function_tmpdir / head_file, tdis=tdis) + head_file = flopy.utils.HeadFile(function_tmpdir / head_file) heads = head_file.get_alldata() assert heads.shape == (nper, 2, 10, 10) head_file.reverse() @@ -511,9 +511,9 @@ def test_binaryfile_reverse_mf6_dis(function_tmpdir): # reverse budget and write to separate file budget_file_rev_path = function_tmpdir / f"{budget_file}_rev" - budget_file = flopy.utils.CellBudgetFile(function_tmpdir / budget_file, tdis=tdis) + budget_file = flopy.utils.CellBudgetFile(function_tmpdir / budget_file) budget_file.reverse(budget_file_rev_path) - budget_file_rev = flopy.utils.CellBudgetFile(budget_file_rev_path, tdis=tdis) + budget_file_rev = flopy.utils.CellBudgetFile(budget_file_rev_path) for kper in range(nper): assert np.allclose(heads[kper], heads_rev[-kper + 1]) @@ -554,7 +554,7 @@ def test_binaryfile_reverse_mf6_disv(function_tmpdir): assert success, pformat(buff) # reverse head file in place and check reversal - head_file = flopy.utils.HeadFile(function_tmpdir / head_file, tdis=tdis) + head_file = flopy.utils.HeadFile(function_tmpdir / head_file) heads = head_file.get_alldata() assert heads.shape == (nper, 2, 1, 100) head_file.reverse() @@ -563,7 +563,7 @@ def test_binaryfile_reverse_mf6_disv(function_tmpdir): # reverse budget and write to separate file budget_file_rev_path = function_tmpdir / f"{budget_file}_rev" - budget_file = flopy.utils.CellBudgetFile(function_tmpdir / budget_file, tdis=tdis) + budget_file = flopy.utils.CellBudgetFile(function_tmpdir / budget_file) budget_file.reverse(budget_file_rev_path) budget_file_rev = flopy.utils.CellBudgetFile(budget_file_rev_path, tdis=tdis) @@ -595,7 +595,7 @@ def test_binaryfile_reverse_mf6_disu(example_data_path, function_tmpdir): # reverse and write to a separate file head_file_rev_path = function_tmpdir / "flow_rev.hds" head_file.reverse(filename=head_file_rev_path) - head_file_rev = HeadFile(head_file_rev_path, tdis=tdis) + head_file_rev = HeadFile(head_file_rev_path) # load budget file file_path = function_tmpdir / "flow.cbc" @@ -604,7 +604,7 @@ def test_binaryfile_reverse_mf6_disu(example_data_path, function_tmpdir): # reverse and write to a separate file budget_file_rev_path = function_tmpdir / "flow_rev.cbc" budget_file.reverse(filename=budget_file_rev_path) - budget_file_rev = CellBudgetFile(budget_file_rev_path, tdis=tdis) + budget_file_rev = CellBudgetFile(budget_file_rev_path) # check that data from both files have the same shape assert head_file.get_alldata().shape == (nper, 1, 1, 121) diff --git a/flopy/utils/binaryfile.py b/flopy/utils/binaryfile/__init__.py similarity index 97% rename from flopy/utils/binaryfile.py rename to flopy/utils/binaryfile/__init__.py index 4aa31260d..194c1140f 100644 --- a/flopy/utils/binaryfile.py +++ b/flopy/utils/binaryfile/__init__.py @@ -19,8 +19,8 @@ import numpy as np import pandas as pd -from ..utils.datafile import Header, LayerFile -from .gridutil import get_lni +from ..datafile import Header, LayerFile +from ..gridutil import get_lni HEAD_TEXT = " HEAD" @@ -663,38 +663,28 @@ def get_max_kper_kstp_tsim(): kstp[header["kper"]] = 0 return kper, kstp, tsim - # get max period and time from the head file maxkper, maxkstp, maxtsim = get_max_kper_kstp_tsim() - # if we have tdis, get max period number and simulation time from it - tdis_maxkper, tdis_maxtsim = None, None - if self.tdis is not None: - pd = self.tdis.perioddata.get_data() - if any(pd): - tdis_maxkper = len(pd) - 1 - tdis_maxtsim = sum([p[0] for p in pd]) - # if we have both, check them against each other - if tdis_maxkper is not None: - assert maxkper == tdis_maxkper, ( - f"Max stress period in binary head file ({maxkper}) != " - f"max stress period in provided tdis ({tdis_maxkper})" - ) - assert maxtsim == tdis_maxtsim, ( - f"Max simulation time in binary head file ({maxtsim}) != " - f"max simulation time in provided tdis ({tdis_maxtsim})" - ) + prev_kper = None + perlen = None def reverse_header(header): """Reverse period, step and time fields in the record header""" + nonlocal prev_kper + nonlocal perlen + # reverse kstp and kper headers kstp = header["kstp"] - 1 kper = header["kper"] - 1 header["kstp"] = maxkstp[kper] - kstp + 1 header["kper"] = maxkper - kper + 1 + if kper != prev_kper: + perlen = header["pertim"] + prev_kper = kper + # reverse totim and pertim headers header["totim"] = maxtsim - header["totim"] - perlen = pd[kper][0] header["pertim"] = perlen - header["pertim"] return header @@ -1022,7 +1012,6 @@ def __init__( self.paknamlist_from = [] self.paknamlist_to = [] self.compact = True # compact budget file flag - self.dis = None self.modelgrid = None if "model" in kwargs.keys(): @@ -2237,24 +2226,46 @@ def reverse(self, filename: Optional[os.PathLike] = None): ] ) - # 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") + nrecords = len(self) + target = filename + + def get_max_kper_kstp_tsim(): + header = self.recordarray[-1] + kper = header["kper"] - 1 + tsim = header["totim"] + kstp = {0: 0} + for i in range(len(self) - 1, -1, -1): + header = self.recordarray[i] + if header["kper"] in kstp and header["kstp"] > kstp[header["kper"]]: + kstp[header["kper"]] += 1 + else: + kstp[header["kper"]] = 0 + return kper, kstp, tsim + + maxkper, maxkstp, maxtsim = get_max_kper_kstp_tsim() + prev_kper = None + perlen = None - # extract perioddata - pd = self.tdis.perioddata.get_data() + def reverse_header(header): + """Reverse period, step and time fields in the record header""" - # get maximum period number and total simulation time - nper = len(pd) - kpermx = nper - 1 - tsimtotal = 0.0 - for tpd in pd: - tsimtotal += tpd[0] + nonlocal prev_kper + nonlocal perlen - # get number of records - nrecords = len(self) + # reverse kstp and kper headers + kstp = header["kstp"] - 1 + kper = header["kper"] - 1 + header["kstp"] = maxkstp[kper] - kstp + 1 + header["kper"] = maxkper - kper + 1 - target = filename + if kper != prev_kper: + perlen = header["pertim"] - 1 + prev_kper = kper + + # reverse totim and pertim headers + header["totim"] = maxtsim - header["totim"] + header["pertim"] = perlen - header["pertim"] + return header # if rewriting the same file, write # temp file then copy it into place @@ -2269,18 +2280,7 @@ def reverse(self, filename: Optional[os.PathLike] = None): 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"] + header = reverse_header(header) # Write main header information to backward budget file h = header[ diff --git a/flopy/utils/binaryfile/reverse.py b/flopy/utils/binaryfile/reverse.py new file mode 100644 index 000000000..f69e00f5d --- /dev/null +++ b/flopy/utils/binaryfile/reverse.py @@ -0,0 +1,31 @@ +import argparse +from pathlib import Path + +from flopy.utils.binaryfile import CellBudgetFile, HeadFile + +if __name__ == "__main__": + """Reverse head or budget files.""" + + parser = argparse.ArgumentParser(description="Reverse head or budget files.") + parser.add_argument( + "--infile", + "-i", + type=str, + help="Input file.", + ) + parser.add_argument( + "--outfile", + "-o", + type=str, + help="Output file.", + ) + args = parser.parse_args() + infile = Path(args.infile) + outfile = Path(args.outfile) + suffix = infile.suffix.lower() + if suffix in [".hds", ".hed"]: + HeadFile(infile).reverse(outfile) + elif suffix in [".bud", ".cbc"]: + CellBudgetFile(infile).reverse(outfile) + else: + raise ValueError(f"Unrecognized file suffix: {suffix}")