From 435b3c1aab4402274e69ee33deca09f081b8b917 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Fri, 30 Jan 2026 16:39:35 -0800 Subject: [PATCH 1/3] Feat(dis_properties): add method to get discretization properties * returns dictionary of keyword arguments to build DIS, DISV, or DISU depending on grid * added cell `.area` calculation via shoelace algorithm to grids --- autotest/test_grid.py | 130 ++++++++++++++--------- flopy/discretization/grid.py | 20 ++++ flopy/discretization/structuredgrid.py | 38 +++++++ flopy/discretization/unstructuredgrid.py | 46 ++++++++ flopy/discretization/vertexgrid.py | 27 +++++ 5 files changed, 209 insertions(+), 52 deletions(-) diff --git a/autotest/test_grid.py b/autotest/test_grid.py index 8b6b53a5b..81c92a2da 100644 --- a/autotest/test_grid.py +++ b/autotest/test_grid.py @@ -16,7 +16,7 @@ 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, ModflowGwfdisv, ModflowGwfdisu from flopy.modflow import Modflow, ModflowDis from flopy.utils import import_optional_dependency from flopy.utils.crs import get_authority_crs @@ -1232,57 +1232,6 @@ def test_tocvfd4(): assert iverts[4] == [2, 5, 13, 10, 17, 8, 2] -@requires_pkg("shapely", "scipy") -@requires_exe("triangle") -@pytest.mark.parametrize( - "max_area,domain_size", - [ - (0.1, 2), # Simple: small domain, large max_area -> few cells - (0.02, 5), # More complex: larger domain, smaller max_area -> more cells - ], - ids=["simple_voronoi", "complex_voronoi"], -) -def test_voronoi_hanging_node_check(function_tmpdir, max_area, domain_size): - """ - Addresses github.com/modflowpy/flopy/issues/2427 - - Tests hanging node check works on programmatically generated Voronoi - grids. Grids with clean geometry should successfully converge, while - the presence of floating-point artifacts like duplicate vertices may - prevent convergence. - """ - import warnings - - tri = Triangle(maximum_area=max_area, angle=30, model_ws=function_tmpdir) - poly = np.array( - ((0, 0), (domain_size, 0), (domain_size, domain_size), (0, domain_size)) - ) - tri.add_polygon(poly) - tri.build(verbose=False) - vor = VoronoiGrid(tri) - - vertdict = {} - for icell, iverts_cell in enumerate(vor.iverts): - points = [] - for iv in iverts_cell: - points.append(tuple(vor.verts[iv])) - points.append(points[0]) # close the polygon - vertdict[icell] = points - - with warnings.catch_warnings(record=True) as w: - warnings.simplefilter("always") - verts_skip, iverts_skip = to_cvfd(vertdict, skip_hanging_node_check=True) - assert len(w) == 0 - - with warnings.catch_warnings(record=True) as w: - warnings.simplefilter("always") - verts_check, iverts_check = to_cvfd(vertdict, skip_hanging_node_check=False) - assert len(w) == 0 - - assert len(iverts_skip) == len(iverts_check) - assert len(iverts_skip) == len(vertdict) - - @requires_pkg("shapely") def test_area_centroid_polygon(): pts = [ @@ -1903,3 +1852,80 @@ 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" + ) diff --git a/flopy/discretization/grid.py b/flopy/discretization/grid.py index 4c0b74e50..d68097d7b 100644 --- a/flopy/discretization/grid.py +++ b/flopy/discretization/grid.py @@ -535,6 +535,26 @@ 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.) + return np.ravel(area) + @property def grid_lines(self): raise NotImplementedError("must define grid_lines in child class") diff --git a/flopy/discretization/structuredgrid.py b/flopy/discretization/structuredgrid.py index 18bab69f8..f4efd000b 100644 --- a/flopy/discretization/structuredgrid.py +++ b/flopy/discretization/structuredgrid.py @@ -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 diff --git a/flopy/discretization/unstructuredgrid.py b/flopy/discretization/unstructuredgrid.py index 2d81e7b4a..48a8cc70f 100644 --- a/flopy/discretization/unstructuredgrid.py +++ b/flopy/discretization/unstructuredgrid.py @@ -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 diff --git a/flopy/discretization/vertexgrid.py b/flopy/discretization/vertexgrid.py index ddb8da4b5..757917d21 100644 --- a/flopy/discretization/vertexgrid.py +++ b/flopy/discretization/vertexgrid.py @@ -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 From b99796c53ba50ec38a5c959f9d1f836030aa091c Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Fri, 30 Jan 2026 16:42:40 -0800 Subject: [PATCH 2/3] woof --- autotest/test_grid.py | 77 +++++++++++++++++++++++++++++++++--- flopy/discretization/grid.py | 7 ++-- 2 files changed, 75 insertions(+), 9 deletions(-) diff --git a/autotest/test_grid.py b/autotest/test_grid.py index 81c92a2da..45d5e01c1 100644 --- a/autotest/test_grid.py +++ b/autotest/test_grid.py @@ -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, ModflowGwf, ModflowGwfdis, ModflowGwfdisv, ModflowGwfdisu +from flopy.mf6 import ( + MFSimulation, + ModflowGwf, + ModflowGwfdis, + ModflowGwfdisv, + ModflowGwfdisu, +) from flopy.modflow import Modflow, ModflowDis from flopy.utils import import_optional_dependency from flopy.utils.crs import get_authority_crs @@ -1232,6 +1238,57 @@ def test_tocvfd4(): assert iverts[4] == [2, 5, 13, 10, 17, 8, 2] +@requires_pkg("shapely", "scipy") +@requires_exe("triangle") +@pytest.mark.parametrize( + "max_area,domain_size", + [ + (0.1, 2), # Simple: small domain, large max_area -> few cells + (0.02, 5), # More complex: larger domain, smaller max_area -> more cells + ], + ids=["simple_voronoi", "complex_voronoi"], +) +def test_voronoi_hanging_node_check(function_tmpdir, max_area, domain_size): + """ + Addresses github.com/modflowpy/flopy/issues/2427 + + Tests hanging node check works on programmatically generated Voronoi + grids. Grids with clean geometry should successfully converge, while + the presence of floating-point artifacts like duplicate vertices may + prevent convergence. + """ + import warnings + + tri = Triangle(maximum_area=max_area, angle=30, model_ws=function_tmpdir) + poly = np.array( + ((0, 0), (domain_size, 0), (domain_size, domain_size), (0, domain_size)) + ) + tri.add_polygon(poly) + tri.build(verbose=False) + vor = VoronoiGrid(tri) + + vertdict = {} + for icell, iverts_cell in enumerate(vor.iverts): + points = [] + for iv in iverts_cell: + points.append(tuple(vor.verts[iv])) + points.append(points[0]) # close the polygon + vertdict[icell] = points + + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter("always") + verts_skip, iverts_skip = to_cvfd(vertdict, skip_hanging_node_check=True) + assert len(w) == 0 + + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter("always") + verts_check, iverts_check = to_cvfd(vertdict, skip_hanging_node_check=False) + assert len(w) == 0 + + assert len(iverts_skip) == len(iverts_check) + assert len(iverts_skip) == len(vertdict) + + @requires_pkg("shapely") def test_area_centroid_polygon(): pts = [ @@ -1864,7 +1921,7 @@ def test_structured_mf6_gridprops(example_data_path): 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") + attrs = ("delc", "delr", "top", "botm", "idomain", "xorigin", "yorigin", "angrot") for attr in attrs: v0 = getattr(dis, attr).array v1 = getattr(new_dis, attr).array @@ -1874,6 +1931,7 @@ def test_structured_mf6_gridprops(example_data_path): 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") @@ -1886,7 +1944,14 @@ def test_vertex_mf6_gridprops(example_data_path): new_disv = ModflowGwfdisv(new_gwf, **modelgrid.disv_properties()) attrs = ( - "vertices", "top", "botm", "idomain", "xorigin", "yorigin", "angrot", "cell2d" + "vertices", + "top", + "botm", + "idomain", + "xorigin", + "yorigin", + "angrot", + "cell2d", ) for attr in attrs: v0 = getattr(disv, attr).array @@ -1894,10 +1959,12 @@ def test_vertex_mf6_gridprops(example_data_path): if attr in ("xorigin", "yorigin", "angrot") and v0 is None: v0 = 0 - if attr in("cell2d", "vertices"): + 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" + v0[col], + v1[col], + err_msg=f"{attr} not consistent with valid array data", ) else: np.testing.assert_allclose( diff --git a/flopy/discretization/grid.py b/flopy/discretization/grid.py index d68097d7b..fd74ade4e 100644 --- a/flopy/discretization/grid.py +++ b/flopy/discretization/grid.py @@ -542,17 +542,16 @@ def area(self): """ # 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 - ) + 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.) + area = np.abs(area_x2 / 2.0) return np.ravel(area) @property From 2b1778c7930a1c8b63e33b9558ca6f4d56d375b1 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Fri, 30 Jan 2026 16:44:42 -0800 Subject: [PATCH 3/3] woof --- autotest/test_grid.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/autotest/test_grid.py b/autotest/test_grid.py index 45d5e01c1..3d6874a27 100644 --- a/autotest/test_grid.py +++ b/autotest/test_grid.py @@ -20,8 +20,8 @@ MFSimulation, ModflowGwf, ModflowGwfdis, - ModflowGwfdisv, ModflowGwfdisu, + ModflowGwfdisv, ) from flopy.modflow import Modflow, ModflowDis from flopy.utils import import_optional_dependency