"""
Module to read MODFLOW 6 binary grid files (*.grb) that define the model
grid binary output files. The module contains the MfGrdFile class that can
be accessed by the user.
"""
import warnings
import numpy as np
from ...utils.utils_def import FlopyBinaryData
warnings.simplefilter("always", DeprecationWarning)
def _crs_to_string(crs):
"""Convert a CRS value to a string suitable for writing to a GRB file.
Strings are returned unchanged (pass-through for DIS package / GRB
round-trip cases). Integers are treated as EPSG codes. Objects with
pyproj-compatible methods (``to_epsg``, ``to_authority``, ``to_wkt``) are
converted in that priority order so the shortest interoperable form is
preferred.
Parameters
----------
crs : str, int, or pyproj.CRS
CRS value to convert.
Returns
-------
str
CRS as a string ready for writing to the GRB CRS field.
"""
if isinstance(crs, str):
return crs
if isinstance(crs, int):
return f"EPSG:{crs}"
# Duck-type pyproj.CRS or any compatible object — no hard dependency.
if hasattr(crs, "to_epsg"):
epsg = crs.to_epsg()
if epsg is not None:
return f"EPSG:{epsg}"
if hasattr(crs, "to_authority"):
authority = crs.to_authority()
if authority is not None:
return ":".join(authority)
if hasattr(crs, "to_wkt"):
return crs.to_wkt()
return str(crs)
[docs]class MfGrdFile(FlopyBinaryData):
"""
The MfGrdFile class.
Parameters
----------
filename : str
Name of the MODFLOW 6 binary grid file
precision : string
'single' or 'double'. Default is 'double'.
verbose : bool
Write information to standard output. Default is False.
Attributes
----------
Methods
-------
See Also
--------
Notes
-----
The MfGrdFile class provides simple ways to retrieve data from binary
MODFLOW 6 binary grid files (.grb). The binary grid file contains data
that can be used for post processing MODFLOW 6 model results. For
example, the ia and ja arrays for a model grid.
Examples
--------
>>> import flopy
>>> gobj = flopy.utils.MfGrdFile('test.dis.grb')
"""
def __init__(self, filename, precision="double", verbose=False):
"""
Class constructor.
"""
# Call base class init
super().__init__()
# set attributes
self.precision = precision
self.verbose = verbose
self._initial_len = 50
self._recorddict = {}
self._datadict = {}
self._recordkeys = []
self.filename = filename
if self.verbose:
print(f"\nProcessing binary grid file: {filename}")
# open the grb file
self.file = open(filename, "rb")
# grid type
line = self.read_text(self._initial_len).strip()
t = line.split()
self._grid_type = t[1]
# version
line = self.read_text(self._initial_len).strip()
t = line.split()
self._version = int(t[1])
# ntxt
line = self.read_text(self._initial_len).strip()
t = line.split()
self._ntxt = int(t[1])
# length of text
line = self.read_text(self._initial_len).strip()
t = line.split()
self._lentxt = int(t[1])
# read text strings
for idx in range(self._ntxt):
line = self.read_text(self._lentxt).strip()
if line.startswith("#"):
continue
t = line.split()
key = t[0]
dt = t[1]
if dt == "INTEGER":
dtype = np.int32
elif dt == "SINGLE":
dtype = np.float32
elif dt == "DOUBLE":
dtype = np.float64
elif dt == "CHARACTER":
dtype = str
else:
dtype = None
nd = int(t[3])
if nd > 0:
shp = [int(v) for v in t[4:]]
shp = tuple(shp[::-1])
else:
shp = (0,)
self._recorddict[key] = (dtype, nd, shp)
self._recordkeys.append(key)
if self.verbose:
s = ""
if nd > 0:
s = shp
print(f" File contains data for {key} with shape {s}")
if self.verbose:
print(f"Attempting to read {len(self._recordkeys)} records from {filename}")
for key in self._recordkeys:
if self.verbose:
print(f" Reading {key}")
dt, nd, shp = self._recorddict[key]
# read array data
if nd > 0:
count = 1
for v in shp:
count *= v
if dt == str:
v = self.read_text(nchar=count)
else:
v = self.read_record(count=count, dtype=dt)
# read variable data
else:
if dt == np.int32:
v = self.read_integer()
elif dt == np.float32:
v = self._read_values(dt, 1)[0]
elif dt == np.float64:
v = self._read_values(dt, 1)[0]
self._datadict[key] = v
if self.verbose:
if nd == 0 or dt == str:
print(f" {key} = {v}")
else:
print(f" {key}: min = {v.min()} max = {v.max()}")
# close the file
self.file.close()
# initialize the model grid to None
self._modelgrid = None
# set ia and ja
self._set_iaja()
# internal functions
def _set_iaja(self):
"""
Set ia and ja from _datadict.
"""
self._ia = self._datadict["IA"] - 1
self._ja = self._datadict["JA"] - 1
def _set_modelgrid(self):
"""
Define structured, vertex, or unstructured grid based on MODFLOW 6
discretization type.
Returns
-------
modelgrid : grid
"""
from ...discretization.structuredgrid import StructuredGrid
from ...discretization.unstructuredgrid import UnstructuredGrid
from ...discretization.vertexgrid import VertexGrid
modelgrid = None
idomain = self.idomain
xorigin = self.xorigin
yorigin = self.yorigin
angrot = self.angrot
try:
top = self.top
botm = self.bot
if self._grid_type == "DISV":
nlay, ncpl = self.nlay, self.ncpl
vertices, cell2d = self.cell2d
top = np.ravel(top)
botm.shape = (nlay, ncpl)
modelgrid = VertexGrid(
vertices,
cell2d,
top,
botm,
idomain,
xoff=xorigin,
yoff=yorigin,
angrot=angrot,
)
elif self._grid_type == "DIS":
nlay, nrow, ncol = (
self.nlay,
self.nrow,
self.ncol,
)
delr, delc = self.delr, self.delc
top.shape = (nrow, ncol)
botm.shape = (nlay, nrow, ncol)
modelgrid = StructuredGrid(
delc,
delr,
top,
botm,
idomain=idomain,
xoff=xorigin,
yoff=yorigin,
angrot=angrot,
)
else:
iverts, verts = self.iverts, self.verts
vertc = self.cellcenters
xc, yc = vertc[:, 0], vertc[:, 1]
modelgrid = UnstructuredGrid(
vertices=verts,
iverts=iverts,
xcenters=xc,
ycenters=yc,
top=top,
botm=botm,
idomain=idomain,
xoff=xorigin,
yoff=yorigin,
angrot=angrot,
)
except:
print(f"could not set model grid for {self.file.name}")
self._modelgrid = modelgrid
return
def _build_vertices_cell2d(self):
"""
Build the mf6 vertices and cell2d array to generate a VertexGrid
Returns:
-------
vertices: list
cell2d: list
"""
iverts, verts = self.iverts, self.verts
vertc = self.cellcenters
vertices = [[ix] + list(i) for ix, i in enumerate(verts)]
cell2d = [
[ix] + list(vertc[ix]) + [len(i) - 1] + i[:-1]
for ix, i in enumerate(iverts)
]
return vertices, cell2d
def _get_iverts(self):
"""
Get a list of the vertices that define each model cell.
Returns
-------
iverts : list of lists
List with lists containing the vertex indices for each model cell.
"""
iverts = None
if "IAVERT" in self._datadict:
iverts = []
iavert = self.iavert
javert = self.javert
nsize = iavert.shape[0] - 1
for ivert in range(nsize):
i0 = iavert[ivert]
i1 = iavert[ivert + 1]
iverts.append((javert[i0:i1]).tolist())
if self.verbose:
print(f"returning iverts from {self.file.name}")
return iverts
def _get_verts(self):
"""
Get a list of the x, y pair for each vertex from the data in the
binary grid file.
Returns
-------
verts : np.ndarray
Array with x, y pairs for every vertex used to define the model.
"""
verts = None
if "VERTICES" in self._datadict:
shpvert = self._recorddict["VERTICES"][2]
verts = self._datadict["VERTICES"].reshape(shpvert)
if self._grid_type == "DISU":
# modify verts
verts = [
[idx, verts[idx, 0], verts[idx, 1]] for idx in range(shpvert[0])
]
if self.verbose:
print(f"returning verts from {self.file.name}")
return verts
def _get_cellcenters(self):
"""
Get the cell centers centroids for a MODFLOW 6 GWF model that uses
the DISV or DISU discretization.
Returns
-------
vertc : np.ndarray
Array with x, y pairs of the centroid for every model cell
"""
xycellcenters = None
if "CELLX" in self._datadict:
x = self._datadict["CELLX"]
y = self._datadict["CELLY"]
xycellcenters = np.column_stack((x, y))
if self.verbose:
print(f"returning cell centers from {self.file.name}")
return xycellcenters
# properties
@property
def version(self):
"""
MODFLOW 6 grid file version.
Returns
-------
version : int
"""
return self._version
@property
def grid_type(self):
"""
Grid type defined in the MODFLOW 6 grid file.
Returns
-------
grid_type : str
"""
return self._grid_type
@property
def nlay(self):
"""
Number of layers. None for DISU grids.
Returns
-------
nlay : int
"""
nlay = None
if "NLAY" in self._datadict:
nlay = self._datadict["NLAY"]
return nlay
@property
def nrow(self):
"""
Number of rows. None for DISV and DISU grids.
Returns
-------
nrow : int
"""
nrow = None
if "NROW" in self._datadict:
nrow = self._datadict["NROW"]
return nrow
@property
def ncol(self):
"""
Number of columns. None for DISV and DISU grids.
Returns
-------
ncol : int
"""
ncol = None
if "NCOL" in self._datadict:
ncol = self._datadict["NCOL"]
return ncol
@property
def ncpl(self):
"""
Number of cells per layer. None for DISU grids.
Returns
-------
ncpl : int
"""
ncpl = None
if "NCPL" in self._datadict:
ncpl = self._datadict["NCPL"]
return ncpl
@property
def ncells(self):
"""
Number of cells.
Returns
-------
ncells : int
"""
# disu is the only grid that has the number of cells
# set to nodes. All other grids use NCELLS in grb
if "NCELLS" in self._datadict:
ncells = self._datadict["NCELLS"]
elif "NODES" in self._datadict:
ncells = self._datadict["NODES"]
else:
ncells = None
return ncells
@property
def nodes(self):
"""
Number of nodes.
Returns
-------
nodes : int
"""
nodes = self.ncells
return nodes
@property
def shape(self):
"""
Shape of the model grid (tuple).
Returns
-------
shape : tuple
"""
if self._grid_type == "DIS":
shape = (self.nlay, self.nrow, self.ncol)
elif self._grid_type == "DIS2D":
shape = (self.nrow, self.ncol)
elif self._grid_type == "DISV":
shape = (self.nlay, self.ncpl)
elif self._grid_type == "DISV2D":
shape = (self.ncells,)
elif self._grid_type == "DISV1D":
shape = (self.ncells,)
elif self._grid_type == "DISU":
shape = (self.nodes,)
else:
shape = None
return shape
@property
def xorigin(self):
"""
x-origin of the model grid. None if not defined in the
MODFLOW 6 grid file.
Returns
-------
xorigin : float
"""
if "XORIGIN" in self._datadict:
xorigin = self._datadict["XORIGIN"]
else:
xorigin = None
return xorigin
@property
def yorigin(self):
"""
y-origin of the model grid. None if not defined in the
MODFLOW 6 grid file.
Returns
-------
yorigin : float
"""
if "YORIGIN" in self._datadict:
yorigin = self._datadict["YORIGIN"]
else:
yorigin = None
return yorigin
@property
def angrot(self):
"""
Model grid rotation angle. None if not defined in the
MODFLOW 6 grid file.
Returns
-------
angrot : float
"""
if "ANGROT" in self._datadict:
angrot = self._datadict["ANGROT"]
else:
angrot = None
return angrot
@property
def idomain(self):
"""
IDOMAIN for the model grid. None if not defined in the
MODFLOW 6 grid file.
Returns
-------
idomain : ndarray of ints
"""
if "IDOMAIN" in self._datadict:
idomain = self._datadict["IDOMAIN"]
else:
idomain = None
return idomain
@property
def delr(self):
"""
Cell size in the row direction (y-direction). None if not
defined in the MODFLOW 6 grid file.
Returns
-------
delr : ndarray of floats
"""
delr = None
if "DELR" in self._datadict:
delr = self._datadict["DELR"]
return delr
@property
def delc(self):
"""
Cell size in the column direction (x-direction). None if not
defined in the MODFLOW 6 grid file.
Returns
-------
delc : ndarray of floats
"""
delc = None
if "DELC" in self._datadict:
delc = self._datadict["DELC"]
return delc
@property
def top(self):
"""
Top of the model cells in the upper model layer for DIS and
DISV grids. Top of the model cells for DISU grids.
Returns
-------
top : ndarray of floats
"""
top = None
if "TOP" in self._datadict:
top = self._datadict["TOP"]
return top
@property
def bot(self):
"""
Bottom of the model cells.
Returns
-------
bot : ndarray of floats
"""
bot = None
if "BOTM" in self._datadict:
bot = self._datadict["BOTM"]
elif "BOT" in self._datadict:
bot = self._datadict["BOT"]
return bot
@property
def nja(self):
"""
Number of non-zero entries JA vector array.
Returns
-------
nja : int
"""
return self._datadict["NJA"]
@property
def ia(self):
"""
index array that defines indexes for `.ja`. Each ia value is the
starting position of data for a cell. [ia[n]:ia[n+1]] would give you
all data for a cell. ia[n] is also the location of data for the
diagonal position. See `.ja` property documentation
for an example of getting a cell's number and connected cells
Returns
-------
ia : ndarray of ints
"""
return np.array(self._ia, dtype=int)
@property
def ja(self):
"""
Flat jagged connection array for a model. `.ja` for a cell includes the
cell number and the cell number for all connected cells. Indexes for
cells are stored in the `.ia` variable.
Returns
-------
ja : ndarray of ints
Examples
--------
>>> from flopy.mf6.utils import MfGrdFile
>>> grb = MfGrdFile("my_model.dis.grb")
>>> ia = grb.ia
>>> ja = grb.ja
>>> # get connections for node 0
>>> ja_node0 = ja[ia[0]:ia[1]]
>>> node = ja_node0[0]
>>> connections = ja_node0[1:]
"""
return self._ja
@property
def iavert(self):
"""
index array that defines indexes for `.javart`. Each ia value is the
starting position of data for a cell. [iavert[n]:iavert[n+1]] would
give you all data for a cell. See `.javert` property documentation for
an example of getting cell number and it's vertex numbers.
Alternatively, the `.iverts` property can be used to get this
information
Returns
-------
iavert : ndarray of ints or None for structured grids
"""
if "IAVERT" in self._datadict:
iavert = self._datadict["IAVERT"] - 1
else:
iavert = None
return iavert
@property
def javert(self):
"""
Flat jagged array of vertex numbers that comprise all of the cells
Returns
-------
javerts : ndarray of ints or None for structured grids
Examples
--------
>>> from flopy.mf6.utils import MfGrdFile
>>> grb = MfGrdFile("my_model.dis.grb")
>>> iavert = self.iavert
>>> javert = self.javert
>>> # get vertex numbers for node 0
>>> vertnums = javert[iavert[0]:iavert[1]]
"""
if "JAVERT" in self._datadict:
javert = self._datadict["JAVERT"] - 1
else:
javert = None
return javert
@property
def iverts(self):
"""
Vertex numbers comprising each cell for every cell in model grid.
Returns
-------
iverts : list of lists of ints
"""
return self._get_iverts()
@property
def verts(self):
"""
x,y location of each vertex that defines the model grid.
Returns
-------
verts : ndarray of floats
"""
return self._get_verts()
@property
def cellcenters(self):
"""
Cell centers (x,y).
Returns
-------
cellcenters : ndarray of floats
"""
return self._get_cellcenters()
@property
def modelgrid(self):
"""
Model grid object.
Returns
-------
modelgrid : StructuredGrid, VertexGrid, UnstructuredGrid
"""
if self._modelgrid is None:
self._set_modelgrid()
return self._modelgrid
@property
def cell2d(self):
"""
cell2d data for a DISV grid. None for DIS and DISU grids.
Returns
-------
cell2d : list of lists
"""
if self._grid_type in ("DISV", "DISV2D", "DISV1D"):
vertices, cell2d = self._build_vertices_cell2d()
else:
vertices, cell2d = None, None
return vertices, cell2d
[docs] def export(
self,
filename,
*,
precision=None,
version=None,
crs=None,
verbose=False,
):
"""
Export the binary grid file to a new file.
Parameters
----------
filename : str or PathLike
Path to output .grb file
precision : str, optional
'single' or 'double'. If None, uses the precision from the
original file (default None)
version : int, optional
Grid file version to write. If None (default), version 2 is
written when a CRS string is available, otherwise version 1.
Pass ``version=1`` to force version 1 even when a CRS is set.
crs : str, optional
CRS user input string to write (e.g. ``"EPSG:26916"`` or OGC
WKT). If None, uses the CRS from the source file if present.
Providing this argument with a version 1 source file upgrades
the output to version 2.
verbose : bool, optional
Print progress messages (default False)
All parameters except ``filename`` are keyword-only, so adding a
new parameter here in the future (e.g. a dedicated parameter and
property for a new GRB field, mirroring ``crs``/``.crs``) can
never silently change the meaning of an existing positional call.
Examples
--------
>>> from flopy.mf6.utils import MfGrdFile
>>> grb = MfGrdFile('model.dis.grb')
>>> grb.export('model_copy.dis.grb')
>>> # Round-trip a v2 file preserving CRS
>>> grb2 = MfGrdFile('model_v2.dis.grb')
>>> grb2.export('model_v2_copy.dis.grb')
>>> # Upgrade a v1 file to v2 with an explicit CRS
>>> grb.export('model_v2.dis.grb', crs='EPSG:26916')
"""
if precision is None:
precision = self.precision
if isinstance(crs, str) and crs.strip() == "":
crs = None
raw_crs = crs if crs is not None else self.crs
effective_crs = _crs_to_string(raw_crs) if raw_crs is not None else None
if effective_crs is not None and effective_crs.strip() == "":
# A blank CRS carries no usable information. Treat it the
# same as no CRS rather than writing an empty crs string.
effective_crs = None
if version is None:
version = 2 if effective_crs is not None else 1
float_type = "SINGLE" if precision.lower() == "single" else "DOUBLE"
# Build var_list and data_dict from the parsed record structure so that
# array dimensions come directly from the source file rather than being
# recomputed. This avoids shape mismatches and naturally handles any
# grid type.
var_list = []
data_dict = {}
for key in self._recordkeys:
if key == "CRS":
continue # handled by the version check below
dt, nd, shp = self._recorddict[key]
if dt == np.int32:
dtype_str = "INTEGER"
elif dt in (np.float32, np.float64):
dtype_str = float_type
elif dt == str:
dtype_str = "CHARACTER"
else:
dtype_str = float_type
# shp is stored in reversed (numpy) order; the write loop reverses
# it back to Fortran (definition-line) order via dims[::-1].
dims = list(shp) if nd > 0 else []
var_list.append((key, dtype_str, nd, dims))
data_dict[key] = self._datadict[key]
if version >= 2:
if effective_crs is None:
raise ValueError(
"version=2 requires a CRS string. Provide crs= or use a "
"version 2 source file."
)
_MF6_CRS_MAXLEN = 5000
if len(effective_crs) > _MF6_CRS_MAXLEN:
raise ValueError(
f"CRS string length {len(effective_crs)} exceeds the "
f"MODFLOW 6 maximum of {_MF6_CRS_MAXLEN} characters."
)
var_list.append(("CRS", "CHARACTER", 1, [len(effective_crs)]))
data_dict["CRS"] = effective_crs
ntxt = len(var_list)
lentxt = 100
if verbose:
print(f"Writing binary grid file: {filename}")
print(f" Grid type: {self.grid_type}")
print(f" Version: {version}")
print(f" Number of variables: {ntxt}")
# Create writer with appropriate precision
writer = FlopyBinaryData()
writer.precision = precision
with open(filename, "wb") as f:
writer.file = f
# Write text header lines and definition lines in MODFLOW 6 format:
# content left-aligned, space-padded to (nchar-1), then newline at
# position (nchar-1). This matches what MF6 writes and is required
# for Fortran readers that parse these as fixed-length text records.
def write_line(text, nchar):
writer.file.write(text.encode("ascii").ljust(nchar - 1) + b"\n")
header_len = 50
write_line(f"GRID {self.grid_type}", header_len)
write_line(f"VERSION {version}", header_len)
write_line(f"NTXT {ntxt}", header_len)
write_line(f"LENTXT {lentxt}", header_len)
# Write variable definition lines (100 chars each)
for name, dtype_str, ndim, dims in var_list:
if ndim == 0:
line = f"{name} {dtype_str} NDIM {ndim}"
else:
dims_str = " ".join(
str(d) for d in dims[::-1]
) # Reverse for Fortran order
line = f"{name} {dtype_str} NDIM {ndim} {dims_str}"
write_line(line, lentxt)
# Write binary data for each variable
for name, dtype_str, ndim, dims in var_list:
if name not in data_dict:
raise ValueError(
f"Required variable '{name}' not found in grid file"
)
value = data_dict[name]
if verbose:
if ndim == 0:
print(f" Writing {name} = {value}")
else:
if hasattr(value, "min"):
print(
f" Writing {name}: min = {value.min()} max = {value.max()}"
)
else:
print(f" Writing {name}")
# Write scalar or array data
if dtype_str == "CHARACTER":
writer.write_text(value, len(value))
elif ndim == 0:
# Scalar value
if dtype_str == "INTEGER":
writer.write_integer(int(value))
elif dtype_str in ("DOUBLE", "SINGLE"):
writer.write_real(float(value))
else:
# Array data
arr = np.asarray(value)
if dtype_str == "INTEGER":
arr = arr.astype(np.int32)
elif dtype_str == "DOUBLE":
arr = arr.astype(np.float64)
elif dtype_str == "SINGLE":
arr = arr.astype(np.float32)
writer.write_record(np.ravel(arr, order="C"), dtype=arr.dtype)
if verbose:
print(f"Successfully wrote {filename}")
@property
def crs(self):
"""
CRS user input string (version 2 GRB file only).
Returns
-------
crs : str or None
"""
crs = None
if "CRS" in self._datadict:
crs = self._datadict["CRS"]
return crs