Skip to content

Commit

Permalink
feat(HeadFile): add reverse() method
Browse files Browse the repository at this point in the history
  • Loading branch information
@wpbonelli
wpbonelli committed Jun 20, 2023
1 parent 10d3659 commit fe789f4
Show file tree
Hide file tree
Showing 2 changed files with 141 additions and 29 deletions.
64 changes: 53 additions & 11 deletions autotest/test_binaryfile.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
CellBudgetFile,
HeadFile,
HeadUFile,
UcnFile,
Util2d,
)
from flopy.utils.binaryfile import (
Expand Down Expand Up @@ -563,13 +564,14 @@ def test_cellbudgetfile_reverse_mf2005(example_data_path, function_tmpdir):


def test_cellbudgetfile_reverse_mf6(example_data_path, function_tmpdir):
cbc_name = "flow_adj"
mf6_model_path = example_data_path / "mf6" / "test006_gwf3"
cbc_fname = mf6_model_path / "expected_output" / f"{cbc_name}.cbc"
f = CellBudgetFile(cbc_fname)
model_path = example_data_path / "mf6" / "test006_gwf3"
file_stem = "flow_adj"
file_path = model_path / "expected_output" / f"{file_stem}.cbc"
f = CellBudgetFile(file_path)
assert isinstance(f, CellBudgetFile)

rf_name = f"{cbc_name}_rev.cbc"
# reverse the file
rf_name = f"{file_stem}_rev.cbc"
f.reverse(function_tmpdir / rf_name)
rf = CellBudgetFile(function_tmpdir / rf_name)
assert isinstance(rf, CellBudgetFile)
Expand All @@ -578,13 +580,11 @@ def test_cellbudgetfile_reverse_mf6(example_data_path, function_tmpdir):
nrecords = f.get_nrecords()
assert nrecords == rf.get_nrecords()

# check that the data is reversed
# check that the data are reversed
for idx in range(nrecords - 1, -1, -1):
# check headers
# todo: figure out disagreement
# f_header = f.recordarray[idx]
# rf_header = rf.recordarray[idx]
# assert np.array_equal(f_header, rf_header)
# todo check headers
f_header = f.recordarray[idx]
rf_header = rf.recordarray[idx]

# check data
f_data = f.get_data(idx=idx)[0]
Expand All @@ -594,7 +594,49 @@ def test_cellbudgetfile_reverse_mf6(example_data_path, function_tmpdir):
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])


def test_headfile_reverse_mf6(example_data_path, function_tmpdir):
model_path = example_data_path / "mf6" / "test006_gwf3"
file_stem = "flow_adj"
file_path = model_path / "expected_output" / f"{file_stem}.hds"
f = HeadFile(file_path)
assert isinstance(f, HeadFile)

# reverse the file
rf_name = f"{file_stem}_rev.hds"
f.reverse(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 that the data are reversed
for idx in range(len(f_data) - 1, -1, -1):
# todo check headers
# todo: figure out disagreement
f_header = f.recordarray[idx]
rf_header = rf.recordarray[idx]

# check data
f_data = f.get_data(idx=idx)[0]
rf_data = rf.get_data(idx=len(f_data) - 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])
106 changes: 88 additions & 18 deletions flopy/utils/binaryfile.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -636,6 +636,73 @@ def __init__(
)
super().__init__(filename, precision, verbose, kwargs)

def reverse(self, filename: os.PathLike):
"""
Write a new binary file with the records in reverse order.
Parameters
----------
filename : str or PathLike
Path of the new reversed binary file to create.
"""

# 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),
]
)

# get maximum period number and total simulation time
kpermx = max(self.recordarray["kper"])
tsimtotal = max(self.get_times())

# get total number of records
nrecords = len(self.get_alldata())

# reverse kstpkper and times
kstpkper = list(reversed(self.kstpkper.copy()))
times = list(reversed(self.times.copy()))

# 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"], header["kper"])
kstpmx = next(iter([kk for kk in kstpkper if kk[1] == kper]))[
0
]
(header["kstp"], header["kper"]) = (
kstpmx - kstp + 1,
kpermx - kper + 1,
)

# reverse totim and pertim in the header array
header["totim"] = tsimtotal - header["totim"]
perlen = times[kper - 1] - times[kper - 2]
header["pertim"] = perlen - header["pertim"]

# write header information to the backward head file
h = np.array(header, dtype=dt)
h.tofile(fbin)
# load data
data = self.get_data(idx=idx)[0][0]
data = np.array(data, dtype=np.float64)
# write data
data.tofile(fbin)


class UcnFile(BinaryLayerFile):
"""
Expand Down Expand Up @@ -1964,13 +2031,6 @@ def reverse(self, filename: os.PathLike):
Path of the new reversed binary cell budget file to create.
"""

# get maximum period number and total simulation time
kpermx = self.nper - 1
tsimtotal = max(self.times)

# get number of records
nrecords = self.get_nrecords()

# header array formats
dt1 = np.dtype(
[
Expand All @@ -1995,29 +2055,39 @@ def reverse(self, filename: os.PathLike):
]
)

# reverse kstpker
kstpker = list(reversed(self.kstpkper))
times = list(reversed(self.times))
# get maximum period number and total simulation time
kpermx = self.nper - 1
tsimtotal = max(self.times)

# get number of records
nrecords = self.get_nrecords()

# Open backward budget file
# reverse kstpkper
kstpkper = list(reversed(self.kstpkper.copy()))
times = list(reversed(self.times.copy()))

# open backward budget file
with open(filename, "wb") as fbin:
# Loop over budget file records in reverse order
# loop over budget file records in reverse order
for idx in range(nrecords - 1, -1, -1):
# Load header array
# load header array
header = self.recordarray[idx]
# Replace kstp and kper in the header array with their
# backward countparts

# reverse kstp and kper in the header array
(kstp, kper) = (header["kstp"], header["kper"])
kstpmx = next(iter([kk for kk in kstpker if kk[1] == kper]))[0]
kstpmx = next(iter([kk for kk in kstpkper if kk[1] == kper]))[
0
]
(header["kstp"], header["kper"]) = (
kstpmx - kstp,
kpermx - kper,
)
# Replace totim and pertim in the header array with their
# backward counterparts

# reverse totim and pertim in the header array
header["totim"] = tsimtotal - header["totim"]
perlen = times[kper - 1] - times[kper - 2]
header["pertim"] = perlen - header["pertim"]

# Write main header information to backward budget file
h = header[
[
Expand Down

0 comments on commit fe789f4

Please sign in to comment.