"""Read MODFLOW 6 binary grid (.grb) and budget (.cbc) files into a FlowStep.
Node numbering (checked on MF6 6.x nightly with IDOMAIN = 0 cells on DIS and DISU):
- The .grb file stores IA/JA in *user* node numbering for DIS, DISV and DISU. Cells removed by
IDOMAIN (<= 0) keep an empty row (no diagonal either) and no other row points at them.
- FLOW-JA-FACE has NJA entries aligned with those IA/JA arrays, so it needs no remapping.
- List budget records (IMETH = 6) give the 1-based user node in ID1 and the bound number in ID2.
- Array budget records (IMETH = 1, e.g. STO-SS) hold one value per user node.
So every array is in user numbering; cells with IDOMAIN <= 0 are flagged inactive.
Sign convention: every budget term is positive INTO the cell, including storage (positive =
release from storage). MF6 writes the per-cell residual (sum of all flows into the cell) in the
diagonal of FLOW-JA-FACE; it is not a connection and is ignored here. DATA-* records (DATA-SPDIS,
DATA-SAT) are not flows and are skipped.
Cell area: DELR x DELC (DIS); vertex polygons (DISV, and DISU when vertices were written);
otherwise AREA from the DISU input file read with FloPy, or None if that is not available.
Saturated thickness (for travel times): cell top and bottom and ICELLTYPE come from the .grb
file (DIS/DISV: the top of layer k > 0 is the bottom of layer k - 1). Confined cells
(ICELLTYPE = 0) are fully saturated. Convertible cells use DATA-SAT x (top - bot) when the step
saved DATA-SAT (NPF SAVE_SATURATION); otherwise, if a head file is open, clip(head, bot, top) -
bot, with dry cells (head <= bot, or |head| >= 1e29 for HDRY/HNOFLO) at 0; otherwise None.
Water mover (checked on test model 7, DRN -> MVR -> SFR): a provider package with MOVER writes
two disjoint aquifer-side records, e.g. DRN (the part of the drain outflow that leaves the model)
and DRN-TO-MVR (the part handed to the mover). Both are outflows from the cell and together they
equal the drain's head-dependent flow, so both are kept as sink terms; the cell budget closes
against MF6's own residual with both, and would not close without the -TO-MVR record. Where the
mover sends the water (here into the stream, which may leak it back into the aquifer through
the SFR record) is not followed in a single-model step: gwsource.network routes it when the
receiving package is read as package nodes (read_network / open_simulation with package_nodes=
and the MVR budget file). Records of one package share its name (paknam2); the second and
later records of a package in a step are named '<package>:<text>' (e.g. 'DRN:DRN-TO-MVR') and
keep the package as their parent.
GWF-GWF exchanges (checked with models split by FloPy's Mf6Splitter): with SAVE_FLOWS in the
exchange file, each model's budget has one FLOW-JA-FACE record with IMETH = 6 per exchange,
PAKNAM2 = the exchange name (e.g. 'GWF-GWF_1'), MODELNAM2 = the other model, ID1 = the cell
here, ID2 = the cell in the other model, q positive into the cell here; the other model's record
is its exact negative. Read for one model alone, these are boundary terms: a package with text
'GWF-GWF', the exchange name, partner = the other model and bound = the cell over there.
gwsource.network turns them into edges when both models are read together.
IFLOWFACE (checked on MF6 6.9 dev, RCH with AUXILIARY IFLOWFACE): a package's auxiliary
variables are saved as extra columns of its list budget record (FloPy names the column by the
AUXTXT, here 'IFLOWFACE'). MF6 PRT reads IFLOWFACE from the same budget file, so this column is
what PRT uses. FlowStep.term_iflowface holds it per term (0 for records without it).
Boundnames are not in the budget file. ID2 of a list record is the bound's position (1-based)
in the package's list for the current stress period, so read_boundnames() loads the list input
of every package with BOUNDNAMES through FloPy, per stress period (a period without a PERIOD
block reuses the last list, as in MF6). SFR, LAK, MAW and UZF boundnames come from their
PACKAGEDATA (ID2 = reach, lake, well or UZF cell; also the bound of their package-node terms).
"""
from dataclasses import dataclass
from pathlib import Path
import flopy
import numpy as np
from .flowstep import FlowStep, Package
COMMENT_CHARS = ("#", "!", "//")
DRY_HEAD = 1.0e29
EXCHANGE_TEXT = "GWF-GWF"
PACKAGEDATA_TYPES = ("sfr", "lak", "maw", "uzf") # boundnames by reach, lake, well or UZF cell
IFLOWFACE = "IFLOWFACE"
[docs]
@dataclass
class ModelFiles:
"""Paths of the files of one GWF model that a FlowStep is read from."""
model: str
grb_file: Path
budget_file: Path
dis_file: Path
dis_type: str
nam_file: Path | None = None
tdis_file: Path | None = None
head_file: Path | None = None
[docs]
@dataclass
class GridData:
"""Connectivity and geometry read from a binary grid file."""
grid_type: str
ia: np.ndarray
ja: np.ndarray
idomain: np.ndarray
shape: tuple[int, ...]
area: np.ndarray | None
top: np.ndarray | None = None
bot: np.ndarray | None = None
icelltype: np.ndarray | None = None
# ---------------------------------------------------------------- name-file parsing
def _block_lines(path: Path, block: str) -> list[list[str]]:
"""Return the tokenised lines inside 'BEGIN block' ... 'END block' of an MF6 input file."""
lines, inside = [], False
for raw in Path(path).read_text(encoding="utf-8").splitlines():
tokens = raw.split()
if not tokens or tokens[0].startswith(COMMENT_CHARS):
continue
key = tokens[0].upper()
if key == "BEGIN" and len(tokens) > 1 and tokens[1].upper() == block.upper():
inside = True
elif key == "END" and inside:
return lines
elif inside:
lines.append(tokens)
return lines
def _find_model_entry(ws: Path, model: str | None) -> tuple[str, Path]:
"""Return (model name, name file) of a GWF model listed in mfsim.nam."""
entries = [t for t in _block_lines(ws / "mfsim.nam", "MODELS") if t[0].upper() == "GWF6"]
assert entries, f"no GWF6 model in {ws / 'mfsim.nam'}"
if model is None:
assert len(entries) == 1, "several GWF models; pass model= to choose one"
return entries[0][2], ws / entries[0][1]
match = [t for t in entries if t[2].lower() == model.lower()]
assert match, f"model {model!r} not found in mfsim.nam"
return match[0][2], ws / match[0][1]
def _fileout(path: Path, keyword: str) -> str | None:
"""Return the file name of a '<keyword> FILEOUT <name>' option in an input file, if present."""
for tokens in _block_lines(path, "OPTIONS"):
if len(tokens) >= 3 and tokens[0].upper() == keyword and tokens[1].upper() == "FILEOUT":
return tokens[2]
return None
[docs]
def find_model_files(ws: str | Path, model: str | None = None) -> ModelFiles:
"""Locate the binary grid and budget files of a GWF model from its name and OC files."""
ws = Path(ws)
name, nam_file = _find_model_entry(ws, model)
packages = {t[0].upper(): t[1] for t in _block_lines(nam_file, "PACKAGES")}
dis_type = next(k for k in ("DIS6", "DISV6", "DISU6") if k in packages)
dis_file = ws / packages[dis_type]
assert "OC6" in packages, f"no OC package in {nam_file}"
budget = _fileout(ws / packages["OC6"], "BUDGET")
assert budget is not None, "OC does not save a budget file (BUDGET FILEOUT)"
grb = _fileout(dis_file, "GRB6") or f"{packages[dis_type]}.grb"
timing = {t[0].upper(): t[1] for t in _block_lines(ws / "mfsim.nam", "TIMING")}
tdis = ws / timing["TDIS6"] if "TDIS6" in timing else None
head = _fileout(ws / packages["OC6"], "HEAD")
head_file = ws / head if head is not None else None
return ModelFiles(
name, ws / grb, ws / budget, dis_file, dis_type[:-1], nam_file, tdis, head_file
)
def _option_value(path: Path | None, keyword: str) -> str | None:
"""Return the value of '<keyword> <value>' in an input file's OPTIONS block, if present."""
if path is None or not Path(path).is_file():
return None
for tokens in _block_lines(path, "OPTIONS"):
if tokens[0].upper() == keyword and len(tokens) > 1:
return tokens[1].lower()
return None
[docs]
def read_units(files: ModelFiles) -> dict[str, str]:
"""Return {'length': ..., 'time': ...} from DIS LENGTH_UNITS and TDIS TIME_UNITS."""
return {
"length": _option_value(files.dis_file, "LENGTH_UNITS") or "unknown",
"time": _option_value(files.tdis_file, "TIME_UNITS") or "unknown",
}
# ---------------------------------------------------------------- binary grid
[docs]
def polygon_areas(vertices: np.ndarray, iavert: np.ndarray, javert: np.ndarray) -> np.ndarray:
"""Return plan areas of polygons given vertex xy and 0-based CSR vertex lists (shoelace)."""
x, y = vertices[javert, 0], vertices[javert, 1]
nxt = np.arange(1, len(javert) + 1)
nxt[iavert[1:] - 1] = iavert[:-1] # last vertex of each polygon wraps to its first
cross = x * y[nxt] - x[nxt] * y
return 0.5 * np.abs(np.add.reduceat(cross, iavert[:-1]))
def _vertex_areas(grb: flopy.mf6.utils.MfGrdFile) -> np.ndarray | None:
"""Return cell-polygon areas from the grid file's vertices, or None if it has none."""
if grb.iavert is None or grb.verts is None:
return None
verts = np.asarray(grb.verts, dtype=float)[:, -2:] # FloPy prefixes DISU vertices with an id
return polygon_areas(verts, np.asarray(grb.iavert), np.asarray(grb.javert))
def _grid_area(grb: flopy.mf6.utils.MfGrdFile) -> np.ndarray | None:
"""Return the plan area of every user node, or None if the grid file has no geometry."""
if grb.grid_type == "DIS":
cell = np.outer(grb.delc, grb.delr).ravel()
return np.tile(cell, grb.nlay)
areas = _vertex_areas(grb)
if areas is not None and grb.grid_type == "DISV":
return np.tile(areas, grb.nlay)
return areas
def _grid_shape(grb: flopy.mf6.utils.MfGrdFile) -> tuple[int, ...]:
"""Return the array shape of the grid: (nlay, nrow, ncol), (nlay, ncpl) or (nodes,)."""
if grb.grid_type == "DIS":
return (int(grb.nlay), int(grb.nrow), int(grb.ncol))
if grb.grid_type == "DISV":
return (int(grb.nlay), int(grb.ncpl))
return (int(grb.nodes),)
def _cell_top_bot(grb: flopy.mf6.utils.MfGrdFile) -> tuple[np.ndarray, np.ndarray]:
"""Return the top and bottom elevation of every user node."""
bot = np.asarray(grb.bot, dtype=float).ravel()
top = np.asarray(grb.top, dtype=float).ravel()
if top.size < bot.size: # DIS/DISV: model top, then each layer's top is the bottom above
top = np.concatenate([top, bot[: bot.size - top.size]])
return top, bot
[docs]
def read_grid(grb_file: str | Path) -> GridData:
"""Read connectivity, IDOMAIN, shape, geometry and ICELLTYPE from a binary grid file."""
grb = flopy.mf6.utils.MfGrdFile(str(grb_file), verbose=False)
assert grb.grid_type in ("DIS", "DISV", "DISU"), f"unsupported grid type {grb.grid_type}"
idomain = np.ravel(grb.idomain) if grb.idomain is not None else np.ones(grb.nodes, int)
top, bot = _cell_top_bot(grb)
icelltype = grb._datadict.get("ICELLTYPE") # no public attribute in FloPy
return GridData(
grid_type=grb.grid_type,
ia=np.asarray(grb.ia, dtype=np.int64),
ja=np.asarray(grb.ja, dtype=np.int64),
idomain=np.asarray(idomain, dtype=int),
shape=_grid_shape(grb),
area=_grid_area(grb),
top=top,
bot=bot,
icelltype=None if icelltype is None else np.asarray(icelltype, dtype=int).ravel(),
)
[docs]
def saturated_thickness(
top: np.ndarray,
bot: np.ndarray,
icelltype: np.ndarray | None,
active: np.ndarray,
saturation: np.ndarray | None = None,
head: np.ndarray | None = None,
) -> np.ndarray:
"""Return per-cell saturated thickness from DATA-SAT (preferred) or heads; 0 if inactive.
Cells with ICELLTYPE = 0 are fully saturated; without ICELLTYPE every cell is convertible.
"""
assert saturation is not None or head is not None, "need saturation or head"
full = np.maximum(top - bot, 0.0)
if saturation is not None:
partial = np.clip(saturation, 0.0, 1.0) * full
else:
dry = (np.abs(head) >= DRY_HEAD) | (head <= bot)
partial = np.where(dry, 0.0, np.clip(head, bot, top) - bot)
confined = np.zeros(top.size, bool) if icelltype is None else icelltype == 0
return np.where(active, np.where(confined, full, partial), 0.0)
[docs]
def read_disu_area(ws: str | Path, model: str) -> np.ndarray:
"""Return AREA from a model's DISU input, read with FloPy."""
sim = flopy.mf6.MFSimulation.load(sim_ws=str(ws), load_only=["disu"], verbosity_level=0)
gwf = sim.get_model(model)
return np.asarray(gwf.disu.area.array, dtype=float).ravel()
# ---------------------------------------------------------------- boundnames
def _name_file_packages(nam_file: Path) -> list[tuple[str, Path, str]]:
"""Return (ftype without '6', file, package name) of each entry of a model name file."""
entries = []
counts: dict[str, int] = {}
for tokens in _block_lines(nam_file, "PACKAGES"):
ftype = tokens[0].upper().removesuffix("6")
counts[ftype] = counts.get(ftype, 0) + 1
pname = tokens[2] if len(tokens) > 2 else f"{ftype}-{counts[ftype]}"
entries.append((ftype, nam_file.parent / tokens[1], pname.upper()))
return entries
def _has_boundnames(path: Path) -> bool:
"""Return True if an input file sets the BOUNDNAMES option."""
return any(t[0].upper() == "BOUNDNAMES" for t in _block_lines(path, "OPTIONS"))
def _names_array(records) -> np.ndarray:
"""Return the boundname column of a FloPy recarray as an object array ('' if missing)."""
if records is None or "boundname" not in (records.dtype.names or ()):
return np.zeros(0, dtype=object) if records is None else np.full(len(records), "", object)
return np.array([str(b) if b is not None else "" for b in records["boundname"]], dtype=object)
def _package_boundnames(package) -> dict[int, np.ndarray]:
"""Return {first kper: names} for a FloPy package with list data or advanced package data."""
if package.package_type in PACKAGEDATA_TYPES:
return {0: _names_array(package.packagedata.get_data())}
spd = getattr(package, "stress_period_data", None)
assert spd is not None, (
f"boundnames of {package.package_type.upper()} packages are not supported yet"
)
return {int(k): _names_array(v) for k, v in (spd.get_data() or {}).items()}
[docs]
def read_boundnames(files: ModelFiles, ws: str | Path) -> BoundnameTable:
"""Load, through FloPy, the boundnames of every package of a model that has BOUNDNAMES."""
entries = [e for e in _name_file_packages(files.nam_file) if _has_boundnames(e[1])]
if not entries:
return BoundnameTable({})
load_only = sorted({ftype.lower() for ftype, _, _ in entries})
sim = flopy.mf6.MFSimulation.load(sim_ws=str(ws), load_only=load_only, verbosity_level=0)
gwf = sim.get_model(files.model)
periods = {}
for _, _, pname in entries:
package = gwf.get_package(pname.lower())
assert package is not None, f"package {pname} not loaded by FloPy"
periods[pname] = _package_boundnames(package)
return BoundnameTable(periods)
# ---------------------------------------------------------------- budget file
[docs]
def available_steps(budget_file: str | Path) -> list[tuple[int, int]]:
"""Return the 0-based (kstp, kper) of every time step saved in a budget file."""
cbc = flopy.utils.CellBudgetFile(str(budget_file))
return [(int(k), int(p)) for k, p in cbc.get_kstpkper()]
def _step_record_indices(cbc: flopy.utils.CellBudgetFile, kstpkper: tuple[int, int]) -> list[int]:
"""Return the record indices of one (0-based) time step."""
kstp, kper = kstpkper
rec = cbc.recordarray
idx = np.flatnonzero((rec["kstp"] == kstp + 1) & (rec["kper"] == kper + 1))
assert idx.size, f"time step {kstpkper} not in budget file"
return [int(i) for i in idx]
def _decode(value) -> str:
"""Decode and strip a budget-header text field."""
return value.decode().strip() if isinstance(value, bytes) else str(value).strip()
def _array_term(data: np.ndarray, active: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
"""Return (0-based cell, flow) for an IMETH=1 array record, one entry per active cell."""
flow = np.asarray(data, dtype=float).ravel()
assert flow.size == active.size, "array budget record does not match the number of nodes"
cells = np.flatnonzero(active)
return cells, flow[cells]
def _saturation_record(data: np.ndarray, n_nodes: int) -> np.ndarray:
"""Return per-node saturation from a DATA-SAT list record (0 where not listed)."""
saturation = np.zeros(n_nodes)
saturation[np.asarray(data["node"]) - 1] = np.asarray(data["sat"], dtype=float)
return saturation
def _record_package(parent: str, text: str, used: set[str]) -> Package:
"""Return the Package of a record; a second record of the same package gets '<pkg>:<text>'."""
name = parent if parent not in used else f"{parent}:{text}"
assert name not in used, f"duplicate budget record {name}"
used.add(name)
return Package(name=name, text=text, parent=parent if name != parent else "")
def _read_terms(cbc, indices: list[int], active: np.ndarray) -> dict:
"""Collect external term entries and the FLOW-JA-FACE array from one step's records."""
out = {"cell": [], "flow": [], "package": [], "bound": [], "packages": [], "flowja": None}
out["iflowface"] = []
out["saturation"] = None
used: set[str] = set()
for idx in indices:
header = cbc.recordarray[idx]
text, imeth = _decode(header["text"]), int(header["imeth"])
data = cbc.get_record(idx, full3D=False)
if text == "FLOW-JA-FACE" and imeth == 1:
out["flowja"] = np.asarray(data, dtype=float).ravel()
continue
if text == "FLOW-JA-FACE": # GWF-GWF exchange: ID1 = node here, ID2 = node in partner
package = Package(
name=_decode(header["paknam2"]),
text=EXCHANGE_TEXT,
partner=_decode(header["modelnam2"]),
)
assert package.name not in used, f"duplicate exchange record {package.name}"
used.add(package.name)
_append_terms(out, package, data["node"] - 1, data["q"], data["node2"])
continue
if text == "DATA-SAT":
out["saturation"] = _saturation_record(data, active.size)
continue
if text.startswith("DATA-"):
continue
if imeth == 6:
parent = _decode(header["paknam2"]) or text
cells, flow, bound = data["node"] - 1, data["q"], data["node2"]
iflowface = _aux_column(data, IFLOWFACE)
else:
parent = text
cells, flow = _array_term(data, active)
bound, iflowface = cells + 1, None
_append_terms(out, _record_package(parent, text, used), cells, flow, bound, iflowface)
assert out["flowja"] is not None, "FLOW-JA-FACE not saved for this time step"
return out
def _aux_column(data: np.ndarray, name: str) -> np.ndarray | None:
"""Return a list record's auxiliary column (name matched case-insensitively), or None."""
found = [c for c in data.dtype.names if c.upper() == name]
return np.asarray(data[found[0]]) if found else None
def _append_terms(out: dict, package: Package, cells, flow, bound, iflowface=None) -> None:
"""Append one record's term entries (0-based cells, flows, bound numbers, IFLOWFACE) to out."""
out["packages"].append(package)
out["cell"].append(np.asarray(cells, dtype=np.int64))
out["flow"].append(np.asarray(flow, dtype=float))
out["bound"].append(np.asarray(bound, dtype=np.int64))
n = len(out["cell"][-1])
out["package"].append(np.full(n, len(out["packages"]) - 1, dtype=np.int64))
out["iflowface"].append(np.zeros(n) if iflowface is None else np.rint(iflowface))
def _concat(parts: list[np.ndarray], dtype) -> np.ndarray:
"""Concatenate a possibly empty list of arrays."""
return np.concatenate(parts).astype(dtype) if parts else np.zeros(0, dtype=dtype)
def _term_boundnames(
table: BoundnameTable, packages: list[Package], term_package, term_bound, kper: int
) -> np.ndarray:
"""Return the boundname of every term entry (looked up by the record's parent package)."""
out = np.full(len(term_bound), "", dtype=object)
for p, package in enumerate(packages):
sel = term_package == p
out[sel] = table.names(package.parent or package.name, kper, term_bound[sel])
return out
[docs]
@dataclass
class StepReader:
"""Grid and budget file opened once; reads a FlowStep for any saved time step.
The grid (IA/JA, IDOMAIN, areas) is read once and shared by every FlowStep it returns.
The budget and head files stay open until close() (or the end of a with block); open
handles stop Windows from deleting the model folder.
"""
grid: GridData
cbc: flopy.utils.CellBudgetFile
area: np.ndarray | None = None
boundnames: BoundnameTable | None = None
notes: dict | None = None
heads: flopy.utils.HeadFile | None = None
[docs]
def close(self) -> None:
"""Close the budget file and the head file."""
self.cbc.close()
if self.heads is not None:
self.heads.close()
def __enter__(self) -> "StepReader":
return self
def __exit__(self, *exc) -> None:
self.close()
@property
def steps(self) -> list[tuple[int, int]]:
"""0-based (kstp, kper) of every saved time step."""
return [(int(k), int(p)) for k, p in self.cbc.get_kstpkper()]
[docs]
def step_length(self, kstpkper: tuple[int, int]) -> float:
"""Return the length (DELT) of a saved time step, from the budget record header."""
return float(self.cbc.recordarray[_step_record_indices(self.cbc, kstpkper)[0]]["delt"])
[docs]
def thickness(
self, kstpkper: tuple[int, int], saturation: np.ndarray | None, active: np.ndarray
) -> tuple[np.ndarray | None, str]:
"""Return (saturated thickness, its source) from DATA-SAT, else heads, else (None, '')."""
grid = self.grid
if grid.top is None or (saturation is None and self.heads is None):
return None, ""
head = None
if saturation is None:
head = np.asarray(self.heads.get_data(kstpkper=tuple(kstpkper)), float).ravel()
thick = saturated_thickness(grid.top, grid.bot, grid.icelltype, active, saturation, head)
return thick, "DATA-SAT" if saturation is not None else "heads"
[docs]
def read(self, kstpkper: tuple[int, int] = (0, 0)) -> FlowStep:
"""Read one time step into a FlowStep."""
active = self.grid.idomain > 0
indices = _step_record_indices(self.cbc, kstpkper)
terms = _read_terms(self.cbc, indices, active)
assert terms["flowja"].size == self.grid.ja.size, "FLOW-JA-FACE length does not match NJA"
term_package = _concat(terms["package"], np.int64)
term_bound = _concat(terms["bound"], np.int64)
names = None
if self.boundnames is not None:
names = _term_boundnames(
self.boundnames, terms["packages"], term_package, term_bound, kstpkper[1]
)
header = self.cbc.recordarray[indices[0]]
thickness, thickness_source = self.thickness(kstpkper, terms["saturation"], active)
return FlowStep(
ia=self.grid.ia,
ja=self.grid.ja,
flowja=terms["flowja"],
term_cell=_concat(terms["cell"], np.int64),
term_flow=_concat(terms["flow"], float),
term_package=term_package,
term_bound=term_bound,
packages=terms["packages"],
term_boundname=names,
term_iflowface=_concat(terms["iflowface"], np.int64),
area=self.area if self.area is not None else self.grid.area,
saturated_thickness=thickness,
active=active,
grid_shape=self.grid.shape,
grid_type=self.grid.grid_type,
kstpkper=tuple(int(k) for k in kstpkper),
totim=float(header["totim"]),
notes={
**(self.notes or {}),
"delt": float(header["delt"]),
"thickness_from": thickness_source,
},
)
[docs]
def open_files(
grb_file: str | Path,
budget_file: str | Path,
area: np.ndarray | None = None,
head_file: str | Path | None = None,
) -> StepReader:
"""Open a binary grid file and a budget file (and optionally a head file) once.
The head file is only used for the saturated thickness when a step has no DATA-SAT.
"""
notes = {"grb_file": str(grb_file), "budget_file": str(budget_file)}
heads = None
head_file = Path(head_file) if head_file is not None else None
if head_file is not None and head_file.is_file() and head_file.stat().st_size > 0:
heads = flopy.utils.HeadFile(str(head_file))
notes["head_file"] = str(head_file)
return StepReader(
grid=read_grid(grb_file),
cbc=flopy.utils.CellBudgetFile(str(budget_file)),
area=None if area is None else np.asarray(area, dtype=float),
notes=notes,
heads=heads,
)
[docs]
def open_model(
ws: str | Path,
model: str | None = None,
area: np.ndarray | None = None,
boundnames: bool = False,
) -> StepReader:
"""Open a GWF model's grid and budget files once; optionally load its boundnames (FloPy).
For DISU grids without vertices the cell area is read from the DISU input (AREA).
"""
files = find_model_files(ws, model)
reader = open_files(files.grb_file, files.budget_file, area=area, head_file=files.head_file)
if reader.area is None and reader.grid.area is None and files.dis_type == "DISU":
reader.area = read_disu_area(ws, files.model)
if boundnames:
reader.boundnames = read_boundnames(files, ws)
reader.notes.update({"model": files.model, "units": read_units(files)})
return reader
[docs]
def read_step_files(
grb_file: str | Path,
budget_file: str | Path,
kstpkper: tuple[int, int] = (0, 0),
area: np.ndarray | None = None,
) -> FlowStep:
"""Read one time step from a binary grid file and a budget file into a FlowStep."""
with open_files(grb_file, budget_file, area=area) as reader:
return reader.read(kstpkper)
[docs]
def read_step(
ws: str | Path,
model: str | None = None,
kstpkper: tuple[int, int] = (0, 0),
area: np.ndarray | None = None,
boundnames: bool = False,
) -> FlowStep:
"""Read one time step of a GWF model in a simulation folder into a FlowStep.
For DISU grids without vertices the cell area is read from the DISU input (AREA).
boundnames=True also fills FlowStep.term_boundname from the package input (FloPy).
"""
with open_model(ws, model, area=area, boundnames=boundnames) as reader:
return reader.read(kstpkper)