feat(dis_properties): add grid method to get discretization properties

.docs/code.rst

@@ -196,7 +196,6 @@ Contents:
    ./source/flopy.mf6.utils.mfobservation.rst
    ./source/flopy.mf6.utils.output_util.rst
    ./source/flopy.mf6.utils.postprocessing.rst
-   ./source/flopy.mf6.utils.reference.rst
    ./source/flopy.mf6.utils.lakpak_utils.rst
    ./source/flopy.mf6.utils.model_splitter.rst
 

.docs/md/optional_dependencies.md

@@ -11,7 +11,6 @@ Dependencies for optional features are listed below. These may be installed with
 | `.export(*.tif)`                                                                     | **rasterio**                                                             |
 | `.export_array(*.asc)` in `flopy.export.utils`                                       | **scipy.ndimage**                                                        |
 | `.resample_to_grid()` in `flopy.utils.rasters`                                       | **scipy.interpolate**                                                    |
-| `.interpolate()` in `flopy.mf6.utils.reference` `StructuredSpatialReference` class   | **scipy.interpolate**                                                    |
 | `.get_authority_crs()` in `flopy.utils.crs`                                          | **pyproj** >= 2.2.0                                                      |
 | `.generate_classes()` in `flopy.mf6.utils`                                           | [**modflow-devtools**](https://github.com/MODFLOW-ORG/modflow-devtools)  |
 | `GridIntersect()` in `flopy.utils.gridintersect`                                     | **shapely**                                                              |

autotest/test_grid.py

@@ -17,7 +17,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
@@ -1898,3 +1904,258 @@ 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"
+        )
+
+
+def test_structured_mf2005_gridprops(example_data_path):
+    mf = Modflow.load("freyberg.nam", model_ws=example_data_path / "freyberg")
+    dis = mf.dis
+    modelgrid = mf.modelgrid
+    modelgrid.set_coord_info(0, 0, 0)
+
+    new_model = Modflow()
+    new_dis = ModflowDis(new_model, **modelgrid.dis_properties(mf2005=True))
+    attrs = ("delc", "delr", "top", "botm", "nlay", "nrow", "ncol")
+    for attr in attrs:
+        v0 = getattr(dis, attr)
+        v1 = getattr(new_dis, attr)
+        if hasattr(v0, "array"):
+            v0 = v0.array
+            v1 = v1.array
+
+        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"
+        )
+
+
+def test_unstructured_mf6_gridprops2():
+    nnodes = 2
+    top = np.ones((nnodes,))
+    botm = np.zeros((nnodes,))
+    area = np.full((nnodes,), 10)
+    idomain = np.ones((nnodes,), dtype=int)
+    iac = [2, 2]
+    ja = [1, 2, 2, 1]
+    ihc = [0, 1, 0, 1]
+    cl12 = [
+        0,
+        10,
+        0,
+        10,
+    ]
+    hwva = [
+        0,
+        100,
+        0,
+        100,
+    ]
+
+    vertices = [
+        [0, 0, 0],
+        [1, 0, 10],
+        [2, 10, 10],
+        [3, 10, 0],
+        [4, 10, 20],
+        [5, 20, 20],
+    ]
+
+    cell2d = [[0, 5, 5, 5, 0, 1, 2, 3, 0], [1, 15, 5, 5, 3, 2, 4, 5, 3]]
+    xoff = 100
+    yoff = 100
+    angrot = 10
+
+    sim = MFSimulation()
+    gwf = ModflowGwf(sim, modelname="usg_test2")
+    disu = ModflowGwfdisu(
+        gwf,
+        xorigin=xoff,
+        yorigin=yoff,
+        angrot=angrot,
+        nodes=nnodes,
+        nja=len(ja),
+        nvert=len(vertices),
+        top=top,
+        bot=botm,
+        area=area,
+        idomain=idomain,
+        iac=iac,
+        ja=ja,
+        ihc=ihc,
+        cl12=cl12,
+        hwva=hwva,
+        vertices=vertices,
+        cell2d=cell2d,
+    )
+    modelgrid = gwf.modelgrid
+
+    sim2 = MFSimulation()
+    gwf2 = ModflowGwf(sim2)
+    disu2 = ModflowGwfdisu(gwf2, cl12=cl12, hwva=hwva, **modelgrid.disu_properties())
+
+    attrs = (
+        "top",
+        "bot",
+        "iac",
+        "ja",
+        "nodes",
+        "cl12",
+        "hwva",
+        "ihc",
+        "cell2d",
+        "vertices",
+        "xorigin",
+        "yorigin",
+        "angrot",
+    )
+    for attr in attrs:
+        v0 = getattr(disu, attr).array
+        v1 = getattr(disu2, attr).array
+        if attr in ("cell2d", "vertices"):
+            for col in v0.dtype.names:
+                np.testing.assert_allclose(
+                    v0[col],
+                    v1[col],
+                    err_msg=f"{attr} column: {col} not "
+                    f"consistent with valid array data",
+                )
+        else:
+            np.testing.assert_allclose(
+                v0, v1, err_msg=f"{attr} not consistent with valid array data"
+            )
+
+
+def test_area():
+    import random
+
+    nlay = 1
+    nrow = 1
+    ncol = 1
+    dy = random.random() * 10
+    dx = random.random() * 10
+    valid_area = dx * dy
+    delc = np.full((nrow,), dy)
+    delr = np.full((ncol,), dx)
+    top = np.ones((nrow, ncol))
+    botm = np.zeros((nlay, nrow, ncol), dtype=int)
+    sgrid = StructuredGrid(delc=delc, delr=delr, nlay=1, top=top, botm=botm)
+    cell_area = sgrid.area
+    np.testing.assert_allclose(
+        [
+            valid_area,
+        ],
+        cell_area,
+        err_msg="shoelace algorithm not returning valid area within tolerance",
+    )
+
+    # triangle test
+    x1 = random.random() * 10
+    x2 = x1 / 2
+    y2 = random.random() * 10
+    verts = np.array([[0, 0, 0], [1, x1, 0], [2, x2, y2]])
+    # a = 0.5 * b * h
+    valid_area = 0.5 * x1 * y2
+
+    xc = np.mean(verts.T[1])
+    yc = np.mean(verts.T[2])
+    cell2d = [
+        (0, xc, yc, 4, 0, 1, 2, 0),
+    ]
+    nlay = 1
+    top = np.ones((len(cell2d),))
+    botm = np.zeros((nlay, len(cell2d)))
+
+    vgrid = VertexGrid(vertices=verts, cell2d=cell2d, nlay=nlay, top=top, botm=botm)
+    cell_area = vgrid.area
+    np.testing.assert_allclose(
+        [
+            valid_area,
+        ],
+        cell_area,
+        err_msg="shoelace algorithm not returning valid area within tolerance",
+    )

flopy/discretization/grid.py

@@ -411,29 +411,6 @@ def laycbd(self):
         else:
             return self._laycbd
 
-    @property
-    def cell_area(self):
-        """
-        Use shoelace algorithm for non-self-intersecting polygons to
-        calculate area.
-
-        Returns
-        -------
-            area : np.ndarray
-                numpy array of cell areas in L^2
-        """
-        from ..plot.plotutil import UnstructuredPlotUtilities
-
-        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 cell_thickness(self):
         """
@@ -558,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.0)
+        return np.ravel(area)
+
     @property
     def grid_lines(self):
         raise NotImplementedError("must define grid_lines in child class")

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