From 34cff25b47d7ced14923dd0f3d4ad2709d041446 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Wed, 5 Aug 2026 13:49:44 -0700 Subject: [PATCH 01/10] update(Raster, ModelTime): Deprecation warnings maintenance: * replace numpy .shape = with .reshape() in Raster * replace .utcfromtimestamp with fromtimestamp(timestamp, timezone.utc) in ModelTime --- flopy/discretization/modeltime.py | 5 +++-- flopy/utils/rasters.py | 2 +- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/flopy/discretization/modeltime.py b/flopy/discretization/modeltime.py index 2056c0110..8ef7a93eb 100644 --- a/flopy/discretization/modeltime.py +++ b/flopy/discretization/modeltime.py @@ -1,6 +1,6 @@ import calendar from dataclasses import dataclass, field -from datetime import datetime, timedelta +from datetime import datetime, timedelta, timezone from difflib import SequenceMatcher import numpy as np @@ -479,7 +479,8 @@ def parse_datetime( elif isinstance(datetime_obj, np.datetime64): unix_time_0 = datetime(1970, 1, 1) ts = (datetime_obj - np.datetime64(unix_time_0)) / np.timedelta64(1, "s") - datetime_obj = datetime.utcfromtimestamp(ts) + datetime_obj = datetime.fromtimestamp(ts, tz=timezone.utc) + datetime_obj = datetime_obj.replace(tzinfo=None) elif isinstance(datetime_obj, pd.Timestamp): datetime_obj = datetime_obj.to_pydatetime() elif isinstance(datetime_obj, datetime): diff --git a/flopy/utils/rasters.py b/flopy/utils/rasters.py index 10489e308..caae365f1 100644 --- a/flopy/utils/rasters.py +++ b/flopy/utils/rasters.py @@ -588,7 +588,7 @@ def resample_to_grid( data = np.where(np.isnan(data), extrapolate, data) # step 4: return grid to user in shape provided - data.shape = data_shape + data = data.reshape(data_shape) # step 5: re-apply nodata values data[np.isnan(data)] = self.nodatavals[0] From bfa371d64ca048a839f30dd2120806c8f601bc86 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Wed, 5 Aug 2026 14:02:33 -0700 Subject: [PATCH 02/10] maintenance `.shape =` calls replaced with `arr = arr.reshape()` * numpy will be deprecating `.shape =` reshaping, migrating to reshape convention to prevent DeprecationWarnings --- flopy/discretization/structuredgrid.py | 4 ++-- flopy/discretization/vertexgrid.py | 2 +- flopy/mf6/utils/binaryfile_utils.py | 4 ++-- flopy/mf6/utils/binarygrid_util.py | 6 +++--- flopy/mf6/utils/model_splitter.py | 2 +- flopy/plot/plotutil.py | 10 +++++----- flopy/utils/postprocessing.py | 10 +++++----- flopy/utils/zonbud.py | 2 +- 8 files changed, 20 insertions(+), 20 deletions(-) diff --git a/flopy/discretization/structuredgrid.py b/flopy/discretization/structuredgrid.py index 18bab69f8..7aa6bcc40 100644 --- a/flopy/discretization/structuredgrid.py +++ b/flopy/discretization/structuredgrid.py @@ -2057,8 +2057,8 @@ def from_binary_grid_file(cls, file_path, verbose=False): 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) + top = top.reshape((nrow, ncol)) + botm = botm.reshape((nlay, nrow, ncol)) return cls( delc, delr, diff --git a/flopy/discretization/vertexgrid.py b/flopy/discretization/vertexgrid.py index ddb8da4b5..66f23dd57 100644 --- a/flopy/discretization/vertexgrid.py +++ b/flopy/discretization/vertexgrid.py @@ -847,7 +847,7 @@ def from_binary_grid_file(cls, file_path, verbose=False): nlay, ncpl = grb_obj.nlay, grb_obj.ncpl top = np.ravel(grb_obj.top) botm = grb_obj.bot - botm.shape = (nlay, ncpl) + botm = botm.reshape((nlay, ncpl)) vertices, cell2d = grb_obj.cell2d return cls( diff --git a/flopy/mf6/utils/binaryfile_utils.py b/flopy/mf6/utils/binaryfile_utils.py index 91bcd0abd..941ba31ab 100644 --- a/flopy/mf6/utils/binaryfile_utils.py +++ b/flopy/mf6/utils/binaryfile_utils.py @@ -375,9 +375,9 @@ def _reshape_binary_data(data, dtype=None): return data elif dtype == "V": nodes = len(data[0][0][0]) - data.shape = (time, -1, nodes) + data = data.reshape((time, -1, nodes)) elif dtype == "U": - data.shape = (time, -1) + data = data.reshape((time, -1)) else: err = "Invalid dtype flag supplied, valid are dtype='U', dtype='V'" raise Exception(err) diff --git a/flopy/mf6/utils/binarygrid_util.py b/flopy/mf6/utils/binarygrid_util.py index 88e9b565c..c74cc20f9 100644 --- a/flopy/mf6/utils/binarygrid_util.py +++ b/flopy/mf6/utils/binarygrid_util.py @@ -238,7 +238,7 @@ def _set_modelgrid(self): nlay, ncpl = self.nlay, self.ncpl vertices, cell2d = self.cell2d top = np.ravel(top) - botm.shape = (nlay, ncpl) + botm = botm.reshape((nlay, ncpl)) modelgrid = VertexGrid( vertices, cell2d, @@ -258,8 +258,8 @@ def _set_modelgrid(self): ) delr, delc = self.delr, self.delc - top.shape = (nrow, ncol) - botm.shape = (nlay, nrow, ncol) + top.reshape((nrow, ncol)) + botm = botm.reshape((nlay, nrow, ncol)) modelgrid = StructuredGrid( delc, delr, diff --git a/flopy/mf6/utils/model_splitter.py b/flopy/mf6/utils/model_splitter.py index 72bd0f3a9..8b333cb5e 100644 --- a/flopy/mf6/utils/model_splitter.py +++ b/flopy/mf6/utils/model_splitter.py @@ -833,7 +833,7 @@ def reconstruct_array(self, arrays): new_array[new_nodes] = array[old_nodes] - new_array.shape = shape + new_array = new_array.reshape(shape) return new_array def reconstruct_recarray(self, recarrays): diff --git a/flopy/plot/plotutil.py b/flopy/plot/plotutil.py index 23f6d00f3..fdce65021 100644 --- a/flopy/plot/plotutil.py +++ b/flopy/plot/plotutil.py @@ -1482,11 +1482,11 @@ def saturated_thickness(head, top, botm, laytyp, mask_values=None): head = np.copy(head) nlay, nrow, ncol = head.shape ncpl = nrow * ncol - head.shape = (nlay, ncpl) - top.shape = (ncpl,) - botm.shape = (nlay, ncpl) + head = head.reshape((nlay, ncpl)) + top = top.reshape((ncpl,)) + botm = botm.reshape((nlay, ncpl)) if laytyp.ndim == 3: - laytyp.shape = (nlay, ncpl) + laytyp = laytyp.reshape((nlay, ncpl)) else: nrow, ncol = None, None @@ -1531,7 +1531,7 @@ def saturated_thickness(head, top, botm, laytyp, mask_values=None): sat_thk = np.where(laytyp != 0, sat_thk_unconf, sat_thk_conf) if nrow is not None and ncol is not None: - sat_thk.shape = (nlay, nrow, ncol) + sat_thk = sat_thk.reshape((nlay, nrow, ncol)) return sat_thk diff --git a/flopy/utils/postprocessing.py b/flopy/utils/postprocessing.py index ec9416626..656f7d0f4 100644 --- a/flopy/utils/postprocessing.py +++ b/flopy/utils/postprocessing.py @@ -801,7 +801,7 @@ def get_specific_discharge( modelgrid = model.modelgrid if head is not None: - head.shape = modelgrid.shape + head = head.reshape(modelgrid.shape) if isinstance(vectors, (list, tuple)): classical_budget = True @@ -857,7 +857,7 @@ def get_specific_discharge( head, mask=[model.hdry, model.hnoflo] ) - saturated_thickness.shape = modelgrid.shape + saturated_thickness = saturated_thickness.reshape(modelgrid.shape) # inform modelgrid of no-flow and dry cells modelgrid = model.modelgrid @@ -927,9 +927,9 @@ def get_specific_discharge( qx[idx] = spdis["qx"] qy[idx] = spdis["qy"] qz[idx] = spdis["qz"] - qx.shape = modelgrid.shape - qy.shape = modelgrid.shape - qz.shape = modelgrid.shape + qx = qx.reshape(modelgrid.shape) + qy = qy.reshape(modelgrid.shape) + qz = qz.reshape(modelgrid.shape) # set no-flow and dry cells to NaN if head is not None and position == "centers": diff --git a/flopy/utils/zonbud.py b/flopy/utils/zonbud.py index b8dd2a8ae..ee87abb93 100644 --- a/flopy/utils/zonbud.py +++ b/flopy/utils/zonbud.py @@ -2677,7 +2677,7 @@ def _read_zb_csv2(fname, add_prefix=True, aliases=None): array = np.genfromtxt(foo, delimiter=",").T if len(array) != len(dtype): array = array[:-1] - array.shape = (len(dtype), -1) + array = array.reshape((len(dtype), -1)) data = {name[0]: list(array[ix]) for ix, name in enumerate(dtype)} data["KPER"] = list(np.array(data["KPER"]) - 1) data["KSTP"] = list(np.array(data["KSTP"]) - 1) From 1e8a843c33bcb4f6447c3fc02080588213092865 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Wed, 5 Aug 2026 15:14:54 -0700 Subject: [PATCH 03/10] feat(dis_properties): add grid method to get discretization properties * returns dictionary of keyword arguments to build DIS, DISV, and DISU depending on grid type * added cell `.area` calculation via shoelace algorithm to `Grid` * remove deprecated flopy.mf6.utils/reference.py which housed "pre-modelgrid" spatial reference support for MF6 models --- .docs/md/optional_dependencies.md | 1 - autotest/test_grid.py | 95 ++- etc/environment.yml | 2 - flopy/discretization/grid.py | 20 + flopy/discretization/structuredgrid.py | 38 + flopy/discretization/unstructuredgrid.py | 51 ++ flopy/discretization/vertexgrid.py | 27 + flopy/mf6/utils/reference.py | 950 ----------------------- 8 files changed, 230 insertions(+), 954 deletions(-) delete mode 100644 flopy/mf6/utils/reference.py diff --git a/.docs/md/optional_dependencies.md b/.docs/md/optional_dependencies.md index 2ce5319c2..6db9b2cb5 100644 --- a/.docs/md/optional_dependencies.md +++ b/.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** | diff --git a/autotest/test_grid.py b/autotest/test_grid.py index 9ff8a1608..35d1cb37c 100644 --- a/autotest/test_grid.py +++ b/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,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" + ) diff --git a/etc/environment.yml b/etc/environment.yml index a3228f805..b5e8038b2 100644 --- a/etc/environment.yml +++ b/etc/environment.yml @@ -55,12 +55,10 @@ dependencies: - pymetis - pyproj - pyshp - - pyvista - rasterio - rasterstats - scipy - shapely>=2.0 - - vtk - xmipy - h5py - scikit-learn diff --git a/flopy/discretization/grid.py b/flopy/discretization/grid.py index 0a50597c9..fb26fd40f 100644 --- a/flopy/discretization/grid.py +++ b/flopy/discretization/grid.py @@ -558,6 +558,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") diff --git a/flopy/discretization/structuredgrid.py b/flopy/discretization/structuredgrid.py index 7aa6bcc40..690cf306e 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..62c6200d1 100644 --- a/flopy/discretization/unstructuredgrid.py +++ b/flopy/discretization/unstructuredgrid.py @@ -591,6 +591,57 @@ 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 + + Parameters + ---------- + mfusg : bool + boolean flag for specifying modflow USG DISU properties + + 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 66f23dd57..7794e8461 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 diff --git a/flopy/mf6/utils/reference.py b/flopy/mf6/utils/reference.py deleted file mode 100644 index 3517f12df..000000000 --- a/flopy/mf6/utils/reference.py +++ /dev/null @@ -1,950 +0,0 @@ -""" -Module spatial referencing for flopy model objects - -.. deprecated:: 3.9 - This module will be removed in FloPy 3.10+. Use - the :mod:`flopy.discretization` module instead. - -""" - -import numpy as np - - -class StructuredSpatialReference: - """ - a simple class to locate the model grid in x-y space - - .. deprecated:: 3.9 - This class will be removed in FloPy 3.10+. Use - :class:`~flopy.discretization.structuredgrid.StructuredGrid` - instead. - - Parameters - ---------- - - delr : numpy ndarray - the model discretization delr vector - - delc : numpy ndarray - the model discretization delc vector - - lenuni : int - the length units flag from the discretization package - - xul : float - the x coordinate of the upper left corner of the grid - - yul : float - the y coordinate of the upper left corner of the grid - - rotation : float - the counter-clockwise rotation (in degrees) of the grid - - proj4_str: str - a PROJ4 string that identifies the grid in space. warning: case - sensitive! - - Attributes - ---------- - xedge : ndarray - array of column edges - - yedge : ndarray - array of row edges - - xgrid : ndarray - numpy meshgrid of xedges - - ygrid : ndarray - numpy meshgrid of yedges - - xcenter : ndarray - array of column centers - - ycenter : ndarray - array of row centers - - xcentergrid : ndarray - numpy meshgrid of column centers - - ycentergrid : ndarray - numpy meshgrid of row centers - - Notes - ----- - - xul and yul can be explicitly (re)set after SpatialReference - instantiation, but only before any of the other attributes and methods are - accessed - - """ - - def __init__( - self, - delr=1.0, - delc=1.0, - lenuni=1, - nlay=1, - xul=None, - yul=None, - rotation=0.0, - proj4_str=None, - **kwargs, - ): - self.delc = np.atleast_1d(np.array(delc)) - self.delr = np.atleast_1d(np.array(delr)) - self.nlay = nlay - self.lenuni = lenuni - self.proj4_str = proj4_str - self._reset() - self.set_spatialreference(xul, yul, rotation) - - @classmethod - def from_namfile_header(cls, namefile): - # check for reference info in the nam file header - header = [] - with open(namefile) as f: - for line in f: - if not line.startswith("#"): - break - header.extend(line.strip().replace("#", "").split(",")) - - xul, yul = None, None - rotation = 0.0 - proj4_str = None - start_datetime = "1/1/1970" - - for item in header: - if "xul" in item.lower(): - try: - xul = float(item.split(":")[1]) - except: - pass - elif "yul" in item.lower(): - try: - yul = float(item.split(":")[1]) - except: - pass - elif "rotation" in item.lower(): - try: - rotation = float(item.split(":")[1]) - except: - pass - elif "proj4_str" in item.lower(): - try: - proj4_str = ":".join(item.split(":")[1:]).strip() - except: - pass - elif "start" in item.lower(): - try: - start_datetime = item.split(":")[1].strip() - except: - pass - - return ( - cls(xul=xul, yul=yul, rotation=rotation, proj4_str=proj4_str), - start_datetime, - ) - - def __setattr__(self, key, value): - reset = True - if key == "delr": - super().__setattr__("delr", np.atleast_1d(np.array(value))) - elif key == "delc": - super().__setattr__("delc", np.atleast_1d(np.array(value))) - elif key == "xul": - super().__setattr__("xul", float(value)) - elif key == "yul": - super().__setattr__("yul", float(value)) - elif key == "rotation": - super().__setattr__("rotation", float(value)) - elif key == "lenuni": - super().__setattr__("lenuni", int(value)) - elif key == "nlay": - super().__setattr__("nlay", int(value)) - else: - super().__setattr__(key, value) - reset = False - if reset: - self._reset() - - def reset(self, **kwargs): - for key, value in kwargs.items(): - setattr(self, key, value) - - def _reset(self): - self._xgrid = None - self._ygrid = None - self._ycentergrid = None - self._xcentergrid = None - - @property - def nrow(self): - return self.delc.shape[0] - - @property - def ncol(self): - return self.delr.shape[0] - - def __eq__(self, other): - if not isinstance(other, StructuredSpatialReference): - return False - if other.xul != self.xul: - return False - if other.yul != self.yul: - return False - if other.rotation != self.rotation: - return False - if other.proj4_str != self.proj4_str: - return False - return True - - @classmethod - def from_gridspec(cls, gridspec_file, lenuni=0): - f = open(gridspec_file, "r") - 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 - f.close() - return cls( - np.array(delr), - np.array(delc), - lenuni, - xul=xul, - yul=yul, - rotation=rot, - ) - - @property - def attribute_dict(self): - return { - "xul": self.xul, - "yul": self.yul, - "rotation": self.rotation, - "proj4_str": self.proj4_str, - } - - def set_spatialreference(self, xul=None, yul=None, rotation=0.0): - """ - set spatial reference - can be called from model instance - """ - - # Set origin and rotation - if xul is None: - self.xul = 0.0 - else: - self.xul = xul - if yul is None: - self.yul = np.add.reduce(self.delc) - else: - self.yul = yul - self.rotation = rotation - self._reset() - - def __repr__(self): - s = f"xul:{self.xul: Date: Wed, 5 Aug 2026 16:07:52 -0700 Subject: [PATCH 04/10] revert environment.yml file --- etc/environment.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/etc/environment.yml b/etc/environment.yml index b5e8038b2..a3228f805 100644 --- a/etc/environment.yml +++ b/etc/environment.yml @@ -55,10 +55,12 @@ dependencies: - pymetis - pyproj - pyshp + - pyvista - rasterio - rasterstats - scipy - shapely>=2.0 + - vtk - xmipy - h5py - scikit-learn From c8afa7b79d76cfa27e6358b09c97e216409c1895 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Tue, 11 Aug 2026 09:38:11 -0700 Subject: [PATCH 05/10] add testing for shoelace algorithm (Grid.area) --- autotest/test_grid.py | 51 ++++++++++++++++++++++++++++++++++++ flopy/discretization/grid.py | 23 ---------------- flopy/utils/voronoi.py | 2 +- 3 files changed, 52 insertions(+), 24 deletions(-) diff --git a/autotest/test_grid.py b/autotest/test_grid.py index 35d1cb37c..4b4bcbf94 100644 --- a/autotest/test_grid.py +++ b/autotest/test_grid.py @@ -1991,3 +1991,54 @@ def test_unstructured_mf6_gridprops(example_data_path): 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", + ) diff --git a/flopy/discretization/grid.py b/flopy/discretization/grid.py index fb26fd40f..a51845e83 100644 --- a/flopy/discretization/grid.py +++ b/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): """ diff --git a/flopy/utils/voronoi.py b/flopy/utils/voronoi.py index e2ef04f49..78d676c0f 100644 --- a/flopy/utils/voronoi.py +++ b/flopy/utils/voronoi.py @@ -385,7 +385,7 @@ def get_disu6_gridprops(self): gridprops["cl12"] = cl12 gridprops["hwva"] = hwva gridprops["angldegx"] = angldegx - gridprops["area"] = ugrid.cell_area + gridprops["area"] = ugrid.area gridprops["nodes"] = len(iac) gridprops["nja"] = len(ja) gridprops["nvert"] = len(gridprops["vertices"]) From 80ea66fbe98e1ed99ec59a5ed8d2a81aece0aace Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Tue, 11 Aug 2026 09:41:02 -0700 Subject: [PATCH 06/10] remove reference.rst from .docs/code.rst --- .docs/code.rst | 1 - 1 file changed, 1 deletion(-) diff --git a/.docs/code.rst b/.docs/code.rst index b9a408ff0..ce360eeb1 100644 --- a/.docs/code.rst +++ b/.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 From 55442fff2a951fe4e07c94072970d56b508d47a1 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Tue, 11 Aug 2026 10:22:49 -0700 Subject: [PATCH 07/10] add test for generating mf2005 DIS --- autotest/test_grid.py | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/autotest/test_grid.py b/autotest/test_grid.py index 4b4bcbf94..d5ebae15e 100644 --- a/autotest/test_grid.py +++ b/autotest/test_grid.py @@ -1927,6 +1927,27 @@ def test_structured_mf6_gridprops(example_data_path): ) +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") From 0d8d3a5cd7174910857906d8e4250abde7251547 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Tue, 11 Aug 2026 12:39:31 -0700 Subject: [PATCH 08/10] add second unstructured grid test for "complete" model grid --- autotest/test_grid.py | 95 ++++++++++++++++++++++++ flopy/discretization/unstructuredgrid.py | 5 +- 2 files changed, 98 insertions(+), 2 deletions(-) diff --git a/autotest/test_grid.py b/autotest/test_grid.py index d5ebae15e..94427f14d 100644 --- a/autotest/test_grid.py +++ b/autotest/test_grid.py @@ -2014,6 +2014,101 @@ def test_unstructured_mf6_gridprops(example_data_path): ) +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 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 diff --git a/flopy/discretization/unstructuredgrid.py b/flopy/discretization/unstructuredgrid.py index 62c6200d1..20e67a6a8 100644 --- a/flopy/discretization/unstructuredgrid.py +++ b/flopy/discretization/unstructuredgrid.py @@ -628,11 +628,12 @@ def disu_properties(self, mfusg=False): dis_props["angrot"] = self.angrot if self.is_valid: - dis_props["vertices"] = self._vertices + dis_props["vertices"] = [[int(i) for i in v] for v in self._vertices] cell2d = [] for ix, iv in enumerate(self._iverts): c2d = tuple( - [ix + 1, self._xc[ix], self._yc[ix], len(iv)] + list(iv) + [ix, self._xc[ix], self._yc[ix], len(iv)] + + [int(i) for i in list(iv)] ) cell2d.append(c2d) dis_props["cell2d"] = cell2d From 35162755b8bd00332418b0658e1c46df574de4ca Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Tue, 11 Aug 2026 12:44:31 -0700 Subject: [PATCH 09/10] woof --- autotest/test_grid.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/autotest/test_grid.py b/autotest/test_grid.py index 94427f14d..829704411 100644 --- a/autotest/test_grid.py +++ b/autotest/test_grid.py @@ -2101,7 +2101,8 @@ def test_unstructured_mf6_gridprops2(): np.testing.assert_allclose( v0[col], v1[col], - err_msg=f"{attr} column: {col} not consistent with valid array data", + err_msg=f"{attr} column: {col} not " + f"consistent with valid array data", ) else: np.testing.assert_allclose( From 32a560971e7a5c2f81769bbc5e44d00dbe751fa0 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Tue, 11 Aug 2026 12:57:03 -0700 Subject: [PATCH 10/10] remove references to `reference.py` --- tach.toml | 8 -------- 1 file changed, 8 deletions(-) diff --git a/tach.toml b/tach.toml index 3fbd2aa00..39660cd8d 100644 --- a/tach.toml +++ b/tach.toml @@ -395,10 +395,6 @@ depends_on = [ "flopy.mf6.utils.binarygrid_util", ] -[[modules]] -path = "flopy.mf6.utils.reference" -depends_on = [] - [[modules]] path = "flopy.mf6.utils.testutils" depends_on = [ @@ -616,10 +612,6 @@ depends_on = [ path = "flopy.utils.recarray_utils" depends_on = [] -[[modules]] -path = "flopy.utils.reference" -depends_on = [] - [[modules]] path = "flopy.utils.sfroutputfile" depends_on = []