Skip to content

Feat(dis_properties): add grid method to get discretization properties #2700

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
jlarsen-usgs wants to merge 3 commits into
base: develop
Choose a base branch
from
Open

Feat(dis_properties): add grid method to get discretization properties #2700

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
435b3c1
Feat(dis_properties): add method to get discretization properties
jlarsen-usgs Jan 31, 2026
b99796c
woof
jlarsen-usgs Jan 31, 2026
2b1778c
woof
jlarsen-usgs Jan 31, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
95 changes: 94 additions & 1 deletion autotest/test_grid.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,13 @@
from autotest.test_dis_cases import case_dis, case_disv
from autotest.test_grid_cases import GridCases
from flopy.discretization import StructuredGrid, UnstructuredGrid, VertexGrid
from flopy.mf6 import MFSimulation
from flopy.mf6 import (
MFSimulation,
ModflowGwf,
ModflowGwfdis,
ModflowGwfdisu,
ModflowGwfdisv,
)
from flopy.modflow import Modflow, ModflowDis
from flopy.utils import import_optional_dependency
from flopy.utils.crs import get_authority_crs
Expand Down Expand Up @@ -1903,3 +1909,90 @@ def test_unstructured_grid_get_node():

with pytest.raises(IndexError, match=r"Node .* out of range"):
ug.get_node(200)


@pytest.mark.mf6
def test_structured_mf6_gridprops(example_data_path):
sim = MFSimulation.load(sim_ws=example_data_path / "mf6-freyberg")
gwf = sim.get_model()
dis = gwf.dis
modelgrid = gwf.modelgrid

new_sim = MFSimulation()
new_gwf = ModflowGwf(new_sim)
new_dis = ModflowGwfdis(new_gwf, **modelgrid.dis_properties())
attrs = ("delc", "delr", "top", "botm", "idomain", "xorigin", "yorigin", "angrot")
for attr in attrs:
v0 = getattr(dis, attr).array
v1 = getattr(new_dis, attr).array
if attr in ("xorigin", "yorigin", "angrot") and v0 is None:
v0 = 0
np.testing.assert_allclose(
v0, v1, err_msg=f"{attr} not consistent with valid array data"
)


@pytest.mark.mf6
def test_vertex_mf6_gridprops(example_data_path):
sim = MFSimulation.load(sim_ws=example_data_path / "mf6" / "test003_gwftri_disv")
gwf = sim.get_model()
disv = gwf.disv
modelgrid = gwf.modelgrid

new_sim = MFSimulation()
new_gwf = ModflowGwf(new_sim)
new_disv = ModflowGwfdisv(new_gwf, **modelgrid.disv_properties())

attrs = (
"vertices",
"top",
"botm",
"idomain",
"xorigin",
"yorigin",
"angrot",
"cell2d",
)
for attr in attrs:
v0 = getattr(disv, attr).array
v1 = getattr(new_disv, attr).array
if attr in ("xorigin", "yorigin", "angrot") and v0 is None:
v0 = 0

if attr in ("cell2d", "vertices"):
for col in v0.dtype.names:
np.testing.assert_allclose(
v0[col],
v1[col],
err_msg=f"{attr} not consistent with valid array data",
)
else:
np.testing.assert_allclose(
v0, v1, err_msg=f"{attr} not consistent with valid array data"
)


@pytest.mark.mf6
def test_unstructured_mf6_gridprops(example_data_path):
sim = MFSimulation.load(sim_ws=example_data_path / "mf6" / "test006_gwf3")
gwf = sim.get_model()
disu = gwf.disu
modelgrid = gwf.modelgrid

new_sim = MFSimulation()
new_gwf = ModflowGwf(new_sim)
dis_props = modelgrid.disu_properties()
dis_props["area"] = disu.area.array
dis_props["cl12"] = disu.cl12.array
new_disu = ModflowGwfdisu(new_gwf, **dis_props)

attrs = ("top", "bot", "iac", "ja", "nodes", "ihc", "xorigin", "yorigin", "angrot")
for attr in attrs:
v0 = getattr(disu, attr).array
v1 = getattr(new_disu, attr).array
if attr in ("xorigin", "yorigin", "angrot") and v0 is None:
v0 = 0

np.testing.assert_allclose(
v0, v1, err_msg=f"{attr} not consistent with valid array data"
)
19 changes: 19 additions & 0 deletions flopy/discretization/grid.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -535,6 +535,25 @@ def xyzextent(self):
np.max(self.xyzvertices[2]),
)

