From 99165c605ef0746c13dae29fee2268805d3f5b98 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Mon, 6 Apr 2026 15:40:44 -0700 Subject: [PATCH 1/7] New feats: updates to Grid, GridIntersect, and VoronoiGrid Updates support development of generalized HFB creation methods * Grid: add cell_area property via the shoelace algorithm * GridIntersect: add "experimental" support for UnstructuredGrid * VoronoiGrid: add support for get_disu6_gridprops() that can build MF6 DISU packages --- flopy/discretization/grid.py | 25 ++++++++++ flopy/utils/gridintersect.py | 27 ++++++++++- flopy/utils/voronoi.py | 90 ++++++++++++++++++++++++++++++++++-- 3 files changed, 138 insertions(+), 4 deletions(-) diff --git a/flopy/discretization/grid.py b/flopy/discretization/grid.py index 4c0b74e50f..8569c43822 100644 --- a/flopy/discretization/grid.py +++ b/flopy/discretization/grid.py @@ -411,6 +411,31 @@ 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.) + return np.ravel(area) + @property def cell_thickness(self): """ diff --git a/flopy/utils/gridintersect.py b/flopy/utils/gridintersect.py index f8c9f46c8a..a32a87b42b 100644 --- a/flopy/utils/gridintersect.py +++ b/flopy/utils/gridintersect.py @@ -106,6 +106,8 @@ def __init__(self, mfgrid, rtree=True, local=False): self.geoms, self.cellids = self._rect_grid_to_geoms_cellids() elif self.mfgrid.grid_type == "vertex": self.geoms, self.cellids = self._vtx_grid_to_geoms_cellids() + elif self.mfgrid.grid_type == "unstructured": + self.geoms, self.cellids = self._usg_grid_to_geoms_cellids() else: raise NotImplementedError( f"Grid type {self.mfgrid.grid_type} not supported" @@ -375,7 +377,30 @@ def _usg_grid_to_geoms_cellids(self): cellids : array_like array of cellids """ - raise NotImplementedError() + warnings.warn( + "UnstructuredGrid intersection is experimental", + category=UserWarning + ) + shapely = import_optional_dependency("shapely") + if self.local: + geoms = [ + shapely.polygons( + list( + zip( + *self.mfgrid.get_local_coords( + *np.array(self.mfgrid.get_cell_vertices(node)).T + ) + ) + ) + ) + for node in range(self.mfgrid.nnodes) + ] + else: + geoms = [ + shapely.polygons(self.mfgrid.get_cell_vertices(node)) + for node in range(self.mfgrid.nnodes) + ] + return np.array(geoms), np.arange(self.mfgrid.nnodes) def _vtx_grid_to_geoms_cellids(self): """internal method, return shapely polygons and cellids for vertex diff --git a/flopy/utils/voronoi.py b/flopy/utils/voronoi.py index c9085e2c46..f2a0dd60d5 100644 --- a/flopy/utils/voronoi.py +++ b/flopy/utils/voronoi.py @@ -4,7 +4,7 @@ import numpy as np from .cvfdutil import get_disv_gridprops -from .geometry import point_in_polygon +from .geometry import point_in_polygon, distance from .utl_import import import_optional_dependency @@ -309,8 +309,92 @@ def get_disu5_gridprops(self): raise NotImplementedError(msg) def get_disu6_gridprops(self): - msg = "This method is not implemented yet." - raise NotImplementedError(msg) + """ + Get a dictionary of arguments that can be passed in to the + flopy.mf6.ModflowGwfdisu class + + Returns + ------- + disu_gridprops : dict + Dictionary of arguments that can be unpacked into the + flopy.mf6.ModflowGwfdisu constructor + + """ + from ..discretization import UnstructuredGrid + + gridprops = self.get_gridprops_unstructuredgrid() + ugrid = UnstructuredGrid( + **gridprops + ) + + xcenters = ugrid.xcellcenters + ycenters = ugrid.ycellcenters + xvertices = ugrid.xvertices + yvertices = ugrid.yvertices + iverts = ugrid.iverts + neighbors = ugrid.neighbors() + iac = [] + ja = [] + ihc = [] + cl12 = [] + hwva = [] + angldegx = [] + for node, nnodes in neighbors.items(): + iac.append(len(nnodes) + 1) + ja.extend([node] + nnodes) + ihc.extend([1, ] * iac[-1]) + # cell center to cell center distance... + cl12.append(0) + hwva.append(0) + angldegx.append(0) + + xc = xcenters[node] + yc = ycenters[node] + nxc = xcenters[nnodes] + nyc = ycenters[nnodes] + dists = distance(nxc, nyc, xc, yc) + cl12.extend(dists) + + angles = np.arctan2(nxc - xc, nyc - yc) * 180. / np.pi + angles = np.where(angles < 0, angles + 360, angles) + angldegx.extend(angles) + + ivrts = iverts[node] + xverts = xvertices[node] + yverts = yvertices[node] + for n in nnodes: + nivrts = iverts[n] + common = set(ivrts) & set(nivrts) + idxs = [ivrts.index(i) for i in common] + hwv = distance( + xverts[idxs[0]], + yverts[idxs[0]], + xverts[idxs[1]], + yverts[idxs[1]] + ) + hwva.append(hwv) + + cell2d = [] + for ix, ivrts in enumerate(iverts): + rec = [ix, xcenters[ix], ycenters[ix], len(ivrts)] + ivrts + cell2d.append(rec) + + gridprops["iac"] = iac + gridprops["ja"] = ja + gridprops["ihc"] = ihc + gridprops["cl12"] = cl12 + gridprops["hwva"] = hwva + gridprops["angldegx"] = angldegx + gridprops["area"] = ugrid.cell_area + gridprops["nodes"] = len(iac) + gridprops["nja"] = len(ja) + gridprops["nvert"] = len(gridprops["vertices"]) + gridprops["cell2d"] = cell2d + gridprops.pop("ncpl") + gridprops.pop("iverts") + gridprops.pop("xcenters") + gridprops.pop("ycenters") + return gridprops def get_gridprops_vertexgrid(self): """ From 81dda33616ecb0e9988641828462c470a9af9ca6 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Mon, 6 Apr 2026 16:43:03 -0700 Subject: [PATCH 2/7] Feat(hfb_util): add hfb builder method * add make_hfb_array and supporting methods to build HFB recarrays from LineString / Grid intersections --- flopy/utils/__init__.py | 1 + flopy/utils/hfb_util.py | 286 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 287 insertions(+) create mode 100644 flopy/utils/hfb_util.py diff --git a/flopy/utils/__init__.py b/flopy/utils/__init__.py index bd15a59b86..0da0e3d107 100644 --- a/flopy/utils/__init__.py +++ b/flopy/utils/__init__.py @@ -30,6 +30,7 @@ get_modflow = get_modflow_module.run_main from .gridintersect import GridIntersect +from .hfb_util import make_hfb_array from .mflistfile import ( Mf6ListBudget, MfListBudget, diff --git a/flopy/utils/hfb_util.py b/flopy/utils/hfb_util.py new file mode 100644 index 0000000000..7fc124c42c --- /dev/null +++ b/flopy/utils/hfb_util.py @@ -0,0 +1,286 @@ +import numpy as np +import pandas as pd +from .geometry import distance +from .gridintersect import GridIntersect +from .geospatial_utils import GeoSpatialUtil + + +def _min_distance_index(x0, x1, y0, y1): + """ + Method to get the index of the minimum distance between points + + Parameters + ---------- + x0 : np.ndarray + array of x coordinates + x1 : float + x coordinate of a point + y0 : np.ndarray + array of y coordinates + y1 : float + y coordinate of a point + + Returns + ------- + int : index location + """ + dist = distance(x0, y0, x1, y1) + idx = np.where(dist == np.min(dist))[0] + return idx[0] + + +def _edge_length_lut(xyverts): + """ + Method to create a look-up table based on index of + grid cell edge lengths + + Parameters + ---------- + xyverts: np.ndarray + numpy array of xy coordinate pairs + + Returns + ------- + dict + """ + xy0 = xyverts.T[:, :-1] + xy1 = xyverts.T[:, 1:] + dist = distance(xy0[0], xy0[1], xy1[0], xy1[1]) + dist_lu = {} + for v, d in enumerate(dist): + if v < len(dist) - 1: + dist_lu[(v, v + 1)] = d + else: + dist_lu[(0, v)] = d + return dist_lu + + +def _minimize_hfb_deviance(idxs, xyverts, pts): + """ + Method to minimumize the distance of the HFB fault trace from + the intersection points of the fault and model grid cell. Used + primarily for tie-breakers where there is not a clear "routing" + option based on vertex to vertex edge length distance + + Parameters + ---------- + idxs : + xyverts : + pts : + + Returns + ------- + float : sum of minimum distances between edge centroids across a potential + fault routing option and the intersection points between the fault line + and grid node + """ + lcs = [] + xyverts = xyverts.T + pts = pts.T + for ix in range(1, len(idxs)): + ixx = (idxs[ix - 1], idxs[ix]) + xc = np.mean(xyverts[0, ixx]) + yc = np.mean(xyverts[1, ixx]) + lcs.append([xc, yc]) + + mins = [] + for xc, yc in lcs: + dists = distance(pts[0], pts[1], xc, yc) + # asq = (pts[0] - xc) ** 2 + # bsq = (pts[1] - yc) ** 2 + # dists = np.sqrt(asq + bsq) + mins.append(np.min(dists)) + + return np.sum(mins) + + +def _edge_neighbors(modelgrid): + node_num = 0 + geoms = [] + node_nums = [] + + for poly in modelgrid.iverts: + poly = [int(i) for i in poly] + if poly[0] == poly[-1]: + poly = poly[:-1] + for v in range(len(poly)): + geoms.append(tuple(sorted([poly[v - 1], poly[v]]))) + node_nums += [node_num] * len(poly) + node_num += 1 + + edge_nodes = {} + for i, item in enumerate(geoms): + if item not in edge_nodes: + edge_nodes[item] = set([node_nums[i],]) + else: + edge_nodes[item].add(node_nums[i]) + + return edge_nodes + + +def make_hfb_array(modelgrid, geom): + """ + Method to make a HFB recarray from geospatial linestring information. + Note: this method was developed to only accept a single linestring at a + time to give the user an opportunity to save unique fault information + for calibration purposes. + + Parameters + ---------- + modelgrid : flopy.discretization.Grid object + FloPy StructuredGrid, VertexGrid, and UnstructuredGrid are supported + geom : geospatial object + geom parameter is a geospatial LineString representation. + Shapely.geometry.LineString, GeoJson, List of vertices, + shapefile.shape types are supported. + + Returns + ------- + np.recarray + numpy recarray of cellid1, cellid2, and hydchr. The hydchr field must + be set by the user, is set to NaN, and is provided for convenience. + """ + gu = GeoSpatialUtil(geom, "LineString") + if gu.shapetype.lower() != "linestring": + raise AssertionError( + f"{gu.shapetype} is not supported, only LineStrings are supported" + ) + + geom = gu.points + + if modelgrid.idomain is not None: + idomain = modelgrid.idomain.ravel() + else: + idomain = np.ones((modelgrid.nnodes,), dtype=int) + edge_set = _edge_neighbors(modelgrid) + vert_ivert = {tuple(i): cnt for cnt, i in enumerate(modelgrid.verts)} + iverts = modelgrid.iverts + xverts, yverts = modelgrid.cross_section_vertices + verts = [] + for node, xv in enumerate(xverts): + yv = yverts[node] + vrts = list(zip(xv, yv)) + if len(vrts) == len(iverts[node]): + if iverts[node][0] != iverts[node][-1]: + vrts.append(vrts[0]) + verts.append(np.array(vrts)) + + + ixs = GridIntersect(modelgrid) + result = ixs.intersect(geom, shapetype="LineString") + result_adj = [] + for record in result: + node = record.cellid + ixshp = record.ixshapes + coords = ixshp.coords.xy + x0, x1 = coords[0][0], coords[0][-1] + y0, y1 = coords[1][0], coords[1][-1] + + xycell = verts[node] + xcell = xycell.T[0][:-1] + ycell = xycell.T[1][:-1] + + vidx0 = _min_distance_index(xcell, x0, ycell, y0) + vidx1 = _min_distance_index(xcell, x1, ycell, y1) + + if vidx0 == vidx1: + continue + + tmp = tuple(sorted([vidx0, vidx1])) + + if tmp[1] - tmp[0] > 1: + nvert = len(xycell) - 1 + elens = _edge_length_lut(xycell) + # construct line routing options + o1 = list(range(tmp[0], tmp[1] + 1)) + o2 = list(range(tmp[1], nvert)) + list(range(0, tmp[0] + 1)) + # calculate routing distance + d1 = np.sum([elens[tuple(sorted([o1[ix - 1], o1[ix]]))] for ix in range(1, len(o1))]) + d2 = np.sum([elens[tuple(sorted([o2[ix - 1], o2[ix]]))] for ix in range(1, len(o2))]) + + # evaluate distance and break ties if necessary + if d1 < d2: + tmp = o1 + elif d2 < d1: + tmp = o2 + else: + om1 = _minimize_hfb_deviance(o1, xycell, np.array(coords)) + om2 = _minimize_hfb_deviance(o2, xycell, np.array(coords)) + if om1 <= om2: + tmp = o1 + else: + tmp = o2 + + edges = [] + for ix in range(1, len(tmp)): + eix0 = tmp[ix - 1] + eix1 = tmp[ix] + + xyv0 = tuple(xycell[eix0]) + xyv1 = tuple(xycell[eix1]) + + iv0 = vert_ivert[xyv0] + iv1 = vert_ivert[xyv1] + + edges.append(tuple(sorted([iv0, iv1]))) + + hfb_neighs = [] + for edge in edges: + for n in edge_set[edge]: + if n == node: + continue + hfb_neighs.append(n) + + res = [int(record.cellid), hfb_neighs] + result_adj.append(res) + + hfb_data = [] + visited = [] + if modelgrid.nlay is not None and modelgrid.grid_type != "unstructured": + for lay in range(modelgrid.nlay): + ncpl_adj = lay * modelgrid.ncpl + for cid0, hfb_neighs in result_adj: + if not idomain[cid0 + ncpl_adj]: + continue + + if modelgrid.grid_type == "structured": + cellid0 = modelgrid.get_lrc(cid0 + ncpl_adj)[0] + else: + cellid0 = (lay, cid0) + + for cid1 in hfb_neighs: + if not idomain[cid1 + ncpl_adj]: + continue + + if modelgrid.grid_type == "structured": + cellid1 = modelgrid.get_lrc(cid1 + ncpl_adj)[0] + else: + cellid1 = (lay, cid1) + + if cellid0 + cellid1 in visited: + continue + elif cellid1 + cellid0 in visited: + continue + + visited.append(cellid0 + cellid1) + hfb_data.append((cellid0, cellid1)) + else: + for cellid0, hfb_neighs in result_adj: + if not idomain[cellid0]: + continue + + for cellid1 in hfb_neighs: + if not idomain[cellid1]: + continue + + if (cellid0, cellid1) in visited: + continue + elif (cellid1, cellid0) in visited: + continue + + visited.append((cellid0, cellid1)) + hfb_data.append(((cellid0,), (cellid1,))) + + df = pd.DataFrame(hfb_data, columns=["cellid1", "cellid2"]) + df["hydchr"] = np.nan + return df.to_records(index=False) \ No newline at end of file From fce60746053bfecb1704387b4259b9833436bbf8 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Mon, 6 Apr 2026 16:44:15 -0700 Subject: [PATCH 3/7] Linting --- flopy/discretization/grid.py | 6 ++---- flopy/utils/gridintersect.py | 3 +-- flopy/utils/hfb_util.py | 17 ++++++++++++----- flopy/utils/voronoi.py | 20 ++++++++++---------- 4 files changed, 25 insertions(+), 21 deletions(-) diff --git a/flopy/discretization/grid.py b/flopy/discretization/grid.py index 8569c43822..0a50597c9e 100644 --- a/flopy/discretization/grid.py +++ b/flopy/discretization/grid.py @@ -425,15 +425,13 @@ def cell_area(self): from ..plot.plotutil import UnstructuredPlotUtilities 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 diff --git a/flopy/utils/gridintersect.py b/flopy/utils/gridintersect.py index a32a87b42b..1304fe41da 100644 --- a/flopy/utils/gridintersect.py +++ b/flopy/utils/gridintersect.py @@ -378,8 +378,7 @@ def _usg_grid_to_geoms_cellids(self): array of cellids """ warnings.warn( - "UnstructuredGrid intersection is experimental", - category=UserWarning + "UnstructuredGrid intersection is experimental", category=UserWarning ) shapely = import_optional_dependency("shapely") if self.local: diff --git a/flopy/utils/hfb_util.py b/flopy/utils/hfb_util.py index 7fc124c42c..a73a0d73ed 100644 --- a/flopy/utils/hfb_util.py +++ b/flopy/utils/hfb_util.py @@ -111,7 +111,11 @@ def _edge_neighbors(modelgrid): edge_nodes = {} for i, item in enumerate(geoms): if item not in edge_nodes: - edge_nodes[item] = set([node_nums[i],]) + edge_nodes[item] = set( + [ + node_nums[i], + ] + ) else: edge_nodes[item].add(node_nums[i]) @@ -165,7 +169,6 @@ def make_hfb_array(modelgrid, geom): vrts.append(vrts[0]) verts.append(np.array(vrts)) - ixs = GridIntersect(modelgrid) result = ixs.intersect(geom, shapetype="LineString") result_adj = [] @@ -195,8 +198,12 @@ def make_hfb_array(modelgrid, geom): o1 = list(range(tmp[0], tmp[1] + 1)) o2 = list(range(tmp[1], nvert)) + list(range(0, tmp[0] + 1)) # calculate routing distance - d1 = np.sum([elens[tuple(sorted([o1[ix - 1], o1[ix]]))] for ix in range(1, len(o1))]) - d2 = np.sum([elens[tuple(sorted([o2[ix - 1], o2[ix]]))] for ix in range(1, len(o2))]) + d1 = np.sum( + [elens[tuple(sorted([o1[ix - 1], o1[ix]]))] for ix in range(1, len(o1))] + ) + d2 = np.sum( + [elens[tuple(sorted([o2[ix - 1], o2[ix]]))] for ix in range(1, len(o2))] + ) # evaluate distance and break ties if necessary if d1 < d2: @@ -283,4 +290,4 @@ def make_hfb_array(modelgrid, geom): df = pd.DataFrame(hfb_data, columns=["cellid1", "cellid2"]) df["hydchr"] = np.nan - return df.to_records(index=False) \ No newline at end of file + return df.to_records(index=False) diff --git a/flopy/utils/voronoi.py b/flopy/utils/voronoi.py index f2a0dd60d5..a6e354747f 100644 --- a/flopy/utils/voronoi.py +++ b/flopy/utils/voronoi.py @@ -318,14 +318,12 @@ def get_disu6_gridprops(self): disu_gridprops : dict Dictionary of arguments that can be unpacked into the flopy.mf6.ModflowGwfdisu constructor - + """ from ..discretization import UnstructuredGrid gridprops = self.get_gridprops_unstructuredgrid() - ugrid = UnstructuredGrid( - **gridprops - ) + ugrid = UnstructuredGrid(**gridprops) xcenters = ugrid.xcellcenters ycenters = ugrid.ycellcenters @@ -342,7 +340,12 @@ def get_disu6_gridprops(self): for node, nnodes in neighbors.items(): iac.append(len(nnodes) + 1) ja.extend([node] + nnodes) - ihc.extend([1, ] * iac[-1]) + ihc.extend( + [ + 1, + ] + * iac[-1] + ) # cell center to cell center distance... cl12.append(0) hwva.append(0) @@ -355,7 +358,7 @@ def get_disu6_gridprops(self): dists = distance(nxc, nyc, xc, yc) cl12.extend(dists) - angles = np.arctan2(nxc - xc, nyc - yc) * 180. / np.pi + angles = np.arctan2(nxc - xc, nyc - yc) * 180.0 / np.pi angles = np.where(angles < 0, angles + 360, angles) angldegx.extend(angles) @@ -367,10 +370,7 @@ def get_disu6_gridprops(self): common = set(ivrts) & set(nivrts) idxs = [ivrts.index(i) for i in common] hwv = distance( - xverts[idxs[0]], - yverts[idxs[0]], - xverts[idxs[1]], - yverts[idxs[1]] + xverts[idxs[0]], yverts[idxs[0]], xverts[idxs[1]], yverts[idxs[1]] ) hwva.append(hwv) From 70a66b86aa930fac9fde8e15d1defdd850171a42 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Tue, 7 Apr 2026 11:47:10 -0700 Subject: [PATCH 4/7] update hfb_util: add perturbation routine for edge cases * catches "mid cell split" and "colinear with cell boundary" edge cases --- autotest/test_hfb_util.py | 9 ++++ flopy/utils/hfb_util.py | 103 +++++++++++++++++++++++++++++++------- flopy/utils/voronoi.py | 2 +- 3 files changed, 95 insertions(+), 19 deletions(-) create mode 100644 autotest/test_hfb_util.py diff --git a/autotest/test_hfb_util.py b/autotest/test_hfb_util.py new file mode 100644 index 0000000000..7119ca817f --- /dev/null +++ b/autotest/test_hfb_util.py @@ -0,0 +1,9 @@ +import numpy as np +import pytest + +from modflow_devtools.markers import requires_exe + + +@pytest.fixture +def test_test(function_tmpdir): + pass diff --git a/flopy/utils/hfb_util.py b/flopy/utils/hfb_util.py index a73a0d73ed..8e776dc9d3 100644 --- a/flopy/utils/hfb_util.py +++ b/flopy/utils/hfb_util.py @@ -1,8 +1,9 @@ import numpy as np import pandas as pd + from .geometry import distance -from .gridintersect import GridIntersect from .geospatial_utils import GeoSpatialUtil +from .gridintersect import GridIntersect def _min_distance_index(x0, x1, y0, y1): @@ -26,7 +27,7 @@ def _min_distance_index(x0, x1, y0, y1): """ dist = distance(x0, y0, x1, y1) idx = np.where(dist == np.min(dist))[0] - return idx[0] + return idx def _edge_length_lut(xyverts): @@ -64,9 +65,12 @@ def _minimize_hfb_deviance(idxs, xyverts, pts): Parameters ---------- - idxs : - xyverts : - pts : + idxs : iterable + list, tuple, or numpy array of polygon vertex indices for test hfb path + xyverts : np.array + numpy array of (x,y) polygon vertices + pts : np.array + numpy array of (x,y) points of intersection between line and polygon Returns ------- @@ -86,15 +90,64 @@ def _minimize_hfb_deviance(idxs, xyverts, pts): mins = [] for xc, yc in lcs: dists = distance(pts[0], pts[1], xc, yc) - # asq = (pts[0] - xc) ** 2 - # bsq = (pts[1] - yc) ** 2 - # dists = np.sqrt(asq + bsq) mins.append(np.min(dists)) return np.sum(mins) +def perturb_intersection_coords(idxs, xyverts, ipt, epsilon=1e-06): + """ + + Parameters + ---------- + idxs : + xyverts : + coords : + cidx : + + Returns + ------- + + """ + xyverts = xyverts.T + xverts = xyverts[0][idxs] + yverts = xyverts[1][idxs] + + if (xverts[0] - xverts[1]) == 0: + # vertical line + new_vrt = [xverts[0], np.mean(yverts) - epsilon] + elif (yverts[0] - yverts[1]) == 0: + # horizontal line + new_vrt = [np.mean(xverts) - epsilon, yverts[0]] + else: + # need to adjust across the line + m = (yverts[1] - yverts[0]) / [xverts[1] - xverts[0]] + if xverts[0] != 0: + vidx = 0 + else: + vidx = 1 + b = yverts[vidx] / (m * xverts[vidx]) + cx = ipt[0] - epsilon + cy = m * cx + b + new_vrt = [cx, cy] + + ipt = new_vrt + return ipt + + def _edge_neighbors(modelgrid): + """ + Method to get a dictionary of unique node edges (by ivert) and the nodes the + edge is shared between + + Parameters + ---------- + modelgrid : flopy.discretization.Grid object + + Returns + ------- + dict : dictionary of {edge iverts : [nodes]} + """ node_num = 0 geoms = [] node_nums = [] @@ -111,11 +164,9 @@ def _edge_neighbors(modelgrid): edge_nodes = {} for i, item in enumerate(geoms): if item not in edge_nodes: - edge_nodes[item] = set( - [ - node_nums[i], - ] - ) + edge_nodes[item] = { + node_nums[i], + } else: edge_nodes[item].add(node_nums[i]) @@ -175,16 +226,32 @@ def make_hfb_array(modelgrid, geom): for record in result: node = record.cellid ixshp = record.ixshapes - coords = ixshp.coords.xy - x0, x1 = coords[0][0], coords[0][-1] - y0, y1 = coords[1][0], coords[1][-1] + # todo: numpy array this and Transpose + coords = np.array(ixshp.coords.xy).T + x0, y0 = coords[0, 0], coords[0, 1] + x1, y1 = coords[-1, 0], coords[-1, 1] xycell = verts[node] xcell = xycell.T[0][:-1] ycell = xycell.T[1][:-1] vidx0 = _min_distance_index(xcell, x0, ycell, y0) + if len(vidx0) > 1: + # perturb line by small epsilon + coords[0] = perturb_intersection_coords(vidx0, xycell, coords[0]) + x0, y0 = coords[0, 0], coords[0, 1] + vidx0 = _min_distance_index(xcell, x0, ycell, y0) + + vidx0 = vidx0[0] + vidx1 = _min_distance_index(xcell, x1, ycell, y1) + if len(vidx1) > 1: + # pertub line by small epsilon + coords[-1] = perturb_intersection_coords(vidx1, xycell, coords[-1]) + x1, y1 = coords[-1, 0], coords[-1, 1] + vidx1 = _min_distance_index(xcell, x1, ycell, y1) + + vidx1 = vidx1[0] if vidx0 == vidx1: continue @@ -211,8 +278,8 @@ def make_hfb_array(modelgrid, geom): elif d2 < d1: tmp = o2 else: - om1 = _minimize_hfb_deviance(o1, xycell, np.array(coords)) - om2 = _minimize_hfb_deviance(o2, xycell, np.array(coords)) + om1 = _minimize_hfb_deviance(o1, xycell, coords) + om2 = _minimize_hfb_deviance(o2, xycell, coords) if om1 <= om2: tmp = o1 else: diff --git a/flopy/utils/voronoi.py b/flopy/utils/voronoi.py index a6e354747f..e2ef04f49d 100644 --- a/flopy/utils/voronoi.py +++ b/flopy/utils/voronoi.py @@ -4,7 +4,7 @@ import numpy as np from .cvfdutil import get_disv_gridprops -from .geometry import point_in_polygon, distance +from .geometry import distance, point_in_polygon from .utl_import import import_optional_dependency From 9c2b1df0bf440279b7697771a4b55f635c0b7749 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Tue, 7 Apr 2026 14:10:11 -0700 Subject: [PATCH 5/7] Added tests for hfb_util --- autotest/test_hfb_util.py | 440 +++++++++++++++++++++++++++++++++++++- 1 file changed, 435 insertions(+), 5 deletions(-) diff --git a/autotest/test_hfb_util.py b/autotest/test_hfb_util.py index 7119ca817f..31f9958dfc 100644 --- a/autotest/test_hfb_util.py +++ b/autotest/test_hfb_util.py @@ -1,9 +1,439 @@ import numpy as np -import pytest - from modflow_devtools.markers import requires_exe +import flopy +from flopy.utils.hfb_util import make_hfb_array +from flopy.utils.triangle import Triangle +from flopy.utils.voronoi import VoronoiGrid + + +def structured_sim(): + lx = 100 + ly = 100 + nlay = 1 + nrow = 10 + ncol = 10 + delc = np.full((nrow,), ly / nrow) + delr = np.full((ncol,), lx / ncol) + top = np.full((nrow, ncol), 10) + botm = np.zeros((nlay, nrow, ncol)) + idomain = np.ones(botm.shape, dtype=int) + + sim = flopy.mf6.MFSimulation(sim_ws="tmp_struct") + ims = flopy.mf6.ModflowIms(sim, complexity="SIMPLE") + tdis = flopy.mf6.ModflowTdis(sim) + + gwf = flopy.mf6.ModflowGwf(sim, modelname="hfb_model") + dis = flopy.mf6.ModflowGwfdis( + gwf, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delc=delc, + delr=delr, + top=top, + botm=botm, + idomain=idomain, + ) + + return sim + + +def vertex_sim(path): + geom = [(0, 0), (100, 0), (100, 100), (0, 100), (0, 0)] + tri = Triangle(angle=30, model_ws=path) + tri.add_polygon(geom) + tri.add_region((5, 5), 0, maximum_area=40) + tri.build() + + vor = VoronoiGrid(tri) + pkg_props = vor.get_disv_gridprops() + ncpl = pkg_props["ncpl"] + nlay = 1 + top = np.full((ncpl,), 10) + botm = np.zeros((nlay, ncpl)) + idomain = np.ones(botm.shape, dtype=int) + + sim = flopy.mf6.MFSimulation(sim_ws="tmp_struct") + ims = flopy.mf6.ModflowIms(sim, complexity="SIMPLE") + tdis = flopy.mf6.ModflowTdis(sim) + + gwf = flopy.mf6.ModflowGwf(sim, modelname="hfb_model") + disv = flopy.mf6.ModflowGwfdisv( + gwf, nlay=1, top=top, botm=botm, idomain=idomain, **pkg_props + ) + + return sim + + +def unstructured_sim(path): + geom = [(0, 0), (100, 0), (100, 100), (0, 100), (0, 0)] + tri = Triangle(angle=30, model_ws=path) + tri.add_polygon(geom) + tri.add_region((5, 5), 0, maximum_area=40) + tri.build() + + vor = VoronoiGrid(tri) + grid_props = vor.get_disu6_gridprops() + + top = np.full((grid_props["nodes"],), 10) + botm = np.zeros((grid_props["nodes"])) + idomain = np.ones(botm.shape, dtype=int) + + sim = flopy.mf6.MFSimulation(sim_ws="tmp_struct") + ims = flopy.mf6.ModflowIms(sim, complexity="SIMPLE") + tdis = flopy.mf6.ModflowTdis(sim) + + gwf = flopy.mf6.ModflowGwf(sim, modelname="hfb_model") + disu = flopy.mf6.ModflowGwfdisu( + gwf, top=top, bot=botm, idomain=idomain, **grid_props + ) + + return sim + + +def test_simple_structured(): + validation = [ + (0, 0, 8, 0, 1, 8), + (0, 0, 8, 0, 0, 9), + (0, 1, 7, 0, 2, 7), + (0, 1, 7, 0, 1, 8), + (0, 2, 6, 0, 3, 6), + (0, 2, 6, 0, 2, 7), + (0, 3, 5, 0, 4, 5), + (0, 3, 5, 0, 3, 6), + (0, 4, 4, 0, 4, 5), + (0, 4, 4, 0, 5, 4), + (0, 5, 3, 0, 5, 4), + (0, 5, 3, 0, 6, 3), + (0, 6, 2, 0, 6, 3), + (0, 6, 2, 0, 7, 2), + (0, 7, 1, 0, 7, 2), + (0, 7, 1, 0, 8, 1), + (0, 8, 0, 0, 8, 1), + (0, 8, 0, 0, 9, 0), + ] + + fault = [(0, 4), (94.0, 100)] + sim = structured_sim() + gwf = sim.get_model() + modelgrid = gwf.modelgrid + hfbs = flopy.utils.make_hfb_array(modelgrid, fault) + for row in hfbs: + cid1 = row.cellid1 + cid2 = row.cellid2 + x = sorted([cid1, cid2]) + test_set = x[0] + x[1] + if test_set not in validation: + raise AssertionError( + f"HFB Line {x[0]} and {x[1]} are outside of validation data set" + ) + + if len(hfbs) != len(validation): + raise AssertionError( + f"HFB data length {len(hfbs)} not equal to validation {len(validation)}" + ) + + +@requires_exe("triangle") +def test_simple_vertex(function_tmpdir): + validation = [ + (0, 53, 0, 92), + (0, 65, 0, 162), + (0, 64, 0, 65), + (0, 69, 0, 71), + (0, 64, 0, 72), + (0, 45, 0, 77), + (0, 53, 0, 77), + (0, 78, 0, 79), + (0, 76, 0, 78), + (0, 78, 0, 162), + (0, 79, 0, 82), + (0, 11, 0, 79), + (0, 69, 0, 80), + (0, 45, 0, 80), + (0, 11, 0, 84), + (0, 69, 0, 84), + (0, 94, 0, 114), + (0, 4, 0, 94), + (0, 4, 0, 96), + (0, 97, 0, 101), + (0, 97, 0, 99), + (0, 98, 0, 99), + (0, 96, 0, 98), + (0, 108, 0, 130), + (0, 108, 0, 109), + (0, 101, 0, 108), + (0, 92, 0, 114), + (0, 107, 0, 130), + (0, 125, 0, 130), + (0, 125, 0, 132), + (0, 149, 0, 171), + (0, 133, 0, 149), + (0, 150, 0, 151), + (0, 151, 0, 171), + (0, 150, 0, 154), + (0, 125, 0, 169), + (0, 119, 0, 169), + (0, 149, 0, 169), + ] + + fault = [(0, 4), (94.0, 100)] + sim = vertex_sim(function_tmpdir) + gwf = sim.get_model() + modelgrid = gwf.modelgrid + hfbs = flopy.utils.make_hfb_array(modelgrid, fault) + for row in hfbs: + cid1 = row.cellid1 + cid2 = row.cellid2 + x = sorted([cid1, cid2]) + test_set = x[0] + x[1] + if test_set not in validation: + raise AssertionError( + f"HFB Line {x[0]} and {x[1]} are outside of validation data set" + ) + + if len(hfbs) != len(validation): + raise AssertionError( + f"HFB data length {len(hfbs)} not equal to validation {len(validation)}" + ) + + +@requires_exe("triangle") +def test_simple_unstructured(function_tmpdir): + validation = [ + (53, 92), + (65, 162), + (64, 65), + (69, 71), + (64, 72), + (45, 77), + (53, 77), + (78, 79), + (76, 78), + (78, 162), + (79, 82), + (11, 79), + (69, 80), + (45, 80), + (11, 84), + (69, 84), + (94, 114), + (4, 94), + (4, 96), + (97, 101), + (97, 99), + (98, 99), + (96, 98), + (108, 130), + (108, 109), + (101, 108), + (92, 114), + (107, 130), + (125, 130), + (125, 132), + (149, 171), + (133, 149), + (150, 151), + (151, 171), + (150, 154), + (125, 169), + (119, 169), + (149, 169), + ] + fault = [(0, 4), (94.0, 100)] + sim = unstructured_sim(function_tmpdir) + gwf = sim.get_model() + modelgrid = gwf.modelgrid + hfbs = flopy.utils.make_hfb_array(modelgrid, fault) + for row in hfbs: + cid1 = row.cellid1 + cid2 = row.cellid2 + x = sorted([cid1, cid2]) + test_set = x[0] + x[1] + if test_set not in validation: + raise AssertionError( + f"HFB Line {x[0]} and {x[1]} are outside of validation data set" + ) + + if len(hfbs) != len(validation): + raise AssertionError( + f"HFB data length {len(hfbs)} not equal to validation {len(validation)}" + ) + + +def test_multisegment_structured(): + validation = [ + (0, 0, 8, 0, 1, 8), + (0, 0, 8, 0, 0, 9), + (0, 1, 7, 0, 1, 8), + (0, 1, 7, 0, 2, 7), + (0, 2, 6, 0, 2, 7), + (0, 3, 6, 0, 4, 6), + (0, 3, 6, 0, 3, 7), + (0, 4, 5, 0, 4, 6), + (0, 4, 5, 0, 5, 5), + (0, 5, 4, 0, 5, 5), + (0, 6, 3, 0, 7, 3), + (0, 6, 3, 0, 6, 4), + (0, 5, 4, 0, 6, 4), + (0, 7, 1, 0, 8, 1), + (0, 6, 2, 0, 7, 2), + (0, 7, 1, 0, 7, 2), + (0, 8, 0, 0, 9, 0), + (0, 8, 0, 0, 8, 1), + ] + + fault = [[0, 11], [55, 45], [89, 100]] + sim = structured_sim() + gwf = sim.get_model() + modelgrid = gwf.modelgrid + hfbs = flopy.utils.make_hfb_array(modelgrid, fault) + for row in hfbs: + cid1 = row.cellid1 + cid2 = row.cellid2 + x = sorted([cid1, cid2]) + test_set = x[0] + x[1] + if test_set not in validation: + raise AssertionError( + f"HFB Line {x[0]} and {x[1]} are outside of validation data set" + ) + + if len(hfbs) != len(validation): + raise AssertionError( + f"HFB data length {len(hfbs)} not equal to validation {len(validation)}" + ) + + +@requires_exe("triangle") +def test_multisegment_vertex(function_tmpdir): + validation = [ + (0, 11, 0, 84), + (0, 11, 0, 79), + (0, 40, 0, 103), + (0, 4, 0, 40), + (0, 40, 0, 114), + (0, 40, 0, 53), + (0, 45, 0, 53), + (0, 45, 0, 77), + (0, 45, 0, 80), + (0, 29, 0, 53), + (0, 64, 0, 72), + (0, 64, 0, 65), + (0, 69, 0, 80), + (0, 69, 0, 71), + (0, 69, 0, 84), + (0, 78, 0, 79), + (0, 76, 0, 78), + (0, 78, 0, 162), + (0, 79, 0, 82), + (0, 97, 0, 110), + (0, 97, 0, 102), + (0, 98, 0, 102), + (0, 98, 0, 103), + (0, 4, 0, 103), + (0, 107, 0, 130), + (0, 107, 0, 108), + (0, 108, 0, 110), + (0, 125, 0, 132), + (0, 125, 0, 130), + (0, 133, 0, 149), + (0, 150, 0, 151), + (0, 150, 0, 154), + (0, 65, 0, 162), + (0, 125, 0, 169), + (0, 119, 0, 169), + (0, 149, 0, 169), + (0, 149, 0, 171), + (0, 151, 0, 171), + ] + + fault = [[0, 11], [55, 45], [89, 100]] + sim = vertex_sim(function_tmpdir) + gwf = sim.get_model() + modelgrid = gwf.modelgrid + hfbs = flopy.utils.make_hfb_array(modelgrid, fault) + for row in hfbs: + cid1 = row.cellid1 + cid2 = row.cellid2 + x = sorted([cid1, cid2]) + test_set = x[0] + x[1] + if test_set not in validation: + raise AssertionError( + f"HFB Line {x[0]} and {x[1]} are outside of validation data set" + ) + + if len(hfbs) != len(validation): + raise AssertionError( + f"HFB data length {len(hfbs)} not equal to validation {len(validation)}" + ) + + +def test_colinear_hfb(): + validation = [ + (0, 0, 4, 0, 0, 5), + (0, 1, 4, 0, 1, 5), + (0, 2, 4, 0, 2, 5), + (0, 3, 4, 0, 3, 5), + (0, 4, 4, 0, 4, 5), + (0, 5, 4, 0, 5, 5), + (0, 6, 4, 0, 6, 5), + (0, 7, 4, 0, 7, 5), + (0, 8, 4, 0, 8, 5), + (0, 9, 4, 0, 9, 5), + ] + + fault = [(50, 0), (50, 100)] + sim = structured_sim() + gwf = sim.get_model() + modelgrid = gwf.modelgrid + hfbs = flopy.utils.make_hfb_array(modelgrid, fault) + for row in hfbs: + cid1 = row.cellid1 + cid2 = row.cellid2 + x = sorted([cid1, cid2]) + test_set = x[0] + x[1] + if test_set not in validation: + raise AssertionError( + f"HFB Line {x[0]} and {x[1]} are outside of validation data set" + ) + + if len(hfbs) != len(validation): + raise AssertionError( + f"HFB data length {len(hfbs)} not equal to validation {len(validation)}" + ) + + +def test_split_cell_hfb(): + validation = [ + (0, 4, 0, 0, 5, 0), + (0, 4, 1, 0, 5, 1), + (0, 4, 2, 0, 5, 2), + (0, 4, 3, 0, 5, 3), + (0, 4, 4, 0, 5, 4), + (0, 4, 5, 0, 5, 5), + (0, 4, 6, 0, 5, 6), + (0, 4, 7, 0, 5, 7), + (0, 4, 8, 0, 5, 8), + (0, 4, 9, 0, 5, 9), + ] + + fault = [(0, 55), (100, 55)] + sim = structured_sim() + gwf = sim.get_model() + modelgrid = gwf.modelgrid + hfbs = flopy.utils.make_hfb_array(modelgrid, fault) + for row in hfbs: + cid1 = row.cellid1 + cid2 = row.cellid2 + x = sorted([cid1, cid2]) + test_set = x[0] + x[1] + if test_set not in validation: + raise AssertionError( + f"HFB Line {x[0]} and {x[1]} are outside of validation data set" + ) -@pytest.fixture -def test_test(function_tmpdir): - pass + if len(hfbs) != len(validation): + raise AssertionError( + f"HFB data length {len(hfbs)} not equal to validation {len(validation)}" + ) From 3c4e1901194ec30e2d157347f87e2e4c7ac59a1a Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Tue, 7 Apr 2026 14:12:30 -0700 Subject: [PATCH 6/7] fix spelling --- flopy/utils/hfb_util.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/flopy/utils/hfb_util.py b/flopy/utils/hfb_util.py index 8e776dc9d3..829e53d50b 100644 --- a/flopy/utils/hfb_util.py +++ b/flopy/utils/hfb_util.py @@ -246,7 +246,7 @@ def make_hfb_array(modelgrid, geom): vidx1 = _min_distance_index(xcell, x1, ycell, y1) if len(vidx1) > 1: - # pertub line by small epsilon + # perturb line by small epsilon coords[-1] = perturb_intersection_coords(vidx1, xycell, coords[-1]) x1, y1 = coords[-1, 0], coords[-1, 1] vidx1 = _min_distance_index(xcell, x1, ycell, y1) From a6cdb00ee9cd2fb035cb0cd4ab1dc277a1eff889 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Mon, 13 Apr 2026 16:57:05 -0700 Subject: [PATCH 7/7] Add notebook documentation for `make_hfb_array` utility --- .../horizontal_flow_barriers_from_linework.py | 305 ++++++++++++++++++ flopy/utils/hfb_util.py | 22 +- 2 files changed, 319 insertions(+), 8 deletions(-) create mode 100644 .docs/Notebooks/horizontal_flow_barriers_from_linework.py diff --git a/.docs/Notebooks/horizontal_flow_barriers_from_linework.py b/.docs/Notebooks/horizontal_flow_barriers_from_linework.py new file mode 100644 index 0000000000..6db2ed3cc0 --- /dev/null +++ b/.docs/Notebooks/horizontal_flow_barriers_from_linework.py @@ -0,0 +1,305 @@ +#!/usr/bin/env python +# coding: utf-8 + +# # Creating HFB input data from shapefile linework +# +# This notebook shows examples of how to create a horizontal flow barrier package (HFB) from shapefile information in FloPy. FloPy has support for creating HFB inputs from LineString like vector data for all discretization types (`DIS`, `DISV`, `DISU`). +# +# This notebook starts off by generating two simple models: +# - a structured (DIS) model +# - a vertex (DISV) model +# +# And then shows how to generate HFB input for each of these models. + +# In[1]: + + +from tempfile import TemporaryDirectory + +import geopandas as gpd +import matplotlib.pyplot as plt +import numpy as np +from shapely.geometry import LineString + +import flopy +from flopy.utils import make_hfb_array + +temp_dir = TemporaryDirectory() +workspace = temp_dir.name + + +# ### Generate a Structured (DIS) model + +# In[2]: + + +def simple_structured_model(): + """ + Method to generate a 10x10 structured grid model + """ + lx = 100 + ly = 100 + nlay = 1 + nrow = 10 + ncol = 10 + delc = np.full((nrow,), ly / nrow) + delr = np.full((ncol,), ly / ncol) + top = np.full((nrow, ncol), 10) + botm = np.zeros((nlay, nrow, ncol)) + idomain = np.ones(botm.shape, dtype=int) + + sim = flopy.mf6.MFSimulation(sim_ws="tmp_struct") + ims = flopy.mf6.ModflowIms(sim, complexity="SIMPLE") + tdis = flopy.mf6.ModflowTdis(sim) + + gwf = flopy.mf6.ModflowGwf(sim, modelname="hfb_model") + dis = flopy.mf6.ModflowGwfdis( + gwf, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delc=delc, + delr=delr, + top=top, + botm=botm, + idomain=idomain, + ) + npf = flopy.mf6.ModflowGwfnpf(gwf) + sto = flopy.mf6.ModflowGwfsto(gwf) + ic = flopy.mf6.ModflowGwfic(gwf, strt=np.full((nlay, nrow, ncol), 9)) + + return sim + + +# ### Generate a voronoi (DISV) model + +# In[3]: + + +def simple_vertex_model(workspace): + """ + Method to generate a voroni vertex grid model + """ + from flopy.utils.triangle import Triangle + from flopy.utils.voronoi import VoronoiGrid + + geom = [(0, 0), (100, 0), (100, 100), (0, 100), (0, 0)] + tri = Triangle(angle=30, model_ws=workspace) + tri.add_polygon(geom) + tri.add_region((5, 5), 0, maximum_area=40) + tri.build() + + vor = VoronoiGrid(tri) + pkg_props = vor.get_disv_gridprops() + ncpl = pkg_props["ncpl"] + nlay = 1 + top = np.full((ncpl,), 10) + botm = np.zeros((nlay, ncpl)) + idomain = np.ones(botm.shape, dtype=int) + + sim = flopy.mf6.MFSimulation(sim_ws="tmp_struct") + ims = flopy.mf6.ModflowIms(sim, complexity="SIMPLE") + tdis = flopy.mf6.ModflowTdis(sim) + + gwf = flopy.mf6.ModflowGwf(sim, modelname="hfb_model") + disv = flopy.mf6.ModflowGwfdisv( + gwf, nlay=1, top=top, botm=botm, idomain=idomain, **pkg_props + ) + npf = flopy.mf6.ModflowGwfnpf(gwf) + sto = flopy.mf6.ModflowGwfsto(gwf) + ic = flopy.mf6.ModflowGwfic(gwf, strt=np.full((nlay, ncpl), 9)) + + return sim + + +# ### Create HFB inputs for a Structured (DIS) model + +# The `make_hfb_array` utility is used to intersect `LineString` like data with `modelgrid` instances to produce HFB data. The `make_hfb_array` utility accepts two parameters: +# - `modelgrid` : a flopy Grid instance +# - `geom` : a geospatial LineString like object that can be: an iterable of points, GeoJSON, Shapely.geometry.LineString, or a shapefile.Shape object. + +# First define a fault line and visualize it with the model + +# In[4]: + + +# define a fault line +fault = [(0, 4), (50, 55), (100, 55)] + +# build simulation +sim = simple_structured_model() +gwf = sim.get_model() +modelgrid = gwf.modelgrid + +fig, ax = plt.subplots(figsize=(5, 5)) +pmv = flopy.plot.PlotMapView(modelgrid=modelgrid, ax=ax) +pmv.plot_grid() +plt.plot(np.array(fault).T[0], np.array(fault).T[1], "r--") +# now create hfb data with the `make_hfb_array` utility + +# In[5]: + + +hfbs = flopy.utils.make_hfb_array(modelgrid, fault) +hfbs[0:5] + + +# notice that the `hydchr` parameter is not filled out; this is left up to the user + +# In[6]: + + +hfbs["hydchr"] = 1e-05 + + +# Now we can build an HFB package and compare it to the fault line that was defined above + +# In[7]: + + +hfb = flopy.mf6.ModflowGwfhfb(gwf, stress_period_data={0: hfbs}) + + +# In[8]: + + +# plot up the results +fig, ax = plt.subplots(figsize=(5, 5)) +pmv = flopy.plot.PlotMapView(modelgrid=modelgrid) +pmv.plot_grid() +pmv.plot_bc(package=hfb, color="b", lw=2) +plt.plot(np.array(fault).T[0], np.array(fault).T[1], "r--") +# ### Create HFB inputs for a Vertex (DISV) model + +# The `make_hfb_array` utility is used to intersect `LineString` like data with `modelgrid` instances to produce HFB data. The `make_hfb_array` utility accepts two parameters: +# - `modelgrid` : a flopy Grid instance +# - `geom` : a geospatial LineString like object that can be: an iterable of points, GeoJSON, Shapely.geometry.LineString, or a shapefile.Shape object. + +# First define a fault line and visualize it with the model + +# In[9]: + + +# define a fault line +fault = [(0, 4), (50, 55), (100, 55)] + +# build simulation +sim = simple_vertex_model(workspace) +gwf = sim.get_model() +modelgrid = gwf.modelgrid + +fig, ax = plt.subplots(figsize=(5, 5)) +pmv = flopy.plot.PlotMapView(modelgrid=modelgrid, ax=ax) +pmv.plot_grid() +plt.plot(np.array(fault).T[0], np.array(fault).T[1], "r--") +# now create hfb data with the `make_hfb_array` utility + +# In[10]: + + +hfbs = flopy.utils.make_hfb_array(modelgrid, fault) +hfbs[0:5] + + +# notice that the `hydchr` parameter is not filled out; this is left up to the user + +# In[11]: + + +hfbs["hydchr"] = 3e-06 + + +# Now we can build an HFB package and compare it to the fault line that was defined above + +# In[12]: + + +hfb = flopy.mf6.ModflowGwfhfb(gwf, stress_period_data={0: hfbs}) + + +# In[13]: + + +# plot up the results +fig, ax = plt.subplots(figsize=(5, 5)) +pmv = flopy.plot.PlotMapView(modelgrid=modelgrid) +pmv.plot_grid() +pmv.plot_bc(package=hfb, color="b", lw=2) +plt.plot(np.array(fault).T[0], np.array(fault).T[1], "r--") +# ### Working with multiple fault geometries +# +# The `make_hfb_array` only accepts a single LineString geometry at a time. This was done by design for a couple reasons. +# 1) This gives the user an opportunity to set the `hydchr` parameter to a unique value for each fault line +# 2) This also gives the user an opportunity to save the boundary condition data (cellid1, cellid2) to an external file for each fault line that can be useful during model calibration with automated software (e.g., PEST++) +# +# As a result of these design descisions recarray data will need to be concatenated when working with multiple LineString geometries. + +# For this example, a geodataframe of two fault lines is created and then used to build a HFB package + +# In[14]: + + +geom1 = LineString([(0, 4), (50, 44)]) +geom2 = LineString([(50, 44), (88, 100)]) +d = {"name": ["fault1", "fault2"], "geometry": [geom1, geom2]} +gdf = gpd.GeoDataFrame(d) +gdf + + +# Visualize the faults on the modelgrid + +# In[15]: + + +# build simulation +sim = simple_structured_model() +gwf = sim.get_model() +modelgrid = gwf.modelgrid + +fig, ax = plt.subplots(figsize=(5, 5)) +pmv = flopy.plot.PlotMapView(modelgrid=modelgrid, ax=ax) +pmv.plot_grid() +gdf.geometry.plot(color="r", linestyle="--", ax=ax) +# now create hfb data with the `make_hfb_array` utility and concatenate it + +# In[16]: + + +hydchrs = {"fault1": 1e-05, "fault2": 1e-06} + +hfb_ras = [] +for name, geom in zip(gdf.name, gdf.geometry): + hfbs = make_hfb_array(modelgrid, geom) + hfbs["hydchr"] = hydchrs[name] + hfb_ras.append(hfbs) + +hfbs = np.concatenate(hfb_ras) +hfbs = hfbs.view(np.recarray) +hfbs + + +# Now create the hfb package from this data and visualize it for comparison with the fault lines + +# In[17]: + + +hfb = flopy.mf6.ModflowGwfhfb(gwf, stress_period_data={0: hfbs}) + + +# In[18]: + + +# plot up the results +fig, ax = plt.subplots(figsize=(5, 5)) +pmv = flopy.plot.PlotMapView(modelgrid=modelgrid) +pmv.plot_grid() +pmv.plot_bc(package=hfb, color="b", lw=2) +gdf.geometry.plot(color="r", linestyle="--", ax=ax) +# In[19]: + + +try: + temp_dir.cleanup() +except PermissionError: + # can occur on windows: https://docs.python.org/3/library/tempfile.html#tempfile.TemporaryDirectory + pass diff --git a/flopy/utils/hfb_util.py b/flopy/utils/hfb_util.py index 829e53d50b..97d6c83983 100644 --- a/flopy/utils/hfb_util.py +++ b/flopy/utils/hfb_util.py @@ -95,19 +95,25 @@ def _minimize_hfb_deviance(idxs, xyverts, pts): return np.sum(mins) -def perturb_intersection_coords(idxs, xyverts, ipt, epsilon=1e-06): +def _perturb_intersection_coords(idxs, xyverts, ipt, epsilon=1e-06): """ + Method to perturb an intersection location by a small amount which is used + to handle colinear intersection and intersection at the midpoint of a cell edge. Parameters ---------- - idxs : - xyverts : - coords : - cidx : + idxs : list + list of vertex indices + xyverts : np.array + numpy array of x,y vertices for a cell + ipt : iterable + intersection point x,y coordinate + epsilon : float + perturbation value Returns ------- - + ipt : list holding x, y coordinate pair of the perturbed intersection """ xyverts = xyverts.T xverts = xyverts[0][idxs] @@ -238,7 +244,7 @@ def make_hfb_array(modelgrid, geom): vidx0 = _min_distance_index(xcell, x0, ycell, y0) if len(vidx0) > 1: # perturb line by small epsilon - coords[0] = perturb_intersection_coords(vidx0, xycell, coords[0]) + coords[0] = _perturb_intersection_coords(vidx0, xycell, coords[0]) x0, y0 = coords[0, 0], coords[0, 1] vidx0 = _min_distance_index(xcell, x0, ycell, y0) @@ -247,7 +253,7 @@ def make_hfb_array(modelgrid, geom): vidx1 = _min_distance_index(xcell, x1, ycell, y1) if len(vidx1) > 1: # perturb line by small epsilon - coords[-1] = perturb_intersection_coords(vidx1, xycell, coords[-1]) + coords[-1] = _perturb_intersection_coords(vidx1, xycell, coords[-1]) x1, y1 = coords[-1, 0], coords[-1, 1] vidx1 = _min_distance_index(xcell, x1, ycell, y1)