import copy
import os.path
from os import PathLike
from typing import Union
import numpy as np
from .grid import CachedData, Grid
[docs]def array_at_verts_basic2d(a):
"""
Computes values at cell vertices on 2d array using neighbor averaging.
Parameters
----------
a : ndarray
Array values at cell centers, could be a slice in any orientation.
Returns
-------
averts : ndarray
Array values at cell vertices, shape (a.shape[0]+1, a.shape[1]+1).
"""
assert a.ndim == 2
shape_verts2d = (a.shape[0] + 1, a.shape[1] + 1)
# create a 3D array of size (nrow+1, ncol+1, 4)
averts3d = np.full(shape_verts2d + (4,), np.nan)
averts3d[:-1, :-1, 0] = a
averts3d[:-1, 1:, 1] = a
averts3d[1:, :-1, 2] = a
averts3d[1:, 1:, 3] = a
# calculate the mean over the last axis, ignoring NaNs
averts = np.nanmean(averts3d, axis=2)
return averts
[docs]def array_at_faces_1d(a, delta):
"""
Interpolate array at cell faces of a 1d grid using linear interpolation.
Parameters
----------
a : 1d ndarray
Values at cell centers.
delta : 1d ndarray
Grid steps.
Returns
-------
afaces : 1d ndarray
Array values interpolated at cell faces, shape as input extended by 1.
"""
# extended array with ghost cells on both sides having zero values
ghost_shape = list(a.shape)
ghost_shape[0] += 2
a_ghost = np.zeros(ghost_shape, dtype=a.dtype)
# extended delta with ghost cells on both sides having zero values
delta_ghost = np.zeros(ghost_shape, dtype=a.dtype)
# fill array with ghost cells
a_ghost[1:-1] = a
a_ghost[0] = a[0]
a_ghost[-1] = a[-1]
# calculate weights
delta_ghost[1:-1] = delta
weight2 = delta_ghost[:-1] / (delta_ghost[:-1] + delta_ghost[1:])
weight1 = 1.0 - weight2
# interpolate
afaces = a_ghost[:-1] * weight1 + a_ghost[1:] * weight2
return afaces
[docs]class StructuredGrid(Grid):
"""
class for a structured model grid
Parameters
----------
delr : float or ndarray
column spacing along a row.
delc : float or ndarray
row spacing along a column.
top : float or ndarray
top elevations of cells in topmost layer
botm : float or ndarray
bottom elevations of all cells
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
----------
nlay
returns the number of model layers
nrow
returns the number of model rows
ncol
returns the number of model columns
delc
returns the delc array
delr
returns the delr array
xyedges
returns x-location points for the edges of the model grid and
y-location points for the edges of the model grid
Methods
-------
get_cell_vertices(i, j)
returns vertices for a single cell at row, column i, j.
"""
def __init__(
self,
delc=None,
delr=None,
top=None,
botm=None,
idomain=None,
lenuni=None,
crs=None,
prjfile=None,
xoff=0.0,
yoff=0.0,
angrot=0.0,
nlay=None,
nrow=None,
ncol=None,
laycbd=None,
**kwargs,
):
super().__init__(
"structured",
top=top,
botm=botm,
idomain=idomain,
lenuni=lenuni,
crs=crs,
prjfile=prjfile,
xoff=xoff,
yoff=yoff,
angrot=angrot,
**kwargs,
)
if delc is not None:
self.__nrow = len(delc)
self.__delc = delc.astype(float)
else:
self.__nrow = nrow
self.__delc = delc
if delr is not None:
self.__ncol = len(delr)
self.__delr = delr.astype(float)
else:
self.__ncol = ncol
self.__delr = delr
if top is not None:
assert self.__nrow * self.__ncol == len(np.ravel(top))
if botm is not None:
if botm.ndim == 3:
assert self.__nrow * self.__ncol == len(np.ravel(botm[0]))
if nlay is not None:
self.__nlay = nlay
else:
if laycbd is not None:
self.__nlay = len(botm) - np.count_nonzero(laycbd)
else:
self.__nlay = len(botm)
elif botm.ndim == 2:
assert botm.shape == (self.__nrow, self.__ncol)
self.__nlay = 1
else:
self.__nlay = nlay
if laycbd is not None:
self._laycbd = laycbd
else:
self._laycbd = np.zeros(self.__nlay or (), dtype=int)
####################
# Properties
####################
@property
def is_valid(self):
if self.__delc is not None and self.__delr is not None:
return True
return False
@property
def is_complete(self):
if self.__delc is not None and self.__delr is not None and super().is_complete:
return True
return False
@property
def nlay(self):
return self.__nlay
@property
def nrow(self):
return self.__nrow
@property
def ncol(self):
return self.__ncol
@property
def ncpl(self):
return self.__nrow * self.__ncol
@property
def nnodes(self):
return self.__nlay * self.__nrow * self.__ncol
@property
def nvert(self):
return (self.__nrow + 1) * (self.__ncol + 1)
@property
def iverts(self):
if self._iverts is None:
self._set_structured_iverts()
return self._iverts
@property
def verts(self):
if self._verts is None:
self._set_structured_verts()
return self._verts
@property
def shape(self):
return self.__nlay, self.__nrow, self.__ncol
@property
def extent(self):
self._copy_cache = False
xyzgrid = self.xyzvertices
self._copy_cache = True
return (
np.min(xyzgrid[0]),
np.max(xyzgrid[0]),
np.min(xyzgrid[1]),
np.max(xyzgrid[1]),
)
@property
def delc(self):
return copy.deepcopy(self.__delc)
@property
def delr(self):
return copy.deepcopy(self.__delr)
@property
def delz(self):
cache_index = "delz"
if (
cache_index not in self._cache_dict
or self._cache_dict[cache_index].out_of_date
):
delz = self.top_botm[:-1, :, :] - self.top_botm[1:, :, :]
self._cache_dict[cache_index] = CachedData(delz)
if self._copy_cache:
return self._cache_dict[cache_index].data
else:
return self._cache_dict[cache_index].data_nocopy
@property
def top_botm_withnan(self):
"""
Same as top_botm array but with NaN where idomain==0 both above and
below a cell.
"""
cache_index = "top_botm_withnan"
if (
cache_index not in self._cache_dict
or self._cache_dict[cache_index].out_of_date
):
is_inactive_above = np.full(self.top_botm.shape, True)
is_inactive_above[:-1, :, :] = self._idomain == 0
is_inactive_below = np.full(self.top_botm.shape, True)
is_inactive_below[1:, :, :] = self._idomain == 0
where_to_nan = np.logical_and(is_inactive_above, is_inactive_below)
top_botm_withnan = np.where(where_to_nan, np.nan, self.top_botm)
self._cache_dict[cache_index] = CachedData(top_botm_withnan)
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
Returns:
[]
2D array
"""
cache_index = "xyzgrid"
if (
cache_index not in self._cache_dict
or self._cache_dict[cache_index].out_of_date
):
xedge = np.concatenate(([0.0], np.add.accumulate(self.__delr)))
length_y = np.add.reduce(self.__delc)
yedge = np.concatenate(
([length_y], length_y - np.add.accumulate(self.delc))
)
xgrid, ygrid = np.meshgrid(xedge, yedge)
zgrid, zcenter = self._zcoords()
if self._has_ref_coordinates:
# transform x and y
pass
xgrid, ygrid = self.get_coords(xgrid, ygrid)
if zgrid is not None:
self._cache_dict[cache_index] = CachedData([xgrid, ygrid, zgrid])
else:
self._cache_dict[cache_index] = CachedData([xgrid, ygrid])
if self._copy_cache:
return self._cache_dict[cache_index].data
else:
return self._cache_dict[cache_index].data_nocopy
@property
def xyedges(self):
"""
Return a list of two 1D numpy arrays: one with the cell edge x
coordinate (size = ncol+1) and the other with the cell edge y
coordinate (size = nrow+1) in model space - not offset or rotated.
"""
cache_index = "xyedges"
if (
cache_index not in self._cache_dict
or self._cache_dict[cache_index].out_of_date
):
xedge = np.concatenate(([0.0], np.add.accumulate(self.__delr)))
length_y = np.add.reduce(self.__delc)
yedge = np.concatenate(
([length_y], length_y - np.add.accumulate(self.delc))
)
self._cache_dict[cache_index] = CachedData([xedge, yedge])
if self._copy_cache:
return self._cache_dict[cache_index].data
else:
return self._cache_dict[cache_index].data_nocopy
@property
def zedges(self):
"""
Return zedges for (column, row)==(0, 0).
"""
cache_index = "zedges"
if (
cache_index not in self._cache_dict
or self._cache_dict[cache_index].out_of_date
):
zedges = np.concatenate((np.array([self.top[0, 0]]), self.botm[:, 0, 0]))
self._cache_dict[cache_index] = CachedData(zedges)
if self._copy_cache:
return self._cache_dict[cache_index].data
else:
return self._cache_dict[cache_index].data_nocopy
@property
def zverts_smooth(self):
"""
Get a unique z of cell vertices using bilinear interpolation of top and
bottom elevation layers.
Returns
-------
zverts : ndarray, shape (nlay+1, nrow+1, ncol+1)
z of cell vertices. NaN values are assigned in accordance with
inactive cells defined by idomain.
"""
cache_index = "zverts_smooth"
if (
cache_index not in self._cache_dict
or self._cache_dict[cache_index].out_of_date
):
zverts_smooth = self.array_at_verts(self.top_botm)
self._cache_dict[cache_index] = CachedData(zverts_smooth)
if self._copy_cache:
return self._cache_dict[cache_index].data
else:
return self._cache_dict[cache_index].data_nocopy
@property
def xycenters(self):
"""
Return a list of two numpy one-dimensional float arrays for center x
and y coordinates in model space - not offset or rotated.
"""
cache_index = "xycenters"
if (
cache_index not in self._cache_dict
or self._cache_dict[cache_index].out_of_date
):
# get x centers
x = np.add.accumulate(self.__delr) - 0.5 * self.delr
# get y centers
Ly = np.add.reduce(self.__delc)
y = Ly - (np.add.accumulate(self.__delc) - 0.5 * self.__delc)
# store in cache
self._cache_dict[cache_index] = CachedData([x, y])
if self._copy_cache:
return self._cache_dict[cache_index].data
else:
return self._cache_dict[cache_index].data_nocopy
@property
def xyzcellcenters(self):
"""
Return a list of three numpy float arrays: two two-dimensional arrays
for center x and y coordinates, and one three-dimensional array for
center z coordinates. Coordinates are given in real-world coordinates.
"""
cache_index = "cellcenters"
if (
cache_index not in self._cache_dict
or self._cache_dict[cache_index].out_of_date
):
# get x centers
x = np.add.accumulate(self.__delr) - 0.5 * self.delr
# get y centers
Ly = np.add.reduce(self.__delc)
y = Ly - (np.add.accumulate(self.__delc) - 0.5 * self.__delc)
x_mesh, y_mesh = np.meshgrid(x, y)
if self.__nlay is not None:
# get z centers
z = np.empty((self.__nlay, self.__nrow, self.__ncol))
z[0, :, :] = (self._top[:, :] + self._botm[0, :, :]) / 2.0
ibs = np.arange(self.__nlay)
quasi3d = [cbd != 0 for cbd in self._laycbd]
if np.any(quasi3d):
ibs[1:] = ibs[1:] + np.cumsum(quasi3d)[: self.__nlay - 1]
for l, ib in enumerate(ibs[1:], 1):
z[l, :, :] = (self._botm[ib - 1, :, :] + self._botm[ib, :, :]) / 2.0
else:
z = None
if self._has_ref_coordinates:
# transform x and y
x_mesh, y_mesh = self.get_coords(x_mesh, y_mesh)
# store in cache
self._cache_dict[cache_index] = CachedData([x_mesh, y_mesh, z])
if self._copy_cache:
return self._cache_dict[cache_index].data
else:
return self._cache_dict[cache_index].data_nocopy
@property
def grid_lines(self):
"""
Get the grid lines as a list
"""
# get edges initially in model coordinates
use_ref_coords = self.use_ref_coords
self.use_ref_coords = False
xyedges = self.xyedges
self.use_ref_coords = use_ref_coords
xmin = xyedges[0][0]
xmax = xyedges[0][-1]
ymin = xyedges[1][-1]
ymax = xyedges[1][0]
lines = []
# Vertical lines
for j in range(self.ncol + 1):
x0 = xyedges[0][j]
x1 = x0
y0 = ymin
y1 = ymax
lines.append([(x0, y0), (x1, y1)])
# horizontal lines
for i in range(self.nrow + 1):
x0 = xmin
x1 = xmax
y0 = xyedges[1][i]
y1 = y0
lines.append([(x0, y0), (x1, y1)])
if self._has_ref_coordinates:
lines_trans = []
for ln in lines:
lines_trans.append([self.get_coords(*ln[0]), self.get_coords(*ln[1])])
return lines_trans
return lines
@property
def is_regular_x(self):
"""
Test whether the grid spacing is regular in the x direction.
"""
cache_index = "is_regular_x"
if (
cache_index not in self._cache_dict
or self._cache_dict[cache_index].out_of_date
):
# relative tolerance to use in test
rel_tol = 1.0e-5
# regularity test in x direction
rel_diff_x = (self.__delr - self.__delr[0]) / self.__delr[0]
is_regular_x = np.count_nonzero(np.abs(rel_diff_x) > rel_tol) == 0
self._cache_dict[cache_index] = CachedData(is_regular_x)
if self._copy_cache:
return self._cache_dict[cache_index].data
else:
return self._cache_dict[cache_index].data_nocopy
@property
def is_regular_y(self):
"""
Test whether the grid spacing is regular in the y direction.
"""
cache_index = "is_regular_y"
if (
cache_index not in self._cache_dict
or self._cache_dict[cache_index].out_of_date
):
# relative tolerance to use in test
rel_tol = 1.0e-5
# regularity test in y direction
rel_diff_y = (self.__delc - self.__delc[0]) / self.__delc[0]
is_regular_y = np.count_nonzero(np.abs(rel_diff_y) > rel_tol) == 0
self._cache_dict[cache_index] = CachedData(is_regular_y)
if self._copy_cache:
return self._cache_dict[cache_index].data
else:
return self._cache_dict[cache_index].data_nocopy
@property
def is_regular_z(self):
"""
Test if the grid spacing is regular in z direction.
"""
cache_index = "is_regular_z"
if (
cache_index not in self._cache_dict
or self._cache_dict[cache_index].out_of_date
):
# relative tolerance to use in test
rel_tol = 1.0e-5
# regularity test in z direction
rel_diff_thick0 = (self.delz[0, :, :] - self.delz[0, 0, 0]) / self.delz[
0, 0, 0
]
failed = np.abs(rel_diff_thick0) > rel_tol
_is_regular_z = np.count_nonzero(failed) == 0
for k in range(1, self.nlay):
rel_diff_zk = (self.delz[k, :, :] - self.delz[0, :, :]) / self.delz[
0, :, :
]
failed = np.abs(rel_diff_zk) > rel_tol
_is_regular_z = _is_regular_z and np.count_nonzero(failed) == 0
self._cache_dict[cache_index] = CachedData(_is_regular_z)
if self._copy_cache:
return self._cache_dict[cache_index].data
else:
return self._cache_dict[cache_index].data_nocopy
@property
def is_regular_xy(self):
"""
Test if the grid spacing is regular and equal in x and y directions.
"""
cache_index = "is_regular_xy"
if (
cache_index not in self._cache_dict
or self._cache_dict[cache_index].out_of_date
):
# relative tolerance to use in test
rel_tol = 1.0e-5
# test if the first delta is equal in x and z
rel_diff_0 = (self.__delc[0] - self.__delr[0]) / self.__delr[0]
first_equal = np.abs(rel_diff_0) <= rel_tol
# combine with regularity tests in x and z directions
is_regular_xy = first_equal and self.is_regular_x and self.is_regular_y
self._cache_dict[cache_index] = CachedData(is_regular_xy)
if self._copy_cache:
return self._cache_dict[cache_index].data
else:
return self._cache_dict[cache_index].data_nocopy
@property
def is_regular_xz(self):
"""
Test if the grid spacing is regular and equal in x and z directions.
"""
cache_index = "is_regular_xz"
if (
cache_index not in self._cache_dict
or self._cache_dict[cache_index].out_of_date
):
# relative tolerance to use in test
rel_tol = 1.0e-5
# test if the first delta is equal in x and z
rel_diff_0 = (self.delz[0, 0, 0] - self.__delr[0]) / self.__delr[0]
first_equal = np.abs(rel_diff_0) <= rel_tol
# combine with regularity tests in x and z directions
is_regular_xz = first_equal and self.is_regular_x and self.is_regular_z
self._cache_dict[cache_index] = CachedData(is_regular_xz)
if self._copy_cache:
return self._cache_dict[cache_index].data
else:
return self._cache_dict[cache_index].data_nocopy
@property
def is_regular_yz(self):
"""
Test if the grid spacing is regular and equal in y and z directions.
"""
cache_index = "is_regular_yz"
if (
cache_index not in self._cache_dict
or self._cache_dict[cache_index].out_of_date
):
# relative tolerance to use in test
rel_tol = 1.0e-5
# test if the first delta is equal in y and z
rel_diff_0 = (self.delz[0, 0, 0] - self.__delc[0]) / self.__delc[0]
first_equal = np.abs(rel_diff_0) <= rel_tol
# combine with regularity tests in x and y directions
is_regular_yz = first_equal and self.is_regular_y and self.is_regular_z
self._cache_dict[cache_index] = CachedData(is_regular_yz)
if self._copy_cache:
return self._cache_dict[cache_index].data
else:
return self._cache_dict[cache_index].data_nocopy
@property
def is_regular(self):
"""
Test if the grid spacing is regular and equal in x, y and z directions.
"""
cache_index = "is_regular"
if (
cache_index not in self._cache_dict
or self._cache_dict[cache_index].out_of_date
):
# relative tolerance to use in test
rel_tol = 1.0e-5
# test if the first delta is equal in x and z
rel_diff_0 = (self.delz[0, 0, 0] - self.__delr[0]) / self.__delr[0]
first_equal = np.abs(rel_diff_0) <= rel_tol
# combine with regularity tests in x, y and z directions
is_regular = first_equal and self.is_regular_z and self.is_regular_xy
self._cache_dict[cache_index] = CachedData(is_regular)
if self._copy_cache:
return self._cache_dict[cache_index].data
else:
return self._cache_dict[cache_index].data_nocopy
@property
def is_rectilinear(self):
"""
Test whether the grid is rectilinear (it is always so in the x and
y directions, but not necessarily in the z direction).
"""
cache_index = "is_rectilinear"
if (
cache_index not in self._cache_dict
or self._cache_dict[cache_index].out_of_date
):
# relative tolerance to use in test
rel_tol = 1.0e-5
# rectilinearity test in z direction
is_rect_z = True
for k in range(self.nlay):
rel_diff_zk = (self.delz[k, :, :] - self.delz[k, 0, 0]) / self.delz[
k, 0, 0
]
failed = np.abs(rel_diff_zk) > rel_tol
is_rect_z = is_rect_z and np.count_nonzero(failed) == 0
self._cache_dict[cache_index] = CachedData(is_rect_z)
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 = (self.xvertices, self.yvertices)
return 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 = [[list(zip(*i))] for i in zip(*self.cross_section_vertices)]
self._cache_dict[cache_index] = CachedData(polys)
else:
polys = self._cache_dict[cache_index].data_nocopy
gdf = super().to_geodataframe(polys)
gdf["row"] = sorted(list(range(1, self.nrow + 1)) * self.ncol)
gdf["col"] = list(range(1, self.ncol + 1)) * self.nrow
if self.idomain is not None:
active = np.sum(
self.idomain.reshape(
(-1, 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 super().is_complete:
return StructuredGrid(
delc=self.delc * factor,
delr=self.delr * factor,
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")
###############
### Methods ###
###############
[docs] def neighbors(self, *args, **kwargs):
"""
Method to get nearest neighbors for a cell
Parameters
----------
*args
lay (int), row (int), column (int)
or
node (int)
**kwargs
k : int
layer number
i : int
row number
j : int
column number
as_nodes : bool
flag to return neighbors as node numbers
method : str
"rook" for shared edge neighbors (default) "queen" for shared
vertex neighbors (for flow accumulation calculations)
reset : bool
flag to re-calculate neighbors, default is False
Returns
-------
list of neighboring cells
"""
nn = None
as_nodes = kwargs.pop("as_nodes", False)
memhog = kwargs.pop("fast", False)
if kwargs:
if "node" in kwargs:
nn = kwargs.pop("node")
as_nodes = True
else:
k = kwargs.pop("k", 0)
i = kwargs.pop("i", None)
j = kwargs.pop("j", None)
if i is None or j is None:
pass
else:
nn = self.get_node([(k, i, j)])[0]
if len(args) > 0:
if len(args) == 1:
nn = args[0]
as_nodes = True
elif len(args) == 2:
k = 0
i, j = args[0:2]
else:
k, i, j = args[0:3]
if nn is None:
nn = self.get_node([(k, i, j)])[0]
else:
as_nodes = True
if memhog:
neighbors = self._fast_neighbors(**kwargs)
else:
neighbors = super().neighbors(nn, **kwargs)
if not as_nodes:
neighbors = self.get_lrc(neighbors)
return neighbors
def _fast_neighbors(self, **kwargs):
"""
Memory intensive and not elegent in any sense, but very fast method to get
neighbors for Structured Grids
Returns
-------
dict
"""
from collections import defaultdict
if self._neighbors is None or kwargs.pop("reset", False):
nrow = self.nrow
ncol = self.ncol
ncpl = nrow * ncol
arr = np.arange(ncpl, dtype=np.int32).reshape((nrow, ncol))
# general case
arr2 = arr[1:-1, 1:-1].ravel()
neighs2 = np.empty((4, (nrow - 2) * (ncol - 2)), dtype=np.int32)
neighs2[0] = arr2 - ncol # u
neighs2[1] = arr2 + 1 # r
neighs2[2] = arr2 + ncol # d
neighs2[3] = arr2 - 1 # l
neighs2 = neighs2.T.tolist()
neighbors = {v: neighs2[ix] for ix, v in enumerate(arr2)}
# ecase 1
arr2 = arr[0, 1:-1].ravel()
neighs2 = np.empty((3, ncol - 2), dtype=np.int32)
# no up
neighs2[0] = arr2 + 1 # r
neighs2[1] = arr2 + ncol # d
neighs2[2] = arr2 - 1 # l
neighs2 = neighs2.T.tolist()
for ix, v in enumerate(arr2):
neighbors[v] = neighs2[ix]
# ecase 2
arr2 = arr[-1, 1:-1].ravel()
neighs2 = np.empty((3, ncol - 2), dtype=np.int32)
neighs2[0] = arr2 - ncol # u
neighs2[1] = arr2 + 1 # r
# no down
neighs2[2] = arr2 - 1 # l
neighs2 = neighs2.T.tolist()
for ix, v in enumerate(arr2):
neighbors[v] = neighs2[ix]
# ecase 3
arr2 = arr[1:-1, 0].ravel()
neighs2 = np.empty((3, nrow - 2), dtype=np.int32)
neighs2[0] = arr2 - ncol # u
neighs2[1] = arr2 + 1 # r
neighs2[2] = arr2 + ncol # d
# no left
neighs2 = neighs2.T.tolist()
for ix, v in enumerate(arr2):
neighbors[v] = neighs2[ix]
# ecase 4
arr2 = arr[1:-1, -1].ravel()
neighs2 = np.empty((3, nrow - 2), dtype=np.int32)
neighs2[0] = arr2 - ncol # u
# no right
neighs2[1] = arr2 + ncol # d
neighs2[2] = arr2 - 1 # l
neighs2 = neighs2.T.tolist()
for ix, v in enumerate(arr2):
neighbors[v] = neighs2[ix]
del arr
# corners
neighbors[0] = [1, ncol] # r, d
neighbors[ncol - 1] = [(ncol - 1) + ncol, ncol - 2] # d, l
neighbors[ncpl - ncol] = [(ncpl - ncol) - ncol, (ncpl - ncol) + 1] # u, r
neighbors[ncpl - 1] = [(ncpl - 1) - ncol, ncpl - 2] # u, l
self._neighbors = neighbors
return self._neighbors
[docs] def intersect(self, x, y, z=None, local=False, forgive=False):
"""
Get the row and column of a point with coordinates x and y
When the point is on the edge of two cells, the cell with the lowest
row or column 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 layer,
row, column) if supplied
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
-------
row : int or ndarray
The row number(s). Returns int for scalar input, ndarray for array input.
col : int or ndarray
The column 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")
# transform x and y to local coordinates
if not local:
x, y = self.get_local_coords(x, y)
# get the cell edges in local coordinates
xe, ye = self.xyedges
# Vectorized row/col calculation
n_points = len(x)
rows = np.full(n_points, np.nan if forgive else -1, dtype=float)
cols = np.full(n_points, np.nan if forgive else -1, dtype=float)
for i in range(n_points):
xi, yi = x[i], y[i]
# Find column
xcomp = xi > xe
if np.all(xcomp) or not np.any(xcomp):
if forgive:
cols[i] = np.nan
else:
raise ValueError(
f"x, y point given is outside of the model area: ({xi}, {yi})"
)
else:
cols[i] = np.asarray(xcomp).nonzero()[0][-1]
# Find row
ycomp = yi < ye
if np.all(ycomp) or not np.any(ycomp):
if forgive:
rows[i] = np.nan
else:
raise ValueError(
f"x, y point given is outside of the model area: ({xi}, {yi})"
)
else:
rows[i] = np.asarray(ycomp).nonzero()[0][-1]
# If either row or col is NaN, set both to NaN
invalid_mask = np.isnan(rows) | np.isnan(cols)
rows[invalid_mask] = np.nan
cols[invalid_mask] = np.nan
# Convert to int where valid
valid_mask = ~invalid_mask
if np.any(valid_mask):
rows[valid_mask] = rows[valid_mask].astype(int)
cols[valid_mask] = cols[valid_mask].astype(int)
if z is None:
# Return results
if is_scalar_input:
row, col = rows[0], cols[0]
if not np.isnan(row) and not np.isnan(col):
row, col = int(row), int(col)
return row, col
else:
return rows.astype(int) if np.all(valid_mask) else rows, cols.astype(
int
) if np.all(valid_mask) else cols
# Handle z-coordinate
lays = np.full(n_points, np.nan if forgive else -1, dtype=float)
for i in range(n_points):
if np.isnan(rows[i]) or np.isnan(cols[i]):
continue
row, col = int(rows[i]), int(cols[i])
zi = z[i]
for layer in range(self.__nlay):
if (
self.top_botm[layer, row, col]
>= zi
>= self.top_botm[layer + 1, row, col]
):
lays[i] = layer
break
if np.isnan(lays[i]) and not forgive:
raise ValueError(
f"point given is outside the model area: ({x[i]}, {y[i]}, {zi})"
)
# Return results
if is_scalar_input:
lay, row, col = lays[0], rows[0], cols[0]
if not np.isnan(lay):
lay, row, col = int(lay), int(row), int(col)
else:
# When x,y are out of bounds: lay=None, row/col keep NaN
lay = None
# row and col already have their NaN values
return lay, row, col
else:
valid_3d = ~np.isnan(lays) & ~np.isnan(rows) & ~np.isnan(cols)
return (
lays.astype(int) if np.all(valid_3d) else lays,
rows.astype(int) if np.all(valid_3d) else rows,
cols.astype(int) if np.all(valid_3d) else cols,
)
def _cell_vert_list(self, i, j):
"""Get vertices for a single cell or sequence of i, j locations."""
self._copy_cache = False
pts = []
xgrid, ygrid = self.xvertices, self.yvertices
pts.append([xgrid[i, j], ygrid[i, j]])
pts.append([xgrid[i + 1, j], ygrid[i + 1, j]])
pts.append([xgrid[i + 1, j + 1], ygrid[i + 1, j + 1]])
pts.append([xgrid[i, j + 1], ygrid[i, j + 1]])
pts.append([xgrid[i, j], ygrid[i, j]])
self._copy_cache = True
if np.isscalar(i):
return pts
else:
vrts = np.array(pts).transpose([2, 0, 1])
return [v.tolist() for v in vrts]
@property
def top_botm(self):
new_top = np.expand_dims(self._top, 0)
return np.concatenate((new_top, self._botm), axis=0)
[docs] def get_cell_vertices(self, *args, **kwargs):
"""
Get a set of cell vertices for a single cell.
Parameters
----------
cellid : int or tuple, optional
Cell identifier. Can be:
- node number (int)
- (row, col) tuple
- (layer, row, col) tuple (layer is ignored, vertices are 2D)
node : int, optional
Node index, mutually exclusive with cellid, i, and j
i, j : int, optional
Row and column index, mutually exclusive with cellid and node
Returns
-------
list
list of (x, y) cell vertex coordinates
Examples
--------
>>> import flopy
>>> import numpy as np
>>> delr, delc = np.array([10.0] * 3), np.array([10.0] * 4)
>>> sg = flopy.discretization.StructuredGrid(delr=delr, delc=delc)
>>> sg.get_cell_vertices(node=0)
[(0.0, 40.0), (10.0, 40.0), (10.0, 30.0), (0.0, 30.0)]
>>> sg.get_cell_vertices(3, 0)
[(0.0, 10.0), (10.0, 10.0), (10.0, 0.0), (0.0, 0.0)]
>>> sg.get_cell_vertices((3, 0))
[(0.0, 10.0), (10.0, 10.0), (10.0, 0.0), (0.0, 0.0)]
>>> sg.get_cell_vertices(cellid=(0, 3, 0))
[(0.0, 10.0), (10.0, 10.0), (10.0, 0.0), (0.0, 0.0)]
"""
if kwargs:
if args:
raise TypeError("mixed positional and keyword arguments not supported")
elif "cellid" in kwargs:
cellid = kwargs.pop("cellid")
# Handle cellid as int, tuple of 2, or tuple of 3
if isinstance(cellid, (tuple, list)):
if len(cellid) == 2:
i, j = cellid
elif len(cellid) == 3:
_, i, j = cellid # ignore layer
else:
raise ValueError(
f"cellid tuple must have 2 or 3 elements, got {len(cellid)}"
)
else:
_, i, j = self.get_lrc(cellid)[0]
elif "node" in kwargs:
_, i, j = self.get_lrc(kwargs.pop("node"))[0]
elif "i" in kwargs and "j" in kwargs:
i = kwargs.pop("i")
j = kwargs.pop("j")
else:
raise TypeError(
"expected cellid, node, or i and j as keyword arguments"
)
if kwargs:
unused = ", ".join(kwargs.keys())
raise TypeError(f"unused keyword arguments: {unused}")
elif len(args) == 0:
raise TypeError("expected one or more arguments")
elif len(args) == 1:
# Single arg could be node number or (row, col) or (layer, row, col) tuple
arg = args[0]
if isinstance(arg, (tuple, list)):
if len(arg) == 2:
i, j = arg
elif len(arg) == 3:
# (layer, row, col) - ignore layer
_, i, j = arg
else:
raise ValueError(
f"cellid tuple must have 2 or 3 elements, got {len(arg)}"
)
else:
# Node number
_, i, j = self.get_lrc(arg)[0]
elif len(args) == 2:
i, j = args
else:
raise TypeError("too many arguments")
self._copy_cache = False
cell_verts = [
(self.xvertices[i, j], self.yvertices[i, j]),
(self.xvertices[i, j + 1], self.yvertices[i, j + 1]),
(self.xvertices[i + 1, j + 1], self.yvertices[i + 1, j + 1]),
(self.xvertices[i + 1, j], self.yvertices[i + 1, j]),
]
self._copy_cache = True
return cell_verts
[docs] def get_lrc(self, nodes):
"""
Get layer, row, column from a list of zero-based
MODFLOW node numbers.
Parameters
----------
nodes : int, list or array_like
Zero-based node number
Returns
-------
list
list of tuples containing the layer (k), row (i),
and column (j) for each node in the input list
Examples
--------
>>> import flopy
>>> sg = flopy.discretization.StructuredGrid(nlay=20, nrow=30, ncol=40)
>>> sg.get_lrc(100)
[(0, 2, 20)]
>>> sg.get_lrc([100, 1000, 10_000])
[(0, 2, 20), (0, 25, 0), (8, 10, 0)]
"""
if isinstance(nodes, int):
nodes = [nodes]
shape = self.shape
if shape[0] is None:
shape = tuple(dim or 1 for dim in shape)
return list(zip(*np.unravel_index(nodes, shape)))
[docs] def get_node(self, lrc_list, cellids=None, node2d=False):
"""
Get node number from a list of zero-based MODFLOW
layer, row, column tuples.
Parameters
----------
lrc_list : tuple of int or list of tuple of int
Zero-based layer, row, column tuples
.. deprecated:: 3.10
This parameter is deprecated and will be
removed in FloPy 3.12. Use cellids instead.
cellids : tuple of int or list of tuple of int
Zero-based layer, row, column tuples
node2d : bool, optional
If True, return 2D node numbers (ignore layer).
If False (default), return 3D node numbers.
Returns
-------
list
list of MODFLOW nodes for each layer (k), row (i),
and column (j) tuple in the input list
Examples
--------
>>> import flopy
>>> sg = flopy.discretization.StructuredGrid(nlay=20, nrow=30, ncol=40)
>>> sg.get_node((0, 2, 20))
[100]
>>> sg.get_node([(0, 2, 20), (0, 25, 0), (8, 10, 0)])
[100, 1000, 10000]
"""
if lrc_list is not None:
if cellids is not None:
raise TypeError("lrc_list and cellids are mutually exclusive")
cellids = lrc_list
if cellids is None:
raise ValueError("Expected a value for cellids")
if not isinstance(cellids, list):
cellids = [cellids]
if node2d:
rc_list = [(row, col) for lay, row, col in cellids]
multi_index = tuple(np.array(rc_list).T)
shape = (self.nrow, self.ncol)
else:
multi_index = tuple(np.array(cellids).T)
shape = self.shape
if shape[0] is None:
shape = tuple(dim or 1 for dim in shape)
return np.ravel_multi_index(multi_index, shape).tolist()
[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)
[docs] def array_at_verts_basic(self, a):
"""
Computes values at cell vertices using neighbor averaging.
Parameters
----------
a : ndarray
Array values at cell centers.
Returns
-------
averts : ndarray
Array values at cell vertices, shape
(a.shape[0]+1, a.shape[1]+1, a.shape[2]+1). NaN values are assigned
in accordance with inactive cells defined by idomain.
"""
assert a.ndim == 3
shape_verts = (a.shape[0] + 1, a.shape[1] + 1, a.shape[2] + 1)
# set to NaN where idomain==0
a[self._idomain == 0] = np.nan
# create a 4D array of size (nlay+1, nrow+1, ncol+1, 8)
averts4d = np.full(shape_verts + (8,), np.nan)
averts4d[:-1, :-1, :-1, 0] = a
averts4d[:-1, :-1, 1:, 1] = a
averts4d[:-1, 1:, :-1, 2] = a
averts4d[:-1, 1:, 1:, 3] = a
averts4d[1:, :-1, :-1, 4] = a
averts4d[1:, :-1, 1:, 5] = a
averts4d[1:, 1:, :-1, 6] = a
averts4d[1:, 1:, 1:, 7] = a
# calculate the mean over the last axis, ignoring NaNs
averts = np.nanmean(averts4d, axis=3)
return averts
[docs] def array_at_verts(self, a):
"""
Interpolate array values at cell vertices.
Parameters
----------
a : ndarray
Array values. Allowed shapes are: (nlay, nrow, ncol),
(nlay, nrow, ncol+1), (nlay, nrow+1, ncol) and
(nlay+1, nrow, ncol).
* When the shape is (nlay, nrow, ncol), input values are
considered at cell centers, and output values are computed by
trilinear interpolation.
* When the shape is extended in one direction, input values are
considered at the center of cell faces in this direction, and
output values are computed by bilinear interpolation in planes
defined by these cell faces.
Returns
-------
averts : ndarray
Array values interpolated at cell vertices, shape
(nlay+1, nrow+1, ncol+1).
Notes
-----
* Output values are smooth (continuous) even if top elevations or
bottom elevations are not constant across layers (i.e., in this
case, vertices of neighboring cells are implicitly merged).
* NaN values are assigned in accordance with inactive cells defined
by idomain.
"""
import scipy.interpolate as interp
# define shapes
shape_ext_x = (self.nlay, self.nrow, self.ncol + 1)
shape_ext_y = (self.nlay, self.nrow + 1, self.ncol)
shape_ext_z = (self.nlay + 1, self.nrow, self.ncol)
shape_verts = (self.nlay + 1, self.nrow + 1, self.ncol + 1)
# get inactive cells
if self._idomain is not None:
inactive = self._idomain == 0
# get local x and y cell center coordinates (1d arrays)
xcenters, ycenters = self.xycenters
# get z center coordinates: make the grid rectilinear if it is not,
# in order to always use RegularGridInterpolator; in most cases this
# will give better results than with the non-structured interpolator
# LinearNDInterpolator (in addition, it will run faster)
zcenters = self.zcellcenters
if self._idomain is not None:
zcenters = np.where(inactive, np.nan, zcenters)
if not self.is_rectilinear or np.count_nonzero(np.isnan(zcenters)) != 0:
zedges = np.nanmean(self.top_botm_withnan, axis=(1, 2))
else:
zedges = self.top_botm_withnan[:, 0, 0]
zcenters = 0.5 * (zedges[1:] + zedges[:-1])
# test grid regularity in z
rel_tol = 1.0e-5
delz = np.diff(zedges)
rel_diff = (delz - delz[0]) / delz[0]
is_regular_z = np.count_nonzero(np.abs(rel_diff) > rel_tol) == 0
# test equality of first grid spacing in x and z, and in y and z
first_equal_xz = np.abs(self.__delr[0] - delz[0]) / delz[0] <= rel_tol
first_equal_yz = np.abs(self.__delc[0] - delz[0]) / delz[0] <= rel_tol
# get output coordinates (i.e. vertices)
xedges, yedges = self.xyedges
xedges = xedges.reshape((1, 1, self.ncol + 1))
xoutput = xedges * np.ones(shape_verts)
yedges = yedges.reshape((1, self.nrow + 1, 1))
youtput = yedges * np.ones(shape_verts)
zoutput = zedges.reshape((self.nlay + 1, 1, 1))
zoutput = zoutput * np.ones(shape_verts)
# indicator of whether basic interpolation is used or not
basic = False
if a.shape == self.shape:
# set array to NaN where inactive
if self._idomain is not None:
inactive = self._idomain == 0
a = np.where(inactive, np.nan, a)
# perform basic interpolation (this will be useful in all cases)
averts_basic = self.array_at_verts_basic(a)
if self.is_regular_xy and is_regular_z and first_equal_xz:
# in this case, basic interpolation is the correct one
averts = averts_basic
basic = True
else:
if self.nlay == 1:
# in this case we need a 2d interpolation in the x, y plane
# flip y coordinates because RegularGridInterpolator
# requires increasing input coordinates
xyinput = (np.flip(ycenters), xcenters)
a = np.squeeze(np.flip(a, axis=[1]))
# interpolate
interp_func = interp.RegularGridInterpolator(
xyinput, a, bounds_error=False, fill_value=np.nan
)
xyoutput = np.empty((youtput[0, :, :].size, 2))
xyoutput[:, 0] = youtput[0, :, :].ravel()
xyoutput[:, 1] = xoutput[0, :, :].ravel()
averts2d = interp_func(xyoutput)
averts2d = averts2d.reshape((1, self.nrow + 1, self.ncol + 1))
averts = averts2d * np.ones(shape_verts)
elif self.nrow == 1:
# in this case we need a 2d interpolation in the x, z plane
# flip z coordinates because RegularGridInterpolator
# requires increasing input coordinates
xzinput = (np.flip(zcenters), xcenters)
a = np.squeeze(np.flip(a, axis=[0]))
# interpolate
interp_func = interp.RegularGridInterpolator(
xzinput, a, bounds_error=False, fill_value=np.nan
)
xzoutput = np.empty((zoutput[:, 0, :].size, 2))
xzoutput[:, 0] = zoutput[:, 0, :].ravel()
xzoutput[:, 1] = xoutput[:, 0, :].ravel()
averts2d = interp_func(xzoutput)
averts2d = averts2d.reshape((self.nlay + 1, 1, self.ncol + 1))
averts = averts2d * np.ones(shape_verts)
elif self.ncol == 1:
# in this case we need a 2d interpolation in the y, z plane
# flip y and z coordinates because RegularGridInterpolator
# requires increasing input coordinates
yzinput = (np.flip(zcenters), np.flip(ycenters))
a = np.squeeze(np.flip(a, axis=[0, 1]))
# interpolate
interp_func = interp.RegularGridInterpolator(
yzinput, a, bounds_error=False, fill_value=np.nan
)
yzoutput = np.empty((zoutput[:, :, 0].size, 2))
yzoutput[:, 0] = zoutput[:, :, 0].ravel()
yzoutput[:, 1] = youtput[:, :, 0].ravel()
averts2d = interp_func(yzoutput)
averts2d = averts2d.reshape((self.nlay + 1, self.nrow + 1, 1))
averts = averts2d * np.ones(shape_verts)
else:
# 3d interpolation
# flip y and z coordinates because RegularGridInterpolator
# requires increasing input coordinates
xyzinput = (np.flip(zcenters), np.flip(ycenters), xcenters)
a = np.flip(a, axis=[0, 1])
# interpolate
interp_func = interp.RegularGridInterpolator(
xyzinput, a, bounds_error=False, fill_value=np.nan
)
xyzoutput = np.empty((zoutput.size, 3))
xyzoutput[:, 0] = zoutput.ravel()
xyzoutput[:, 1] = youtput.ravel()
xyzoutput[:, 2] = xoutput.ravel()
averts = interp_func(xyzoutput)
averts = averts.reshape(shape_verts)
elif a.shape == shape_ext_x:
# set array to NaN where inactive on both side
if self._idomain is not None:
inactive_ext_x = np.full(shape_ext_x, True)
inactive_ext_x[:, :, :-1] = inactive
inactive_ext_x[:, :, 1:] = np.logical_and(
inactive_ext_x[:, :, 1:], inactive
)
a = np.where(inactive_ext_x, np.nan, a)
averts = np.empty(shape_verts, dtype=a.dtype)
averts_basic = np.empty(shape_verts, dtype=a.dtype)
for j in range(self.ncol + 1):
# perform basic interpolation (will be useful in all cases)
averts_basic[:, :, j] = array_at_verts_basic2d(a[:, :, j])
if self.is_regular_y and is_regular_z and first_equal_yz:
# in this case, basic interpolation is the correct one
averts2d = averts_basic[:, :, j]
basic = True
else:
if self.nlay == 1:
# in this case we need a 1d interpolation along y
averts1d = array_at_faces_1d(a[0, :, j], self.__delc)
averts2d = averts1d.reshape((1, self.nrow + 1))
averts2d = averts2d * np.ones((2, self.nrow + 1))
elif self.nrow == 1:
# in this case we need a 1d interpolation along z
delz1d = np.abs(np.diff(self.zverts_smooth[:, 0, j]))
averts1d = array_at_faces_1d(a[:, 0, j], delz1d)
averts2d = averts1d.reshape((self.nlay + 1, 1))
averts2d = averts2d * np.ones((self.nlay + 1, 2))
else:
# 2d interpolation
# flip y and z coordinates because
# RegularGridInterpolator requires increasing input
# coordinates
yzinput = (np.flip(zcenters), np.flip(ycenters))
a2d = np.flip(a[:, :, j], axis=[0, 1])
interp_func = interp.RegularGridInterpolator(
yzinput, a2d, bounds_error=False, fill_value=np.nan
)
yzoutput = np.empty((zoutput[:, :, j].size, 2))
yzoutput[:, 0] = zoutput[:, :, j].ravel()
yzoutput[:, 1] = youtput[:, :, j].ravel()
averts2d = interp_func(yzoutput)
averts2d = averts2d.reshape(zoutput[:, :, j].shape)
averts[:, :, j] = averts2d
elif a.shape == shape_ext_y:
# set array to NaN where inactive on both side
if self._idomain is not None:
inactive_ext_y = np.full(shape_ext_y, True)
inactive_ext_y[:, :-1, :] = inactive
inactive_ext_y[:, 1:, :] = np.logical_and(
inactive_ext_y[:, 1:, :], inactive
)
a = np.where(inactive_ext_y, np.nan, a)
averts = np.empty(shape_verts, dtype=a.dtype)
averts_basic = np.empty(shape_verts, dtype=a.dtype)
for i in range(self.nrow + 1):
# perform basic interpolation (will be useful in all cases)
averts_basic[:, i, :] = array_at_verts_basic2d(a[:, i, :])
if self.is_regular_x and is_regular_z and first_equal_xz:
# in this case, basic interpolation is the correct one
averts2d = averts_basic[:, i, :]
basic = True
else:
if self.nlay == 1:
# in this case we need a 1d interpolation along x
averts1d = array_at_faces_1d(a[0, i, :], self.__delr)
averts2d = averts1d.reshape((1, self.ncol + 1))
averts2d = averts2d * np.ones((2, self.ncol + 1))
elif self.ncol == 1:
# in this case we need a 1d interpolation along z
delz1d = np.abs(np.diff(self.zverts_smooth[:, i, 0]))
averts1d = array_at_faces_1d(a[:, i, 0], delz1d)
averts2d = averts1d.reshape((self.nlay + 1, 1))
averts2d = averts2d * np.ones((self.nlay + 1, 2))
else:
# 2d interpolation
# flip z coordinates because RegularGridInterpolator
# requires increasing input coordinates
xzinput = (np.flip(zcenters), xcenters)
a2d = np.flip(a[:, i, :], axis=[0])
interp_func = interp.RegularGridInterpolator(
xzinput, a2d, bounds_error=False, fill_value=np.nan
)
xzoutput = np.empty((zoutput[:, i, :].size, 2))
xzoutput[:, 0] = zoutput[:, i, :].ravel()
xzoutput[:, 1] = xoutput[:, i, :].ravel()
averts2d = interp_func(xzoutput)
averts2d = averts2d.reshape(zoutput[:, i, :].shape)
averts[:, i, :] = averts2d
elif a.shape == shape_ext_z:
# set array to NaN where inactive on both side
if self._idomain is not None:
inactive_ext_z = np.full(shape_ext_z, True)
inactive_ext_z[:-1, :, :] = inactive
inactive_ext_z[1:, :, :] = np.logical_and(
inactive_ext_z[1:, :, :], inactive
)
a = np.where(inactive_ext_z, np.nan, a)
averts = np.empty(shape_verts, dtype=a.dtype)
averts_basic = np.empty(shape_verts, dtype=a.dtype)
for k in range(self.nlay + 1):
# perform basic interpolation (will be useful in all cases)
averts_basic[k, :, :] = array_at_verts_basic2d(a[k, :, :])
if self.is_regular_xy:
# in this case, basic interpolation is the correct one
averts2d = averts_basic[k, :, :]
basic = True
else:
if self.nrow == 1:
# in this case we need a 1d interpolation along x
averts1d = array_at_faces_1d(a[k, 0, :], self.__delr)
averts2d = averts1d.reshape((1, self.ncol + 1))
averts2d = averts2d * np.ones((2, self.ncol + 1))
elif self.ncol == 1:
# in this case we need a 1d interpolation along y
averts1d = array_at_faces_1d(a[k, :, 0], self.__delc)
averts2d = averts1d.reshape((self.nrow + 1, 1))
averts2d = averts2d * np.ones((self.nrow + 1, 2))
else:
# 2d interpolation
# flip y coordinates because RegularGridInterpolator
# requires increasing input coordinates
xyinput = (np.flip(ycenters), xcenters)
a2d = np.flip(a[k, :, :], axis=[0])
interp_func = interp.RegularGridInterpolator(
xyinput, a2d, bounds_error=False, fill_value=np.nan
)
xyoutput = np.empty((youtput[k, :, :].size, 2))
xyoutput[:, 0] = youtput[k, :, :].ravel()
xyoutput[:, 1] = xoutput[k, :, :].ravel()
averts2d = interp_func(xyoutput)
averts2d = averts2d.reshape(youtput[k, :, :].shape)
averts[k, :, :] = averts2d
if not basic:
# use basic interpolation for remaining NaNs at boundaries
where_nan = np.isnan(averts)
averts[where_nan] = averts_basic[where_nan]
return averts
[docs] def array_at_faces(self, a, direction, withnan=True):
"""
Computes values at the center of cell faces using linear interpolation.
Parameters
----------
a : ndarray
Values at cell centers, shape (nlay, row, ncol).
direction : str, possible values are 'x', 'y' and 'z'
Direction in which values will be interpolated at cell faces.
withnan : bool
If True (default), the result value will be set to NaN where the
cell face sits between inactive cells. If False, not.
Returns
-------
afaces : ndarray
Array values interpolated at cell vertices, shape as input extended
by 1 along the specified direction.
"""
# get the dimension that corresponds to the direction
dir_to_dim = {"x": 2, "y": 1, "z": 0}
dim = dir_to_dim[direction]
# extended array with ghost cells on both sides having zero values
ghost_shape = list(a.shape)
ghost_shape[dim] += 2
a_ghost = np.zeros(ghost_shape, dtype=a.dtype)
# extended delta with ghost cells on both sides having zero values
delta_ghost = np.zeros(ghost_shape, dtype=a.dtype)
# inactive bool array
if withnan and self._idomain is not None:
inactive = self._idomain == 0
if dim == 0:
# fill array with ghost cells
a_ghost[1:-1, :, :] = a
a_ghost[0, :, :] = a[0, :, :]
a_ghost[-1, :, :] = a[-1, :, :]
# calculate weights
delta_ghost[1:-1, :, :] = self.delz
weight2 = delta_ghost[:-1, :, :] / (
delta_ghost[:-1, :, :] + delta_ghost[1:, :, :]
)
weight1 = 1.0 - weight2
# interpolate
afaces = a_ghost[:-1, :, :] * weight1 + a_ghost[1:, :, :] * weight2
# assign NaN where idomain==0 on both sides
if withnan and self._idomain is not None:
inactive_faces = np.full(afaces.shape, True)
inactive_faces[:-1, :, :] = np.logical_and(
inactive_faces[:-1, :, :], inactive
)
inactive_faces[1:, :, :] = np.logical_and(
inactive_faces[1:, :, :], inactive
)
afaces[inactive_faces] = np.nan
elif dim == 1:
# fill array with ghost cells
a_ghost[:, 1:-1, :] = a
a_ghost[:, 0, :] = a[:, 0, :]
a_ghost[:, -1, :] = a[:, -1, :]
# calculate weights
delc = np.reshape(self.delc, (1, self.nrow, 1))
delc_3D = delc * np.ones(a.shape)
delta_ghost[:, 1:-1, :] = delc_3D
weight2 = delta_ghost[:, :-1, :] / (
delta_ghost[:, :-1, :] + delta_ghost[:, 1:, :]
)
weight1 = 1.0 - weight2
# interpolate
afaces = a_ghost[:, :-1, :] * weight1 + a_ghost[:, 1:, :] * weight2
# assign NaN where idomain==0 on both sides
if withnan and self._idomain is not None:
inactive_faces = np.full(afaces.shape, True)
inactive_faces[:, :-1, :] = np.logical_and(
inactive_faces[:, :-1, :], inactive
)
inactive_faces[:, 1:, :] = np.logical_and(
inactive_faces[:, 1:, :], inactive
)
afaces[inactive_faces] = np.nan
elif dim == 2:
# fill array with ghost cells
a_ghost[:, :, 1:-1] = a
a_ghost[:, :, 0] = a[:, :, 0]
a_ghost[:, :, -1] = a[:, :, -1]
# calculate weights
delr = np.reshape(self.delr, (1, 1, self.ncol))
delr_3D = delr * np.ones(a.shape)
delta_ghost[:, :, 1:-1] = delr_3D
weight2 = delta_ghost[:, :, :-1] / (
delta_ghost[:, :, :-1] + delta_ghost[:, :, 1:]
)
weight1 = 1.0 - weight2
# interpolate
afaces = a_ghost[:, :, :-1] * weight1 + a_ghost[:, :, 1:] * weight2
# assign NaN where idomain==0 on both sides
if withnan and self._idomain is not None:
inactive_faces = np.full(afaces.shape, True)
inactive_faces[:, :, :-1] = np.logical_and(
inactive_faces[:, :, :-1], inactive
)
inactive_faces[:, :, 1:] = np.logical_and(
inactive_faces[:, :, 1:], inactive
)
afaces[inactive_faces] = np.nan
return afaces
@property
def cross_section_vertices(self):
"""
Get a set of xvertices and yvertices ordered by node
for plotting cross sections
Returns
-------
xverts, yverts: (np.ndarray, np.ndarray)
"""
xv = self.xyzvertices[0]
yv = self.xyzvertices[1]
xverts, yverts = [], []
for i in range(self.nrow):
for j in range(self.ncol):
xverts.append(
[
xv[i, j],
xv[i + 1, j],
xv[i + 1, j + 1],
xv[i, j + 1],
xv[i, j],
]
)
yverts.append(
[
yv[i, j],
yv[i + 1, j],
yv[i + 1, j + 1],
yv[i, j + 1],
yv[i, j],
]
)
return np.array(xverts), np.array(yverts)
[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.nrow / self.ncol
nplottable = int(nplottable)
return nplottable
[docs] def get_plottable_layer_array(self, a, layer):
# ensure plotarray is 2d and correct shape
required_shape = self.get_plottable_layer_shape()
if a.ndim == 3:
plotarray = a[layer, :, :]
elif a.ndim == 2:
plotarray = a
elif a.ndim == 1:
plotarray = a
if plotarray.shape[0] == self.nrow * self.ncol:
plotarray = plotarray.reshape(required_shape)
elif plotarray.shape[0] == self.nnodes:
plotarray = plotarray.reshape(self.shape)
plotarray = plotarray[layer, :, :]
else:
raise ValueError("Array to plot must be of dimension 1, 2, or 3")
msg = f"{plotarray.shape} /= {required_shape}"
assert plotarray.shape == required_shape, msg
return plotarray
def _set_structured_iverts(self):
"""
Build a list of the vertices that define each model cell and the x, y
pair for each vertex
"""
rowarr = np.repeat(np.arange(self.nrow, dtype=int), self.ncol)
colarr = np.tile(np.arange(self.ncol, dtype=int), self.nrow)
iverts = np.empty((4, self.ncpl), dtype=int)
iverts[0] = rowarr * (self.ncol + 1) + colarr
iverts[1] = rowarr * (self.ncol + 1) + colarr + 1
iverts[2] = (rowarr + 1) * (self.ncol + 1) + colarr + 1
iverts[3] = (rowarr + 1) * (self.ncol + 1) + colarr
self._iverts = iverts.T.tolist()
return
def _build_structured_iverts(self, i, j):
"""
Build list of vertex numbers for a cell
Parameters
----------
i : int
row index
j : int
column index
Returns
-------
iv_list : list
list of vertex numbers for a cell
"""
iv_list = [i * (self.ncol + 1) + j]
iv_list.append(i * (self.ncol + 1) + j + 1)
iv_list.append((i + 1) * (self.ncol + 1) + j + 1)
iv_list.append((i + 1) * (self.ncol + 1) + j)
return iv_list
def _set_structured_verts(self):
"""
Build a list of the vertices that define each model cell and the x, y
pair for each vertex
Returns
-------
verts : np.ndarray
Array with x, y pairs for every vertex used to define the model.
"""
self._verts = np.column_stack(
(self.xvertices.flatten(), self.yvertices.flatten())
)
return
# Importing
# initialize grid from a grb file
[docs] @classmethod
def from_binary_grid_file(cls, file_path, verbose=False):
"""
Instantiate a StructuredGrid 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 : StructuredGrid
"""
from ..mf6.utils.binarygrid_util import MfGrdFile
grb_obj = MfGrdFile(file_path, verbose=verbose)
if grb_obj.grid_type != "DIS":
raise ValueError(
f"Binary grid file ({os.path.basename(file_path)}) "
"is not a structured (DIS) grid."
)
idomain = grb_obj.idomain
xorigin = grb_obj.xorigin
yorigin = grb_obj.yorigin
angrot = grb_obj.angrot
nlay, nrow, ncol = (grb_obj.nlay, grb_obj.nrow, grb_obj.ncol)
delr, delc = grb_obj.delr, grb_obj.delc
top, botm = grb_obj.top, grb_obj.bot
top.shape = (nrow, ncol)
botm.shape = (nlay, nrow, ncol)
return cls(
delc,
delr,
top,
botm,
idomain=idomain,
xoff=xorigin,
yoff=yorigin,
angrot=angrot,
)
[docs] @classmethod
def from_gridspec(cls, file_path: Union[str, PathLike], lenuni=0):
"""
Instantiate a StructuredGrid from grid specification file.
Parameters
----------
file_path: str or PathLike
Path to the grid specification file
lenuni: int
Length unit code
Returns
-------
A StructuredGrid
"""
with open(file_path) as f:
raw = f.readline().strip().split()
nrow = int(raw[0])
ncol = int(raw[1])
raw = f.readline().strip().split()
xul, yul, rot = float(raw[0]), float(raw[1]), float(raw[2])
delr = []
j = 0
while j < ncol:
raw = f.readline().strip().split()
for r in raw:
if "*" in r:
rraw = r.split("*")
for n in range(int(rraw[0])):
delr.append(float(rraw[1]))
j += 1
else:
delr.append(float(r))
j += 1
delc = []
i = 0
while i < nrow:
raw = f.readline().strip().split()
for r in raw:
if "*" in r:
rraw = r.split("*")
for n in range(int(rraw[0])):
delc.append(float(rraw[1]))
i += 1
else:
delc.append(float(r))
i += 1
grd = cls(np.array(delc), np.array(delr), lenuni=lenuni)
xll = grd._xul_to_xll(xul, angrot=rot)
yll = grd._yul_to_yll(yul, angrot=rot)
grd.set_coord_info(xoff=xll, yoff=yll, angrot=rot)
return grd