@property
def area(self):
"""
Returns a numpy array of cell areas calculated using the shoelace algorithm
"""
# irregular_shape_patch
from ..plot.plotutil import UnstructuredPlotUtilities

# when looping through to create determinants, need to start at -1
xverts, yverts = self.cross_section_vertices
xverts, yverts = UnstructuredPlotUtilities.irregular_shape_patch(xverts, yverts)
area_x2 = np.zeros((1, len(xverts)))
for i in range(xverts.shape[-1]):
# calculate the determinant of each line in polygon
area_x2 += xverts[:, i - 1] * yverts[:, i] - yverts[:, i - 1] * xverts[:, i]

area = np.abs(area_x2 / 2.0)
return np.ravel(area)

@property
def grid_lines(self):
raise NotImplementedError("must define grid_lines in child class")
Expand Down
38 changes: 38 additions & 0 deletions flopy/discretization/structuredgrid.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -759,6 +759,44 @@ def map_polygons(self):

return self._polygons

def dis_properties(self, mf2005=False):
"""
Method to get DIS package properties

Parameters
----------
mf2005 : bool
flag to get legacy mf2005/mfnwt discretization package properties
from the modelgrid object

Returns
-------
dict : dictionary of discretization properties that can be used to build a
DIS package
"""
dis_props = {
"delc": self.__delc,
"delr": self.__delr,
"top": self.top,
"botm": self.botm,
"nlay": self.nlay,
"nrow": self.nrow,
"ncol": self.ncol,
}

if mf2005:
if self.is_valid:
dis_props["xul"] = self.xvertices[0, 0]
dis_props["yul"] = self.yvertices[0, 0]
dis_props["rotation"] = self.angrot
else:
dis_props["xorigin"] = self.xoffset
dis_props["yorigin"] = self.yoffset
dis_props["angrot"] = self.angrot
dis_props["idomain"] = self.idomain

return dis_props

def to_geodataframe(self):
"""
Returns a geopandas GeoDataFrame of the model grid
Expand Down
46 changes: 46 additions & 0 deletions flopy/discretization/unstructuredgrid.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -591,6 +591,52 @@ def map_polygons(self):

return copy.copy(self._polygons)

def disu_properties(self, mfusg=False):
"""
Method that returns disu properties from the grid for constructing
DISU packages for MFUSG and MF6. Note: not all required information for
DISU construction is stored in the UnstructuredGrid class, CL12 and HWVA
is not available from the Grid. Please double-check the properties returned
on a case by case basis and fill in where necessary for individual applications

Returns
-------
dict : dictionary of unstructured discretization properties
"""

dis_props = {
"top": self._top,
"bot": self._botm,
"iac": self._iac,
"ja": self._ja,
}
if mfusg:
dis_props["nodelay"] = self.ncpl
dis_props["ivc"] = np.where(self._ihc < 1, 1, 0)

else:
dis_props["nodes"] = self.nnodes
dis_props["ihc"] = self._ihc
dis_props["idomain"] = self.idomain
dis_props["xorigin"] = self.xoffset
dis_props["yorigin"] = self.yoffset
dis_props["angrot"] = self.angrot

if self.is_valid:
dis_props["vertices"] = self._vertices
cell2d = []
for ix, iv in enumerate(self._iverts):
c2d = tuple(
[ix + 1, self._xc[ix], self._yc[ix], len(iv)] + list(iv)
)
cell2d.append(c2d)
dis_props["cell2d"] = cell2d

if self.is_valid:
dis_props["area"] = self.area

return dis_props

def to_geodataframe(self):
"""
Returns a geopandas GeoDataFrame of the model grid
Expand Down
27 changes: 27 additions & 0 deletions flopy/discretization/vertexgrid.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -300,6 +300,33 @@ def map_polygons(self):

return copy.copy(self._polygons)

def disv_properties(self):
"""
Method to get DISV package properties

Returns
-------
dict : dictionary of properties that can be used to build a DISV package
"""
dis_props = {
"nlay": self.nlay,
"ncpl": self.ncpl,
"vertices": self._vertices,
"top": self.top,
"botm": self.botm,
"idomain": self.idomain,
"xorigin": self.xoffset,
"yorigin": self.yoffset,
"angrot": self.angrot,
}

if self._cell2d is None and self._cell1d is not None:
dis_props["cell1d"] = self._cell1d
else:
dis_props["cell2d"] = self._cell2d

return dis_props

def to_geodataframe(self):
"""
Returns a geopandas GeoDataFrame of the model grid
Expand Down
Loading