From 5a5a18372827e30ba8af735b54c6353cd2ef459e Mon Sep 17 00:00:00 2001 From: Joshua Larsen Date: Tue, 11 Jan 2022 19:26:32 -0800 Subject: [PATCH 1/3] update(Grid.intersect): added optional z-coordinate to intersect() method * implemented intersect method for UnstructuredGrid --- autotest/t062_test_intersect.py | 121 ++++++++++++++++++++++- flopy/discretization/structuredgrid.py | 28 +++++- flopy/discretization/unstructuredgrid.py | 75 +++++++++++++- flopy/discretization/vertexgrid.py | 23 ++++- 4 files changed, 237 insertions(+), 10 deletions(-) diff --git a/autotest/t062_test_intersect.py b/autotest/t062_test_intersect.py index 41d43feef6..66df9d022d 100644 --- a/autotest/t062_test_intersect.py +++ b/autotest/t062_test_intersect.py @@ -1,3 +1,4 @@ +import os import flopy import numpy as np import matplotlib.pyplot as plt @@ -94,6 +95,27 @@ def get_vlist(i, j, nrow, ncol): return gwf +# simple functions to load vertices and index lists for unstructured grid +def load_verts(fname): + verts = np.genfromtxt(fname, dtype=[int, float, float], + names=['iv', 'x', 'y']) + verts['iv'] -= 1 # zero based + return verts + + +def load_iverts(fname): + f = open(fname, 'r') + iverts = [] + xc = [] + yc = [] + for line in f: + ll = line.strip().split() + iverts.append([int(i) - 1 for i in ll[4:]]) + xc.append(float(ll[1])) + yc.append(float(ll[2])) + return iverts, np.array(xc), np.array(yc) + + def test_intersection(): ml_dis = dis_model() ml_disv = disv_model() @@ -140,9 +162,11 @@ def test_intersection(): else: print("In real_world coordinates:") try: - row, col = ml_dis.modelgrid.intersect(x, y, local, forgive=forgive) + row, col = ml_dis.modelgrid.intersect( + x, y, local=local, forgive=forgive + ) cell2d_disv = ml_disv.modelgrid.intersect( - x, y, local, forgive=forgive + x, y, local=local, forgive=forgive ) except Exception as e: if not forgive and any( @@ -162,5 +186,98 @@ def test_intersection(): assert all(np.isnan([row, col, cell2d_disv])) +def test_structured_xyz_intersect(): + model_ws = os.path.join( + "..", "examples", "data", "freyberg_multilayer_transient" + ) + fname = "freyberg.nam" + + ml = flopy.modflow.Modflow.load(fname, model_ws=model_ws) + mg = ml.modelgrid + top_botm = ml.modelgrid.top_botm + xc, yc, zc = mg.xyzcellcenters + + for _ in range(10): + k = np.random.randint(0, mg.nlay, 1)[0] + i = np.random.randint(0, mg.nrow, 1)[0] + j = np.random.randint(0, mg.ncol, 1)[0] + x = xc[i, j] + y = yc[i, j] + z = zc[k, i, j] + k2, i2, j2 = ml.modelgrid.intersect(x, y, z) + if (k, i, j) != (k2, i2, j2): + raise AssertionError("Structured grid intersection failed") + + +def test_vertex_xyz_intersect(): + sim_ws = os.path.join("..", "examples", "data", "mf6", "test003_gwfs_disv") + + sim = flopy.mf6.MFSimulation.load(sim_ws=sim_ws) + ml = sim.get_model(list(sim.model_names)[0]) + mg = ml.modelgrid + + xc, yc, zc = mg.xyzcellcenters + for _ in range(10): + icell = np.random.randint(0, mg.ncpl, 1)[0] + lay = np.random.randint(0, mg.nlay, 1)[0] + x = xc[icell] + y = yc[icell] + z = zc[lay, icell] + lay1, icell1 = mg.intersect(x, y, z) + + if (lay, icell) != (lay1, icell1): + raise AssertionError("Vertex grid intersection failed") + + +def test_unstructured_xyz_intersect(): + ws = os.path.join("..", "examples", "data", "unstructured") + # usg example + name = os.path.join(ws, "ugrid_verts.dat") + verts = load_verts(name) + + name = os.path.join(ws, "ugrid_iverts.dat") + iverts, xc, yc = load_iverts(name) + + # create a 3 layer model grid + ncpl = np.array(3 * [len(iverts)]) + nnodes = np.sum(ncpl) + + top = np.ones((nnodes), ) + botm = np.ones((nnodes), ) + + # set top and botm elevations + i0 = 0 + i1 = ncpl[0] + elevs = [100, 0, -100, -200] + for ix, cpl in enumerate(ncpl): + top[i0:i1] *= elevs[ix] + botm[i0:i1] *= elevs[ix + 1] + i0 += cpl + i1 += cpl + + # create the modelgrid + mg = flopy.discretization.UnstructuredGrid(vertices=verts, + iverts=iverts, + xcenters=xc, + ycenters=yc, top=top, + botm=botm, ncpl=ncpl) + + xc, yc, zc = mg.xyzcellcenters + zc = zc[0].reshape(mg.nlay, mg.ncpl[0]) + for _ in range(10): + icell = np.random.randint(0, mg.ncpl[0], 1)[0] + lay = np.random.randint(0, mg.nlay, 1)[0] + x = xc[icell] + y = yc[icell] + z = zc[lay, icell] + icell1 = mg.intersect(x, y, z) + icell = icell + (mg.ncpl[0] * lay) + if icell != icell1: + raise AssertionError("Unstructured grid intersection failed") + + if __name__ == "__main__": test_intersection() + test_structured_xyz_intersect() + test_vertex_xyz_intersect() + test_unstructured_xyz_intersect() diff --git a/flopy/discretization/structuredgrid.py b/flopy/discretization/structuredgrid.py index 0fe6309614..fca33a609d 100644 --- a/flopy/discretization/structuredgrid.py +++ b/flopy/discretization/structuredgrid.py @@ -741,7 +741,7 @@ def map_polygons(self): ############### ### Methods ### ############### - def intersect(self, x, y, local=False, forgive=False): + def intersect(self, x, y, z=None, local=False, forgive=False): """ Get the row and column of a point with coordinates x and y @@ -754,6 +754,9 @@ def intersect(self, x, y, local=False, forgive=False): The x-coordinate of the requested point y : float The y-coordinate of the requested point + z : float + Optional z-coordinate of the requested point (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) @@ -797,7 +800,28 @@ def intersect(self, x, y, local=False, forgive=False): row = np.where(ycomp)[0][-1] if np.any(np.isnan([row, col])): row = col = np.nan - return row, col + if z is not None: + return None, row, col + + if z is None: + return row, col + + lay = np.nan + for layer in range(self.__nlay): + if ( + self.top_botm[layer, row, col] + >= z + >= self.top_botm[layer + 1, row, col] + ): + lay = layer + break + + if np.any(np.isnan([lay, row, col])): + lay = row = col = np.nan + if not forgive: + raise Exception("point given is outside the model area") + + return lay, row, col def _cell_vert_list(self, i, j): """Get vertices for a single cell or sequence of i, j locations.""" diff --git a/flopy/discretization/unstructuredgrid.py b/flopy/discretization/unstructuredgrid.py index 0fc9e1d339..b21861d0fc 100644 --- a/flopy/discretization/unstructuredgrid.py +++ b/flopy/discretization/unstructuredgrid.py @@ -1,7 +1,11 @@ import os import copy import numpy as np + +from matplotlib.path import Path + from .grid import Grid, CachedData +from ..utils.geometry import is_clockwise class UnstructuredGrid(Grid): @@ -506,9 +510,74 @@ def map_polygons(self): return copy.copy(self._polygons) - def intersect(self, x, y, local=False, forgive=False): - x, y = super().intersect(x, y, local, forgive) - raise Exception("Not implemented yet") + def intersect(self, x, y, z=None, local=False, forgive=False): + """ + Get the CELL2D number of a point with coordinates x and y + + When the point is on the edge of two cells, the cell with the lowest + CELL2D number is returned. + + Parameters + ---------- + x : float + The x-coordinate of the requested point + y : float + The y-coordinate of the requested point + z : float, None + optional, z-coordiante of the requested point + 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 + ------- + icell2d : int + The CELL2D number + + """ + if local: + # transform x and y to real-world coordinates + x, y = super().get_coords(x, y) + xv, yv, zv = self.xyzvertices + + if self.grid_varies_by_layer: + ncpl = self.nnodes + else: + ncpl = self.ncpl[0] + + for icell2d in range(ncpl): + xa = np.array(xv[icell2d]) + ya = np.array(yv[icell2d]) + # x and y at least have to be within the bounding box of the cell + if ( + np.any(x <= xa) + and np.any(x >= xa) + and np.any(y <= ya) + and np.any(y >= ya) + ): + if is_clockwise(xa, ya): + radius = -1e-9 + else: + radius = 1e-9 + path = Path(np.stack((xa, ya)).transpose()) + # use a small radius, so that the edge of the cell is included + if path.contains_point((x, y), radius=radius): + if z is None: + return icell2d + + for lay in range(self.nlay): + if lay != 0 and not self.grid_varies_by_layer: + icell2d += self.ncpl[lay - 1] + if zv[0, icell2d] >= z >= zv[1, icell2d]: + return icell2d + + if forgive: + icell2d = np.nan + return icell2d + + raise Exception("point given is outside of the model area") @property def top_botm(self): diff --git a/flopy/discretization/vertexgrid.py b/flopy/discretization/vertexgrid.py index 6eab057007..dc19923219 100644 --- a/flopy/discretization/vertexgrid.py +++ b/flopy/discretization/vertexgrid.py @@ -254,7 +254,7 @@ def map_polygons(self): return copy.copy(self._polygons) - def intersect(self, x, y, local=False, forgive=False): + def intersect(self, x, y, z=None, local=False, forgive=False): """ Get the CELL2D number of a point with coordinates x and y @@ -267,6 +267,9 @@ def intersect(self, x, y, local=False, forgive=False): The x-coordinate of the requested point y : float The y-coordinate of the requested point + z : float, None + optional, z-coordiante of the requested point will return + (lay, icell2d) local: bool (optional) If True, x and y are in local coordinates (defaults to False) forgive: bool (optional) @@ -301,11 +304,25 @@ def intersect(self, x, y, local=False, forgive=False): else: radius = 1e-9 if path.contains_point((x, y), radius=radius): - return icell2d + if z is None: + return icell2d + + for lay in range(self.nlay): + if ( + self.top_botm[lay, icell2d] + >= z + >= self.top_botm[lay + 1, icell2d] + ): + return lay, icell2d + if forgive: icell2d = np.nan + if z is not None: + return np.nan, icell2d + return icell2d - raise Exception("x, y point given is outside of the model area") + + raise Exception("point given is outside of the model area") def get_cell_vertices(self, cellid): """ From b012b9e4a177ed876a9ec11ec992089c2a5294eb Mon Sep 17 00:00:00 2001 From: Joshua Larsen Date: Wed, 12 Jan 2022 09:50:28 -0800 Subject: [PATCH 2/3] update(Grid): Added interface change warnings to intersect method --- flopy/discretization/grid.py | 24 ++++++++++++++++++++++++ flopy/discretization/structuredgrid.py | 5 +++++ flopy/discretization/unstructuredgrid.py | 4 ++++ flopy/discretization/vertexgrid.py | 3 +++ 4 files changed, 36 insertions(+) diff --git a/flopy/discretization/grid.py b/flopy/discretization/grid.py index 837b80a5cc..bb3f699f82 100644 --- a/flopy/discretization/grid.py +++ b/flopy/discretization/grid.py @@ -671,6 +671,30 @@ def intersect(self, x, y, local=False, forgive=False): else: return x, y + def _warn_intersect(self, module, lineno): + """ + Warning for modelgrid intersect() interface change. + + Should be be removed after a couple of releases. Added in 3.3.5 + + Parameters + ---------- + module : str + module name path + lineno : int + line number where warning is called from + + Returns + ------- + None + """ + module = os.path.split(module)[-1] + warning = "The interface 'intersect(self, x, y, local=False, " \ + "forgive=False)' has been deprecated. Use the " \ + "intersect(self, x, y, z=None, local=False, " \ + "forgive=False) interface instead." + warnings.warn_explicit(warning, UserWarning, module, lineno) + def set_coord_info( self, xoff=None, diff --git a/flopy/discretization/structuredgrid.py b/flopy/discretization/structuredgrid.py index fca33a609d..b7153f9368 100644 --- a/flopy/discretization/structuredgrid.py +++ b/flopy/discretization/structuredgrid.py @@ -1,5 +1,6 @@ import copy import os.path +import inspect import numpy as np from .grid import Grid, CachedData @@ -771,6 +772,10 @@ def intersect(self, x, y, z=None, local=False, forgive=False): The column number """ + # trigger interface change warning + frame_info = inspect.getframeinfo(inspect.currentframe()) + self._warn_intersect(frame_info.filename, frame_info.lineno) + # transform x and y to local coordinates x, y = super().intersect(x, y, local, forgive) diff --git a/flopy/discretization/unstructuredgrid.py b/flopy/discretization/unstructuredgrid.py index b21861d0fc..6d0adb3f31 100644 --- a/flopy/discretization/unstructuredgrid.py +++ b/flopy/discretization/unstructuredgrid.py @@ -1,6 +1,7 @@ import os import copy import numpy as np +import inspect from matplotlib.path import Path @@ -537,6 +538,9 @@ def intersect(self, x, y, z=None, local=False, forgive=False): The CELL2D number """ + frame_info = inspect.getframeinfo(inspect.currentframe()) + self._warn_intersect(frame_info.filename, frame_info.lineno) + if local: # transform x and y to real-world coordinates x, y = super().get_coords(x, y) diff --git a/flopy/discretization/vertexgrid.py b/flopy/discretization/vertexgrid.py index dc19923219..558dc10dc9 100644 --- a/flopy/discretization/vertexgrid.py +++ b/flopy/discretization/vertexgrid.py @@ -1,6 +1,7 @@ import os import copy import numpy as np +import inspect from matplotlib.path import Path @@ -282,6 +283,8 @@ def intersect(self, x, y, z=None, local=False, forgive=False): The CELL2D number """ + frame_info = inspect.getframeinfo(inspect.currentframe()) + self._warn_intersect(frame_info.filename, frame_info.lineno) if local: # transform x and y to real-world coordinates From 075795d4345b47e5160935074048612c770ca1f5 Mon Sep 17 00:00:00 2001 From: Joshua Larsen Date: Wed, 12 Jan 2022 09:54:07 -0800 Subject: [PATCH 3/3] Linting --- flopy/discretization/grid.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/flopy/discretization/grid.py b/flopy/discretization/grid.py index bb3f699f82..e09c837fb7 100644 --- a/flopy/discretization/grid.py +++ b/flopy/discretization/grid.py @@ -689,10 +689,12 @@ def _warn_intersect(self, module, lineno): None """ module = os.path.split(module)[-1] - warning = "The interface 'intersect(self, x, y, local=False, " \ - "forgive=False)' has been deprecated. Use the " \ - "intersect(self, x, y, z=None, local=False, " \ - "forgive=False) interface instead." + warning = ( + "The interface 'intersect(self, x, y, local=False, " + "forgive=False)' has been deprecated. Use the " + "intersect(self, x, y, z=None, local=False, " + "forgive=False) interface instead." + ) warnings.warn_explicit(warning, UserWarning, module, lineno) def set_coord_info(