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

@@ -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)

flopy/utils/binaryfile.py → 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[

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}")