import copy
import os
import numpy as np
from matplotlib.path import Path
from ..utils.geometry import is_clockwise, transform
from .grid import CachedData, Grid
[docs]class VertexGrid(Grid):
"""
class for a vertex model grid
Parameters
----------
vertices
list of vertices that make up the grid
cell2d
list of cells and their vertices
top : list or ndarray
top elevations for all cells in the grid.
botm : list or ndarray
bottom elevations for all cells in the grid.
idomain : int or ndarray
ibound/idomain value for each cell
lenuni : int or ndarray
model length units
crs : pyproj.CRS, int, str, optional if `prjfile` is specified
Coordinate reference system (CRS) for the model grid
(must be projected; geographic CRS are not supported).
The value can be anything accepted by
:meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
such as an authority string (eg "EPSG:26916") or a WKT string.
prjfile : str or PathLike, optional if `crs` is specified
ESRI-style projection file with well-known text defining the CRS
for the model grid (must be projected; geographic CRS are not supported).
xoff : float
x coordinate of the origin point (lower left corner of model grid)
in the spatial reference coordinate system
yoff : float
y coordinate of the origin point (lower left corner of model grid)
in the spatial reference coordinate system
angrot : float
rotation angle of model grid, as it is rotated around the origin point
**kwargs : dict, optional
Support deprecated keyword options.
.. deprecated:: 3.5
The following keyword options will be removed for FloPy 3.6:
- ``prj`` (str or PathLike): use ``prjfile`` instead.
- ``epsg`` (int): use ``crs`` instead.
- ``proj4`` (str): use ``crs`` instead.
Properties
----------
vertices
returns list of vertices that make up the grid
cell2d
returns list of cells and their vertices
Methods
-------
get_cell_vertices(cellid)
returns vertices for a single cell at cellid.
"""
def __init__(
self,
vertices=None,
cell2d=None,
top=None,
botm=None,
idomain=None,
lenuni=None,
crs=None,
prjfile=None,
xoff=0.0,
yoff=0.0,
angrot=0.0,
nlay=None,
ncpl=None,
cell1d=None,
**kwargs,
):
super().__init__(
"vertex",
top=top,
botm=botm,
idomain=idomain,
lenuni=lenuni,
crs=crs,
prjfile=prjfile,
xoff=xoff,
yoff=yoff,
angrot=angrot,
**kwargs,
)
self._vertices = vertices
self._cell1d = cell1d
self._cell2d = cell2d
if botm is None:
self._nlay = nlay
self._ncpl = ncpl
else:
self._nlay = None
self._ncpl = None
@property
def is_valid(self):
if self._vertices is not None and (
self._cell2d is not None or self._cell1d is not None
):
return True
return False
@property
def is_complete(self):
if (
self._vertices is not None
and (self._cell2d is not None or self._cell1d is not None)
and super().is_complete
):
return True
return False
@property
def nlay(self):
if self._cell1d is not None:
return 1
elif self._botm is not None:
return len(self._botm)
else:
return self._nlay
@property
def ncpl(self):
if self._cell1d is not None:
return len(self._cell1d)
if self._botm is not None:
if self._botm.ndim == 2: # (nlay, ncpl)
return self._botm.shape[1]
elif self._botm.ndim == 1: # (ncpl,)
return self._botm.shape[0]
if self._cell2d is not None and self._nlay is None:
return len(self._cell2d)
else:
return self._ncpl
@property
def nnodes(self):
return self.nlay * self.ncpl
@property
def nvert(self):
return len(self._vertices)
@property
def iverts(self):
if self._cell2d is not None:
return [list(t)[4:] for t in self.cell2d]
elif self._cell1d is not None:
return [list(t)[3:] for t in self.cell1d]
@property
def cell1d(self):
if self._cell1d is not None:
return [[ivt for ivt in t if ivt is not None] for t in self._cell1d]
@property
def cell2d(self):
if self._cell2d is not None:
return [[ivt for ivt in t if ivt is not None] for t in self._cell2d]
@property
def verts(self):
verts = np.array([list(t)[1:] for t in self._vertices], dtype=float).T
x, y = transform(
verts[0], verts[1], self.xoffset, self.yoffset, self.angrot_radians
)
return np.array(list(zip(x, y)))
@property
def shape(self):
return self.nlay, self.ncpl
@property
def top_botm(self):
new_top = np.expand_dims(self._top, 0)
return np.concatenate((new_top, self._botm), axis=0)
@property
def extent(self):
self._copy_cache = False
xvertices = np.hstack(self.xvertices)
yvertices = np.hstack(self.yvertices)
self._copy_cache = True
return (
np.min(xvertices),
np.max(xvertices),
np.min(yvertices),
np.max(yvertices),
)
@property
def grid_lines(self):
"""
Creates a series of grid line vertices for drawing
a model grid line collection
Returns:
list: grid line vertices
"""
self._copy_cache = False
xgrid = self.xvertices
ygrid = self.yvertices
# close the cell by connecting the last vertex with the first
close_cell = True
if self.cell1d is not None:
close_cell = False
# go through each cell and create a line segment for each face
lines = []
ncpl = len(xgrid)
for icpl in range(ncpl):
xcoords = xgrid[icpl]
ycoords = ygrid[icpl]
npoints = len(xcoords)
for ipoint in range(npoints - 1):
lines.append(
[
(xcoords[ipoint], ycoords[ipoint]),
(xcoords[ipoint + 1], ycoords[ipoint + 1]),
]
)
if close_cell:
lines.append([(xcoords[-1], ycoords[-1]), (xcoords[0], ycoords[0])])
self._copy_cache = True
return lines
@property
def xyzcellcenters(self):
"""
Method to get cell centers and set to grid
"""
cache_index = "cellcenters"
if (
cache_index not in self._cache_dict
or self._cache_dict[cache_index].out_of_date
):
self._build_grid_geometry_info()
if self._copy_cache:
return self._cache_dict[cache_index].data
else:
return self._cache_dict[cache_index].data_nocopy
@property
def xyzvertices(self):
"""
Method to get all grid vertices in a layer, arranged per cell
Returns:
list of size sum(nvertices per cell)
"""
cache_index = "xyzgrid"
if (
cache_index not in self._cache_dict
or self._cache_dict[cache_index].out_of_date
):
self._build_grid_geometry_info()
if self._copy_cache:
return self._cache_dict[cache_index].data
else:
return self._cache_dict[cache_index].data_nocopy
@property
def map_polygons(self):
"""
Get a list of matplotlib Polygon patches for plotting
Returns
-------
list of Polygon objects
"""
cache_index = "xyzgrid"
if (
cache_index not in self._cache_dict
or self._cache_dict[cache_index].out_of_date
):
self.xyzvertices
self._polygons = None
if self._polygons is None:
self._polygons = [
Path(self.get_cell_vertices(nn)) for nn in range(self.ncpl)
]
return copy.copy(self._polygons)
[docs] def to_geodataframe(self):
"""
Returns a geopandas GeoDataFrame of the model grid
Returns
-------
GeoDataFrame
"""
cache_index = "gdf_polys"
if (
cache_index not in self._cache_dict
or self._cache_dict[cache_index].out_of_date
):
polys = [[self.get_cell_vertices(nn)] for nn in range(self.ncpl)]
self._cache_dict[cache_index] = CachedData(polys)
else:
polys = self._cache_dict[cache_index].data_nocopy
featuretype = "Polygon"
if self._cell1d is not None:
featuretype = "multilinestring"
gdf = super().to_geodataframe(polys, featuretype)
if self.idomain is not None:
active = np.sum(
self.idomain.reshape(
(self.nlay, self.ncpl),
),
axis=0,
)
active = np.where(active > 0, 1, 0)
gdf["active"] = active
else:
gdf["active"] = 1
return gdf
[docs] def grid_line_geodataframe(self):
"""
Method to get a GeoDataFrame of grid lines
Returns
-------
GeoDataFrame
"""
gdf = super().to_geodataframe(self.grid_lines, featuretype="LineString")
gdf = gdf.rename(columns={"node": "number"})
return gdf
@property
def geo_dataframe(self):
"""
DEPRECATED -- Use to_geodataframe() instead. Will be removed in 3.11
Returns a geopandas GeoDataFrame of the model grid
Returns
-------
GeoDataFrame
"""
import warnings
warnings.warn(
"geo_dataframe has been deprecated, use to_geodataframe() instead",
DeprecationWarning,
)
return self.to_geodataframe()
[docs] def convert_grid(self, factor):
"""
Method to scale the model grid based on user supplied scale factors
Parameters
----------
factor
Returns
-------
Grid object
"""
if self.is_complete:
return VertexGrid(
vertices=[[i[0], i[1] * factor, i[2] * factor] for i in self._vertices],
cell2d=[
[i[0], i[1] * factor, i[2] * factor] + i[3:] for i in self._cell2d
],
top=self.top * factor,
botm=self.botm * factor,
idomain=self.idomain,
xoff=self.xoffset * factor,
yoff=self.yoffset * factor,
angrot=self.angrot,
)
else:
raise AssertionError("Grid is not complete and cannot be converted")
[docs] def intersect(self, x, y, z=None, local=False, forgive=False):
"""
Get the CELL2D number of a point with coordinates x and y
When the point is on the edge of two cells, the cell with the lowest
CELL2D number is returned.
Supports both scalar and array inputs for vectorized operations.
Parameters
----------
x : float or array-like
The x-coordinate(s) of the requested point(s)
y : float or array-like
The y-coordinate(s) of the requested point(s)
z : float, array-like, or None
optional, z-coordinate(s) of the requested point(s) will return
(lay, icell2d)
local: bool (optional)
If True, x and y are in local coordinates (defaults to False)
forgive: bool (optional)
Forgive x,y arguments that fall outside the model grid and
return NaNs instead (defaults to False - will throw exception)
Returns
-------
icell2d : int or ndarray
The CELL2D number(s). Returns int for scalar input,
ndarray for array input.
lay : int or ndarray (only if z is provided)
The layer number(s). Returns int for scalar input,
ndarray for array input.
"""
# Check if inputs are scalar
x_is_scalar = np.isscalar(x)
y_is_scalar = np.isscalar(y)
z_is_scalar = z is None or np.isscalar(z)
is_scalar_input = x_is_scalar and y_is_scalar and z_is_scalar
# Convert to arrays for uniform processing
x = np.atleast_1d(x)
y = np.atleast_1d(y)
if z is not None:
z = np.atleast_1d(z)
# Validate array shapes
if len(x) != len(y):
raise ValueError("x and y must have the same length")
if z is not None and len(z) != len(x):
raise ValueError("z must have the same length as x and y")
if local:
# transform x and y to real-world coordinates
x, y = super().get_coords(x, y)
xv, yv, zv = self.xyzvertices
# Initialize result arrays
n_points = len(x)
results = np.full(n_points, np.nan if forgive else -1, dtype=float)
if z is not None:
lays = np.full(n_points, np.nan if forgive else -1, dtype=float)
# Process each point
for i in range(n_points):
xi, yi = x[i], y[i]
found = False
# Check each cell
for icell2d in range(self.ncpl):
xa = np.array(xv[icell2d])
ya = np.array(yv[icell2d])
# x and y at least have to be within the bounding box of the cell
if (
np.any(xi <= xa)
and np.any(xi >= xa)
and np.any(yi <= ya)
and np.any(yi >= ya)
):
path = Path(np.stack((xa, ya)).transpose())
# use a small radius, so that the edge of the cell is included
if is_clockwise(xa, ya):
radius = -1e-9
else:
radius = 1e-9
if path.contains_point((xi, yi), radius=radius):
results[i] = icell2d
found = True
if z is not None:
zi = z[i]
for lay in range(self.nlay):
if (
self.top_botm[lay, icell2d]
>= zi
>= self.top_botm[lay + 1, icell2d]
):
lays[i] = lay
break
break
if not found and not forgive:
raise ValueError(
f"point given is outside of the model area: ({xi}, {yi})"
)
# Return results
if z is None:
if is_scalar_input:
result = results[0]
return int(result) if not np.isnan(result) else np.nan
else:
valid_mask = ~np.isnan(results)
return results.astype(int) if np.all(valid_mask) else results
else:
if is_scalar_input:
lay, icell2d = lays[0], results[0]
if not np.isnan(lay) and not np.isnan(icell2d):
return int(lay), int(icell2d)
else:
return np.nan, np.nan
else:
valid_mask = ~np.isnan(lays) & ~np.isnan(results)
return (
lays.astype(int) if np.all(valid_mask) else lays,
results.astype(int) if np.all(valid_mask) else results,
)
[docs] def get_cell_vertices(self, cellid=None, node=None):
"""
Get a set of cell vertices for a single cell.
Parameters
----------
cellid : int or tuple, optional
Cell identifier. Can be:
- cell2d index (int, 0 to ncpl-1)
- node number (int, >= ncpl) - will be converted to cell2d
- (cell2d,) single-element tuple
- (layer, cell2d) tuple (layer is ignored, vertices are 2D)
node : int, optional
Node number, mutually exclusive with cellid
Returns
-------
list
list of (x, y) cell vertex coordinates
Examples
--------
>>> import flopy
>>> from flopy.utils.gridutil import get_disv_kwargs
>>> disv_props = get_disv_kwargs(1, 10, 10, 1.0, 1.0, 1.0, [0.0])
>>> vg = flopy.discretization.VertexGrid(**disv_props)
>>> vg.get_cell_vertices(5) # cell2d index
>>> vg.get_cell_vertices((0, 5)) # (layer, cell2d) tuple
>>> vg.get_cell_vertices(node=105) # node number
>>> vg.get_cell_vertices(cellid=(1, 5)) # explicit cellid kwarg
"""
# Handle arguments
if cellid is not None and node is not None:
raise ValueError("cellid and node are mutually exclusive")
if cellid is None and node is None:
raise TypeError("expected cellid or node argument")
# Use cellid if provided, otherwise use node
if node is not None:
idx = node
else:
idx = cellid
# Handle tuple forms
if isinstance(idx, (tuple, list)):
if len(idx) == 1:
# (cell2d,) or (node,)
idx = idx[0]
elif len(idx) == 2:
# (layer, cell2d) - ignore layer since vertices are 2D
_, idx = idx
else:
raise ValueError(
f"cellid tuple must have 1 or 2 elements, got {len(idx)}"
)
# Convert node to cell2d if necessary
while idx >= self.ncpl:
if idx > self.nnodes:
raise IndexError(
f"node number {idx} exceeds grid node count {self.nnodes}"
)
idx -= self.ncpl
self._copy_cache = False
cell_verts = list(zip(self.xvertices[idx], self.yvertices[idx]))
self._copy_cache = True
return cell_verts
[docs] def get_node(self, cellids, node2d=False):
"""
Get node number from a list of zero-based MODFLOW
(layer, cell2d) tuples.
Parameters
----------
cellid_list : tuple of int or list of tuple of int
Zero-based (layer, cell2d) tuples
node2d : bool, optional
If True, return 2D node numbers (cell2d values).
If False (default), return 3D node numbers.
Returns
-------
list
list of MODFLOW nodes for each (layer, cell2d) tuple
in the input list
Examples
--------
>>> import flopy
>>> vg = flopy.discretization.VertexGrid(nlay=3, ncpl=100, ...)
>>> vg.get_node((0, 5))
[5]
>>> vg.get_node((1, 5))
[105]
>>> vg.get_node([(0, 5), (1, 5)], node2d=True)
[5, 5]
"""
if not isinstance(cellids, list):
cellids = [cellids]
# Validate
for cellid in cellids:
if len(cellid) != 2:
raise ValueError("VertexGrid cellid must be (layer, cell2d) tuple")
if node2d:
return [cell2d for lay, cell2d in cellids]
else:
nodes = []
for lay, cell2d in cellids:
if lay < 0 or lay >= self.nlay:
raise IndexError(f"Layer {lay} out of range [0, {self.nlay})")
if cell2d < 0 or cell2d >= self.ncpl:
raise IndexError(f"Cell2d {cell2d} out of range [0, {self.ncpl})")
nodes.append(lay * self.ncpl + cell2d)
return nodes
[docs] def plot(self, **kwargs):
"""
Plot the grid lines.
Parameters
----------
kwargs : ax, colors. The remaining kwargs are passed into the
the LineCollection constructor.
Returns
-------
lc : matplotlib.collections.LineCollection
"""
from ..plot import PlotMapView
mm = PlotMapView(modelgrid=self)
return mm.plot_grid(**kwargs)
def _build_grid_geometry_info(self):
cache_index_cc = "cellcenters"
cache_index_vert = "xyzgrid"
xcenters = []
ycenters = []
xvertices = []
yvertices = []
if self._cell1d is not None:
zcenters = []
zvertices = []
vertexdict = {v[0]: [v[1], v[2]] for v in self._vertices}
for cell1d in self.cell1d:
cell1d = tuple(cell1d)
xcenters.append(cell1d[1])
ycenters.append(cell1d[2])
zcenters.append(0.0)
vert_number = []
for i in cell1d[3:]:
vert_number.append(int(i))
xcellvert = []
ycellvert = []
zcellvert = []
for ix in vert_number:
xcellvert.append(vertexdict[ix][0])
ycellvert.append(vertexdict[ix][1])
zcellvert.append(0.0)
xvertices.append(xcellvert)
yvertices.append(ycellvert)
zvertices.append(zcellvert)
else:
vertexdict = {v[0]: [v[1], v[2]] for v in self._vertices}
# build xy vertex and cell center info
for cell2d in self.cell2d:
cell2d = tuple(cell2d)
xcenters.append(cell2d[1])
ycenters.append(cell2d[2])
vert_number = []
for i in cell2d[4:]:
vert_number.append(int(i))
xcellvert = []
ycellvert = []
for ix in vert_number:
xcellvert.append(vertexdict[ix][0])
ycellvert.append(vertexdict[ix][1])
xvertices.append(xcellvert)
yvertices.append(ycellvert)
# build z cell centers
zvertices, zcenters = self._zcoords()
if self._has_ref_coordinates:
# transform x and y
xcenters, ycenters = self.get_coords(xcenters, ycenters)
xvertxform = []
yvertxform = []
# vertices are a list within a list
for xcellvertices, ycellvertices in zip(xvertices, yvertices):
xcellvertices, ycellvertices = self.get_coords(
xcellvertices, ycellvertices
)
xvertxform.append(xcellvertices)
yvertxform.append(ycellvertices)
xvertices = xvertxform
yvertices = yvertxform
self._cache_dict[cache_index_cc] = CachedData(
[np.array(xcenters), np.array(ycenters), np.array(zcenters)]
)
self._cache_dict[cache_index_vert] = CachedData(
[xvertices, yvertices, zvertices]
)
[docs] def get_xvertices_for_layer(self, layer):
xgrid = np.array(self.xvertices, dtype=object)
return xgrid
[docs] def get_yvertices_for_layer(self, layer):
ygrid = np.array(self.yvertices, dtype=object)
return ygrid
[docs] def get_xcellcenters_for_layer(self, layer):
xcenters = np.array(self.xcellcenters)
return xcenters
[docs] def get_ycellcenters_for_layer(self, layer):
ycenters = np.array(self.ycellcenters)
return ycenters
[docs] def get_number_plottable_layers(self, a):
"""
Calculate and return the number of 2d plottable arrays that can be
obtained from the array passed (a)
Parameters
----------
a : ndarray
array to check for plottable layers
Returns
-------
nplottable : int
number of plottable layers
"""
nplottable = 0
required_shape = self.get_plottable_layer_shape()
if a.shape == required_shape:
nplottable = 1
else:
nplottable = a.size / self.ncpl
nplottable = int(nplottable)
return nplottable
[docs] def get_plottable_layer_array(self, a, layer):
# ensure plotarray is 1d with length ncpl
required_shape = self.get_plottable_layer_shape()
if a.ndim == 3:
if a.shape[0] == 1:
a = np.squeeze(a, axis=0)
plotarray = a[layer, :]
elif a.shape[1] == 1:
a = np.squeeze(a, axis=1)
plotarray = a[layer, :]
else:
raise ValueError(
"Array has 3 dimensions so one of them must be of size 1 "
"for a VertexGrid."
)
elif a.ndim == 2:
plotarray = a[layer, :]
elif a.ndim == 1:
plotarray = a
if plotarray.shape[0] == self.nnodes:
plotarray = plotarray.reshape(self.nlay, self.ncpl)
plotarray = plotarray[layer, :]
else:
raise ValueError("Array to plot must be of dimension 1 or 2")
msg = f"{plotarray.shape[0]} /= {required_shape}"
assert plotarray.shape == required_shape, msg
return plotarray
# initialize grid from a grb file
[docs] @classmethod
def from_binary_grid_file(cls, file_path, verbose=False):
"""
Instantiate a VertexGrid model grid from a MODFLOW 6 binary
grid (*.grb) file.
Parameters
----------
file_path : str
file path for the MODFLOW 6 binary grid file
verbose : bool
Write information to standard output. Default is False.
Returns
-------
return : VertexGrid
"""
from ..mf6.utils.binarygrid_util import MfGrdFile
grb_obj = MfGrdFile(file_path, verbose=verbose)
if grb_obj.grid_type != "DISV":
raise ValueError(
f"Binary grid file ({os.path.basename(file_path)}) "
"is not a vertex (DISV) grid."
)
idomain = grb_obj.idomain
xorigin = grb_obj.xorigin
yorigin = grb_obj.yorigin
angrot = grb_obj.angrot
nlay, ncpl = grb_obj.nlay, grb_obj.ncpl
top = np.ravel(grb_obj.top)
botm = grb_obj.bot
botm.shape = (nlay, ncpl)
vertices, cell2d = grb_obj.cell2d
return cls(
vertices,
cell2d,
top,
botm,
idomain,
xoff=xorigin,
yoff=yorigin,
angrot=angrot,
)