diff --git a/autotest/test_plot_cross_section.py b/autotest/test_plot_cross_section.py index 9ab96480b0..9fd3770551 100644 --- a/autotest/test_plot_cross_section.py +++ b/autotest/test_plot_cross_section.py @@ -1,7 +1,7 @@ import matplotlib.pyplot as plt import numpy as np import pytest -from matplotlib.collections import PatchCollection +from matplotlib.collections import LineCollection, PatchCollection from modflow_devtools.markers import requires_pkg import flopy @@ -259,3 +259,236 @@ def test_plot_centers(): xmax = np.max(verts[0]) if xmax < center < xmin: raise AssertionError("Cell center not properly drawn on cross-section") + + +@pytest.fixture(params=["dis", "disv", "disu"]) +def hfb_xc_model(request): + """Create a MODFLOW 6 model with HFB for cross section testing. + + Parameters + ---------- + request.param : str + Grid type: "dis" (structured), "disv" (vertex), or "disu" (unstructured) + + Returns + ------- + tuple + (gwf_model, cross_section_line, grid_type) + """ + from flopy.utils.gridutil import get_disu_kwargs, get_disv_kwargs + + grid_type = request.param + + # Create simulation + sim = flopy.mf6.MFSimulation(sim_name=f"test_hfb_xc_{grid_type}") + tdis = flopy.mf6.ModflowTdis(sim) + ims = flopy.mf6.ModflowIms(sim) + gwf = flopy.mf6.ModflowGwf(sim, modelname="test") + + # Create discretization based on grid type + if grid_type == "dis": + # Structured grid with 2 layers + dis = flopy.mf6.ModflowGwfdis( + gwf, + nlay=2, + nrow=10, + ncol=10, + delr=100.0, + delc=100.0, + top=100.0, + botm=[50.0, 0.0], + ) + # HFB cellids for structured grid: (layer, row, col) + # Create barriers along column boundary at row 4 + hfb_data = [ + [(0, 3, 4), (0, 3, 5), 1e-6], + [(0, 4, 4), (0, 4, 5), 1e-6], + [(0, 5, 4), (0, 5, 5), 1e-6], + ] + xc_line = {"row": 4} + + elif grid_type == "disv": + # Vertex grid (regular grid converted to vertex format) + disv_kwargs = get_disv_kwargs(2, 10, 10, 100.0, 100.0, 100.0, [50.0, 0.0]) + disv = flopy.mf6.ModflowGwfdisv(gwf, **disv_kwargs) + # HFB cellids for vertex grid: (layer, cell2d_id) + # For a 10x10 grid, cell2d IDs are row*ncol + col + hfb_data = [ + [(0, 3 * 10 + 4), (0, 3 * 10 + 5), 1e-6], + [(0, 4 * 10 + 4), (0, 4 * 10 + 5), 1e-6], + [(0, 5 * 10 + 4), (0, 5 * 10 + 5), 1e-6], + ] + # For DISV, use line coordinates instead of row + xc_line = {"line": ([0, 450], [1000, 450])} + + elif grid_type == "disu": + # Unstructured grid (regular grid converted to unstructured format) + disu_kwargs = get_disu_kwargs( + 2, 10, 10, 100.0, 100.0, 100.0, [50.0, 0.0], return_vertices=True + ) + disu = flopy.mf6.ModflowGwfdisu(gwf, **disu_kwargs) + # HFB cellids for unstructured grid: (node,) + # For a 10x10 grid, node = layer * nrow * ncol + row * ncol + col + hfb_data = [ + [(0 * 10 * 10 + 3 * 10 + 4,), (0 * 10 * 10 + 3 * 10 + 5,), 1e-6], + [(0 * 10 * 10 + 4 * 10 + 4,), (0 * 10 * 10 + 4 * 10 + 5,), 1e-6], + [(0 * 10 * 10 + 5 * 10 + 4,), (0 * 10 * 10 + 5 * 10 + 5,), 1e-6], + ] + # For DISU, use line coordinates + xc_line = {"line": ([0, 450], [1000, 450])} + + # Add packages + ic = flopy.mf6.ModflowGwfic(gwf, strt=75.0) + npf = flopy.mf6.ModflowGwfnpf(gwf, save_flows=True) + hfb = flopy.mf6.ModflowGwfhfb(gwf, stress_period_data=hfb_data) + + return gwf, xc_line, grid_type + + +def test_cross_section_bc_hfb(hfb_xc_model): + """Test plotting HFB (Horizontal Flow Barrier) in cross sections. + + HFB packages have cellid1/cellid2 fields instead of a single cellid field. + Plot barriers by showing both cells that the barrier affects. + + Tests structured (DIS), vertex (DISV), and unstructured (DISU) grids. + + Addresses issue #2676. + """ + gwf, xc_line, grid_type = hfb_xc_model + + # Create figure and cross section + fig, ax = plt.subplots(1, 1, figsize=(10, 5)) + xc = flopy.plot.PlotCrossSection(model=gwf, line=xc_line, ax=ax) + + # Plot HFB + xc.plot_grid() + hfb_result = xc.plot_bc("HFB", alpha=0.5) + + assert hfb_result is not None, f"HFB plot should return a result ({grid_type})" + assert isinstance(hfb_result, PatchCollection) + assert len(ax.collections) > 0 + + # plt.show() + plt.close(fig) + + +@pytest.fixture(params=["dis", "disv"]) +def vertical_hfb_xc_model(request): + """Create a MODFLOW 6 model with vertical HFBs for cross section testing. + + Vertical HFBs are barriers between vertically stacked cells (different layers). + + Parameters + ---------- + request.param : str + Grid type: "dis" (structured) or "disv" (vertex) + + Returns + ------- + tuple + (gwf_model, cross_section_line, grid_type) + """ + import numpy as np + + from flopy.utils.gridutil import get_disv_kwargs + + grid_type = request.param + + # Create simulation + sim = flopy.mf6.MFSimulation(sim_name=f"test_vhfb_xc_{grid_type}") + tdis = flopy.mf6.ModflowTdis(sim) + ims = flopy.mf6.ModflowIms(sim) + gwf = flopy.mf6.ModflowGwf(sim, modelname="test") + + # Create discretization based on grid type + if grid_type == "dis": + # Structured grid with 3 layers + dis = flopy.mf6.ModflowGwfdis( + gwf, + nlay=3, + nrow=10, + ncol=10, + delr=100.0, + delc=100.0, + top=100.0, + botm=[50.0, 0.0, -50.0], + ) + # Vertical HFB cellids: barriers between layers at same row/col + # Add vertical barriers along row 4 + vhfb_data = [ + [(0, 4, 3), (1, 4, 3), 1e-6], + [(0, 4, 4), (1, 4, 4), 1e-6], + [(0, 4, 5), (1, 4, 5), 1e-6], + ] + # Also add horizontal barriers for comparison + hhfb_data = [ + [(1, 4, 4), (1, 4, 5), 1e-6], + ] + hfb_data = vhfb_data + hhfb_data + xc_line = {"row": 4} + + elif grid_type == "disv": + # Vertex grid with 3 layers + disv_kwargs = get_disv_kwargs( + 3, 10, 10, 100.0, 100.0, 100.0, np.array([50.0, 0.0, -50.0]) + ) + disv = flopy.mf6.ModflowGwfdisv(gwf, **disv_kwargs) + # Vertical HFB cellids for vertex grid: (layer, cell2d_id) + # Cross section at y=450 intersects row 5 (cells 50-59), not row 4 + vhfb_data = [ + [(0, 5 * 10 + 3), (1, 5 * 10 + 3), 1e-6], + [(0, 5 * 10 + 4), (1, 5 * 10 + 4), 1e-6], + [(0, 5 * 10 + 5), (1, 5 * 10 + 5), 1e-6], + ] + # Horizontal barriers + hhfb_data = [ + [(1, 5 * 10 + 4), (1, 5 * 10 + 5), 1e-6], + ] + hfb_data = vhfb_data + hhfb_data + # For DISV, use line coordinates + xc_line = {"line": ([0, 450], [1000, 450])} + + # Add packages + ic = flopy.mf6.ModflowGwfic(gwf, strt=75.0) + npf = flopy.mf6.ModflowGwfnpf(gwf, save_flows=True) + hfb = flopy.mf6.ModflowGwfhfb(gwf, stress_period_data=hfb_data) + + return gwf, xc_line, grid_type + + +def test_cross_section_vertical_hfb(vertical_hfb_xc_model): + """Test plotting vertical HFB (barriers between vertically stacked cells). + + Vertical HFBs should be rendered as lines at layer interfaces in cross sections, + while horizontal HFBs are rendered as patches. + + Tests structured (DIS) and vertex (DISV) grids. + """ + gwf, xc_line, grid_type = vertical_hfb_xc_model + + # Create figure and cross section + fig, ax = plt.subplots(1, 1, figsize=(10, 5)) + xc = flopy.plot.PlotCrossSection(model=gwf, line=xc_line, ax=ax) + xc.plot_grid() + + # Plot HFB + hfb_result = xc.plot_bc("HFB", alpha=0.5) + + assert hfb_result is not None + + # Result should be a list containing both PatchCollection (horizontal) + # and LineCollection (vertical) + if isinstance(hfb_result, list): + # Should have both types + assert len(hfb_result) == 2, f"Expected 2 collections ({grid_type})" + has_patches = any(isinstance(c, PatchCollection) for c in hfb_result) + has_lines = any(isinstance(c, LineCollection) for c in hfb_result) + assert has_patches and has_lines + else: + # If only one type, it should be LineCollection (only vertical barriers in view) + # or PatchCollection (only horizontal barriers in view) + assert isinstance(hfb_result, (PatchCollection, LineCollection)) + + # plt.show() + plt.close(fig) diff --git a/autotest/test_plot_map_view.py b/autotest/test_plot_map_view.py index 641fa5817a..e9b5fe16f7 100644 --- a/autotest/test_plot_map_view.py +++ b/autotest/test_plot_map_view.py @@ -312,3 +312,270 @@ def test_plot_centers(): vert = tuple(vert) if vert not in xycenters: raise AssertionError("center location not properly plotted") + + +@pytest.fixture(params=["dis", "disv", "disu"]) +def hfb_model(request): + """Create a MODFLOW 6 model with HFB for testing different grid types. + + Parameters + ---------- + request.param : str + Grid type: "dis" (structured), "disv" (vertex), or "disu" (unstructured) + + Returns + ------- + tuple + (gwf_model, expected_barrier_count, grid_type) + """ + from flopy.utils.gridutil import get_disu_kwargs, get_disv_kwargs + + grid_type = request.param + + # Create simulation + sim = flopy.mf6.MFSimulation(sim_name=f"test_hfb_{grid_type}") + tdis = flopy.mf6.ModflowTdis(sim) + ims = flopy.mf6.ModflowIms(sim) + gwf = flopy.mf6.ModflowGwf(sim, modelname="test") + + # Create discretization based on grid type + if grid_type == "dis": + # Structured grid + dis = flopy.mf6.ModflowGwfdis( + gwf, + nlay=1, + nrow=10, + ncol=10, + delr=100.0, + delc=100.0, + ) + # HFB cellids for structured grid: (layer, row, col) + hfb_data = [ + [(0, 3, 4), (0, 3, 5), 1e-6], + [(0, 4, 4), (0, 4, 5), 1e-6], + [(0, 5, 4), (0, 5, 5), 1e-6], + [(0, 6, 4), (0, 6, 5), 1e-6], + [(0, 7, 4), (0, 7, 5), 1e-6], + ] + chd_spd = [[(0, 0, 0), 100.0], [(0, 9, 9), 95.0]] + expected_barriers = 5 + + elif grid_type == "disv": + # Vertex grid (regular grid converted to vertex format) + disv_kwargs = get_disv_kwargs(1, 10, 10, 100.0, 100.0, 100.0, [0.0]) + disv = flopy.mf6.ModflowGwfdisv(gwf, **disv_kwargs) + # HFB cellids for vertex grid: (layer, cell2d_id) + # For a 10x10 grid, cell2d IDs are row*ncol + col + # Barriers between cells in column 4 and 5 for rows 3-7 + hfb_data = [ + [(0, 3 * 10 + 4), (0, 3 * 10 + 5), 1e-6], + [(0, 4 * 10 + 4), (0, 4 * 10 + 5), 1e-6], + [(0, 5 * 10 + 4), (0, 5 * 10 + 5), 1e-6], + [(0, 6 * 10 + 4), (0, 6 * 10 + 5), 1e-6], + [(0, 7 * 10 + 4), (0, 7 * 10 + 5), 1e-6], + ] + chd_spd = [[(0, 0), 100.0], [(0, 99), 95.0]] + expected_barriers = 5 + + elif grid_type == "disu": + # Unstructured grid (regular grid converted to unstructured format) + disu_kwargs = get_disu_kwargs( + 1, 10, 10, 100.0, 100.0, 100.0, [0.0], return_vertices=True + ) + disu = flopy.mf6.ModflowGwfdisu(gwf, **disu_kwargs) + # HFB cellids for unstructured grid: (node,) + # For a 10x10 grid, node = layer * nrow * ncol + row * ncol + col + # Barriers between cells in column 4 and 5 for rows 3-7 + hfb_data = [ + [(3 * 10 + 4,), (3 * 10 + 5,), 1e-6], + [(4 * 10 + 4,), (4 * 10 + 5,), 1e-6], + [(5 * 10 + 4,), (5 * 10 + 5,), 1e-6], + [(6 * 10 + 4,), (6 * 10 + 5,), 1e-6], + [(7 * 10 + 4,), (7 * 10 + 5,), 1e-6], + ] + chd_spd = [[(0,), 100.0], [(99,), 95.0]] + expected_barriers = 5 + + # Add packages + ic = flopy.mf6.ModflowGwfic(gwf, strt=100.0) + npf = flopy.mf6.ModflowGwfnpf(gwf, save_flows=True) + chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=chd_spd) + hfb = flopy.mf6.ModflowGwfhfb(gwf, stress_period_data=hfb_data) + + return gwf, expected_barriers, grid_type + + +def test_plot_bc_hfb(hfb_model): + """Test plotting HFB (Horizontal Flow Barrier) boundary conditions. + + HFB packages have cellid1/cellid2 fields instead of a single cellid field, + representing barriers between pairs of cells. This test verifies that HFB + can be plotted as lines on the shared faces between cells. + + Tests structured (DIS), vertex (DISV), and unstructured (DISU) grids. + + Addresses issue #2676. + """ + gwf, expected_barriers, grid_type = hfb_model + + # Create map view and plot both CHD and HFB + fig = plt.figure(figsize=(8, 8)) + ax = fig.add_subplot(1, 1, 1, aspect="equal") + mapview = flopy.plot.PlotMapView(model=gwf, ax=ax) + + # Plot grid lines + grid_result = mapview.plot_grid() + assert isinstance(grid_result, LineCollection), ( + f"Expected LineCollection for grid ({grid_type})" + ) + + # Plot CHD (should work as before) + chd_result = mapview.plot_bc("CHD") + assert chd_result is not None, f"CHD plot should return a result ({grid_type})" + + # Plot HFB (the new functionality) - uses default orange color + hfb_result = mapview.plot_bc("HFB", linewidth=3) + + assert isinstance(hfb_result, LineCollection) + + # Verify that the correct number of barrier segments were plotted + segments = hfb_result.get_segments() + assert len(segments) == expected_barriers + + # Verify each segment has 2 points (start and end of the barrier line) + for seg in segments: + assert len(seg) == 2 + + # plt.show() + plt.close(fig) + + +@pytest.fixture(params=["dis", "disv"]) +def vertical_hfb_model(request): + """Create a MODFLOW 6 model with vertical HFBs for testing. + + Vertical HFBs are barriers between vertically stacked cells (different layers). + + Returns + ------- + tuple + (gwf_model, expected_barrier_count, grid_type) + """ + from flopy.utils.gridutil import get_disv_kwargs + + grid_type = request.param + + # Create simulation + sim = flopy.mf6.MFSimulation(sim_name=f"test_vhfb_{grid_type}") + tdis = flopy.mf6.ModflowTdis(sim) + ims = flopy.mf6.ModflowIms(sim) + gwf = flopy.mf6.ModflowGwf(sim, modelname="test") + + # Create discretization based on grid type + if grid_type == "dis": + # Structured grid with 3 layers + dis = flopy.mf6.ModflowGwfdis( + gwf, + nlay=3, + nrow=10, + ncol=10, + delr=100.0, + delc=100.0, + top=100.0, + botm=[50.0, 0.0, -50.0], + ) + # Vertical HFB cellids: barriers between layers at same row/col + # Add some vertical barriers between layer 0 and 1 + vhfb_data = [ + [(0, 3, 4), (1, 3, 4), 1e-6], + [(0, 4, 4), (1, 4, 4), 1e-6], + [(0, 5, 4), (1, 5, 4), 1e-6], + ] + # Also add a few horizontal barriers for comparison + hhfb_data = [ + [(1, 3, 4), (1, 3, 5), 1e-6], + [(1, 4, 4), (1, 4, 5), 1e-6], + ] + hfb_data = vhfb_data + hhfb_data + chd_spd = [[(0, 0, 0), 100.0], [(2, 9, 9), 95.0]] + expected_vertical_barriers = 3 + expected_horizontal_barriers = 2 + + elif grid_type == "disv": + # Vertex grid with 3 layers + import numpy as np + + disv_kwargs = get_disv_kwargs( + 3, 10, 10, 100.0, 100.0, 100.0, np.array([50.0, 0.0, -50.0]) + ) + disv = flopy.mf6.ModflowGwfdisv(gwf, **disv_kwargs) + # Vertical HFB cellids for vertex grid: (layer, cell2d_id) + vhfb_data = [ + [(0, 3 * 10 + 4), (1, 3 * 10 + 4), 1e-6], + [(0, 4 * 10 + 4), (1, 4 * 10 + 4), 1e-6], + [(0, 5 * 10 + 4), (1, 5 * 10 + 4), 1e-6], + ] + # Horizontal barriers + hhfb_data = [ + [(1, 3 * 10 + 4), (1, 3 * 10 + 5), 1e-6], + [(1, 4 * 10 + 4), (1, 4 * 10 + 5), 1e-6], + ] + hfb_data = vhfb_data + hhfb_data + chd_spd = [[(0, 0), 100.0], [(2, 99), 95.0]] + expected_vertical_barriers = 3 + expected_horizontal_barriers = 2 + + # Add packages + ic = flopy.mf6.ModflowGwfic(gwf, strt=100.0) + npf = flopy.mf6.ModflowGwfnpf(gwf, save_flows=True) + chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=chd_spd) + hfb = flopy.mf6.ModflowGwfhfb(gwf, stress_period_data=hfb_data) + + return gwf, expected_vertical_barriers, expected_horizontal_barriers, grid_type + + +def test_plot_bc_vertical_hfb(vertical_hfb_model): + """Test plotting vertical HFB (barriers between vertically stacked cells). + + Vertical HFBs should be rendered as patches (full cells) in map view, + while horizontal HFBs are rendered as lines. + + Tests structured (DIS) and vertex (DISV) grids. + """ + gwf, expected_vbarriers, expected_hbarriers, grid_type = vertical_hfb_model + + # Test on layer 0 - should show vertical barriers as patches + fig = plt.figure(figsize=(8, 8)) + ax = fig.add_subplot(1, 1, 1, aspect="equal") + mapview = flopy.plot.PlotMapView(model=gwf, layer=0, ax=ax) + mapview.plot_grid() + + # Plot HFB + hfb_result = mapview.plot_bc("HFB", linewidth=3) + + # Result should be a list containing both LineCollection (horizontal) + # and PatchCollection (vertical) + assert hfb_result is not None + + # If both types exist, result is a list + if isinstance(hfb_result, list): + # Should have both line and patch collections + assert len(hfb_result) == 2, f"Expected 2 collections ({grid_type})" + has_lines = any(isinstance(c, LineCollection) for c in hfb_result) + has_patches = any(isinstance(c, PatchCollection) for c in hfb_result) + assert has_lines and has_patches + + # plt.show() + plt.close(fig) + + # Test on layer 1 - should show both vertical barriers and horizontal barriers + fig = plt.figure(figsize=(8, 8)) + ax = fig.add_subplot(1, 1, 1, aspect="equal") + mapview = flopy.plot.PlotMapView(model=gwf, layer=1, ax=ax) + mapview.plot_grid() + + hfb_result = mapview.plot_bc("HFB", linewidth=3) + assert hfb_result is not None + + # plt.show() + plt.close(fig) diff --git a/flopy/plot/crosssection.py b/flopy/plot/crosssection.py index 2480fa647a..fbf0e4d58a 100644 --- a/flopy/plot/crosssection.py +++ b/flopy/plot/crosssection.py @@ -11,7 +11,12 @@ from ..utils import geometry, import_optional_dependency from ..utils.geospatial_utils import GeoSpatialUtil from . import plotutil -from .plotutil import to_mp7_endpoints, to_mp7_pathlines +from .plotutil import ( + get_shared_face_3d, + is_vertical_barrier, + to_mp7_endpoints, + to_mp7_pathlines, +) warnings.simplefilter("always", PendingDeprecationWarning) @@ -776,6 +781,119 @@ def plot_grid(self, **kwargs): ax.add_collection(col) return col + def _cellid_to_node(self, cellid): + """ + Convert a cellid tuple to a node number. + + Parameters + ---------- + cellid : tuple + Cell identifier + + Returns + ------- + int + Node number + """ + if len(cellid) == 3: + # Structured grid: (layer, row, col) + layer, row, col = cellid + return layer * self.mg.nrow * self.mg.ncol + row * self.mg.ncol + col + elif len(cellid) == 2: + # Vertex grid: (layer, cell2d) + layer, cell2d = cellid + return layer * self._ncpl + cell2d + else: + # Unstructured grid: (node,) + return cellid[0] + + def _plot_vertical_hfb_lines(self, color=None, **kwargs): + """ + Plot vertical HFBs as lines at layer interfaces. + + Parameters + ---------- + color : str + Color for the lines + **kwargs : dict + Keyword arguments (linewidth, etc.) + + Returns + ------- + LineCollection or None + """ + from matplotlib.collections import LineCollection + + if ( + not hasattr(self, "_vertical_hfbs_to_plot") + or not self._vertical_hfbs_to_plot + ): + return None + + ax = kwargs.pop("ax", self.ax) + line_segments = [] + + for cellid1, cellid2 in self._vertical_hfbs_to_plot: + # Get the 2D cell identifier (row, col for DIS or cell2d for DISV) + if len(cellid1) == 3: + # Structured grid + node_2d = cellid1[1] * self.mg.ncol + cellid1[2] + elif len(cellid1) == 2: + # Vertex grid + node_2d = cellid1[1] + else: + # Unstructured - skip for now + continue + + # Check if this cell intersects the cross section + if node_2d not in self.xypts: + continue + + # Determine the layer interface elevation + # The interface is between the two layers - use the top of the lower layer + lower_layer = max(cellid1[0], cellid2[0]) + + # Get the interface elevation at this cell + if lower_layer < len(self.elev): + interface_elev = self.elev[lower_layer, node_2d] + else: + continue + + # Get the x-coordinates along the cross section for this cell + # Need to find this cell in projpts - convert both cellids to nodes + node1 = self._cellid_to_node(cellid1) + node2 = self._cellid_to_node(cellid2) + + # Get x-coordinates from either node's projection + xcoords = [] + for node in [node1, node2]: + if node in self.projpts: + poly_verts = self.projpts[node] + # Extract x-coordinates (first element of each vertex) + xs = [v[0] for v in poly_verts] + xcoords.extend(xs) + + if len(xcoords) >= 2: + # Create a line segment at the interface elevation + x_min = min(xcoords) + x_max = max(xcoords) + line_segments.append([(x_min, interface_elev), (x_max, interface_elev)]) + + # Clear the stored vertical HFBs + self._vertical_hfbs_to_plot = [] + + if not line_segments: + return None + + # Create LineCollection + if "linewidth" not in kwargs and "lw" not in kwargs: + kwargs["linewidth"] = 2 + + lc = LineCollection(line_segments, colors=color, **kwargs) + ax.add_collection(lc) + + return lc + def plot_bc(self, name=None, package=None, kper=0, color=None, head=None, **kwargs): """ Plot boundary conditions locations for a specific boundary @@ -825,6 +943,7 @@ def plot_bc(self, name=None, package=None, kper=0, color=None, head=None, **kwar p = [p] idx = np.array([]) + for pp in p: if pp.package_type in ("lak", "sfr", "maw", "uzf"): t = plotutil.advanced_package_bc_helper(pp, self.mg, kper) @@ -836,7 +955,44 @@ def plot_bc(self, name=None, package=None, kper=0, color=None, head=None, **kwar if mflist is None: return - t = np.array([list(i) for i in mflist["cellid"]], dtype=int).T + if "cellid" in mflist.dtype.names: + t = np.array([list(i) for i in mflist["cellid"]], dtype=int).T + elif ( + "cellid1" in mflist.dtype.names + and "cellid2" in mflist.dtype.names + ): + # Barrier packages (e.g., HFB) sit at interfaces between cells. + # Separate horizontal and vertical barriers: + # - Horizontal barriers: plot both affected cells as patches + # - Vertical barriers: will plot as lines at layer interface + cellids = [] + vertical_hfbs = [] + for entry in mflist: + cellid1 = tuple(entry["cellid1"]) + cellid2 = tuple(entry["cellid2"]) + + if is_vertical_barrier(self.mg, cellid1, cellid2): + # Store vertical HFBs for line plotting + vertical_hfbs.append((cellid1, cellid2)) + else: + # Horizontal barriers - plot both cells as patches + cellids.append(list(entry["cellid1"])) + cellids.append(list(entry["cellid2"])) + + # Store vertical HFBs for later processing + if not hasattr(self, "_vertical_hfbs_to_plot"): + self._vertical_hfbs_to_plot = [] + self._vertical_hfbs_to_plot.extend(vertical_hfbs) + + if cellids: + t = np.array(cellids, dtype=int).T + else: + continue + else: + raise ValueError( + f"Package {pp.package_type} has unexpected cellid fields. " + f"Available fields: {mflist.dtype.names}" + ) if len(idx) == 0: idx = np.copy(t) @@ -886,6 +1042,16 @@ def plot_bc(self, name=None, package=None, kper=0, color=None, head=None, **kwar plotarray, masked_values=[0], head=head, cmap=cmap, norm=norm, **kwargs ) + # Plot vertical HFBs as lines at layer interfaces + line_collection = self._plot_vertical_hfb_lines(color=c, **kwargs) + + # Return both patches and lines if both exist + if line_collection is not None: + if patches is not None: + return [patches, line_collection] + else: + return line_collection + return patches def plot_centers( diff --git a/flopy/plot/map.py b/flopy/plot/map.py index 8a08f708dd..0cee045341 100644 --- a/flopy/plot/map.py +++ b/flopy/plot/map.py @@ -10,7 +10,12 @@ from ..utils import geometry from . import plotutil -from .plotutil import to_mp7_endpoints, to_mp7_pathlines +from .plotutil import ( + get_shared_face, + get_shared_face_3d, + is_vertical_barrier, + to_mp7_pathlines, +) warnings.simplefilter("always", PendingDeprecationWarning) @@ -439,6 +444,120 @@ def plot_grid(self, **kwargs): ax = self._set_axes_limits(ax) return collection + def _plot_barrier_bc(self, barrier_data, color=None, name=None, **kwargs): + """ + Plot barrier-type boundary conditions (e.g., HFB) as lines or patches. + + Horizontal barriers (between horizontally adjacent cells) are plotted as + lines on shared faces. Vertical barriers (between vertically stacked cells) + are plotted as full cell patches. + + Parameters + ---------- + barrier_data : list of tuples + List of (cellid1, cellid2) tuples representing barriers between cells + color : string + matplotlib color string + name : string + Package name for color lookup + **kwargs : dictionary + keyword arguments passed to LineCollection or PatchCollection + + Returns + ------- + matplotlib.collections.LineCollection or PatchCollection + """ + ax = kwargs.pop("ax", self.ax) + + # Separate horizontal and vertical barriers + horizontal_line_segments = [] + vertical_cell_indices = [] + + for cellid1, cellid2 in barrier_data: + # Only plot barriers on the current layer (for layered grids) + # For DISU (len==1), plot all barriers since there's no layer filtering + if len(cellid1) >= 2: # DIS or DISV (has layer info) + if cellid1[0] != self.layer and cellid2[0] != self.layer: + continue + + # Check if this is a vertical barrier + if is_vertical_barrier(self.mg, cellid1, cellid2): + # For vertical barriers, plot the cells on the current layer + if len(cellid1) >= 2: + if cellid1[0] == self.layer: + if len(cellid1) == 3: + vertical_cell_indices.append( + [self.layer, cellid1[1], cellid1[2]] + ) + else: + vertical_cell_indices.append([self.layer, cellid1[1]]) + if cellid2[0] == self.layer: + if len(cellid2) == 3: + vertical_cell_indices.append( + [self.layer, cellid2[1], cellid2[2]] + ) + else: + vertical_cell_indices.append([self.layer, cellid2[1]]) + else: + # Horizontal barrier - plot as line on shared face + shared_face = get_shared_face(self.mg, cellid1, cellid2) + if shared_face is not None: + horizontal_line_segments.append(shared_face) + + # Determine color + if color is None: + key = name[:3].upper() if name else "HFB" + c = plotutil.bc_color_dict.get(key, None) + if c is None: + c = plotutil.bc_color_dict["default"] + else: + c = color + + collections = [] + + # Plot horizontal barriers as lines + if horizontal_line_segments: + lc_kwargs = dict(kwargs.items()) + if "linewidth" not in lc_kwargs and "lw" not in lc_kwargs: + lc_kwargs["linewidth"] = 2 + + lc = LineCollection(horizontal_line_segments, colors=c, **lc_kwargs) + ax.add_collection(lc) + collections.append(lc) + + # Plot vertical barriers as patches + if vertical_cell_indices: + idx = np.array(vertical_cell_indices, dtype=int).T + pc_kwargs = { + k: v for k, v in kwargs.items() if k not in ["linewidth", "lw"] + } + + # Create a plot array with 1s for cells to plot + plotarray = np.zeros(self.mg.shape, dtype=int) + if len(self.mg.shape) > 1: + plotarray[tuple(idx)] = 1 + else: + plotarray[idx] = 1 + + # Mask the plot array + plotarray = np.ma.masked_equal(plotarray, 0) + + # Set the colormap + cmap = matplotlib.colors.ListedColormap(["none", c]) + bounds = [0, 1, 2] + norm = matplotlib.colors.BoundaryNorm(bounds, cmap.N) + + # Plot using plot_array + pc = self.plot_array(plotarray, cmap=cmap, norm=norm, **pc_kwargs) + if pc is not None: + collections.append(pc) + + if not collections: + return None + + ax = self._set_axes_limits(ax) + return collections if len(collections) > 1 else collections[0] + def plot_bc( self, name=None, @@ -496,6 +615,9 @@ def plot_bc( p = [p] idx = np.array([]) + is_barrier_package = False + barrier_faces = [] + for pp in p: if pp.package_type in ("lak", "sfr", "maw", "uzf"): t = plotutil.advanced_package_bc_helper(pp, self.mg, kper) @@ -508,13 +630,37 @@ def plot_bc( return if boundname is not None: mflist = mflist[mflist["boundname"] == boundname] - t = np.array([list(i) for i in mflist["cellid"]], dtype=int).T + + if "cellid" in mflist.dtype.names: + t = np.array([list(i) for i in mflist["cellid"]], dtype=int).T + elif ( + "cellid1" in mflist.dtype.names + and "cellid2" in mflist.dtype.names + ): + # if this is a barrier-type package (e.g. HFB), + # it has cellid1, cellid2 instead of cellid. + # collect the shared faces for plotting. + is_barrier_package = True + for entry in mflist: + barrier_faces.append( + (tuple(entry["cellid1"]), tuple(entry["cellid2"])) + ) + continue + else: + raise ValueError( + f"Package {pp.package_type} has unexpected cellid fields. " + f"Available fields: {mflist.dtype.names}" + ) if len(idx) == 0: idx = np.copy(t) else: idx = np.append(idx, t, axis=1) + # Handle barrier packages differently + if is_barrier_package: + return self._plot_barrier_bc(barrier_faces, color, name, **kwargs) + else: # modflow-2005 structured and unstructured grid if p.package_type in ("uzf", "lak"): diff --git a/flopy/plot/plotutil.py b/flopy/plot/plotutil.py index 8c1bfec419..f80183556f 100644 --- a/flopy/plot/plotutil.py +++ b/flopy/plot/plotutil.py @@ -30,6 +30,7 @@ "SFR": "teal", "UZF": "peru", "LAK": "royalblue", + "HFB": "orange", } @@ -2954,3 +2955,158 @@ def to_prt_pathlines( return df else: return ret + + +def get_shared_face(mg, cellid1, cellid2) -> list | None: + """ + Get the coordinates of the shared face between two cells. + + Parameters + ---------- + mg : Grid + Model grid + cellid1 : tuple + First cell ID + cellid2 : tuple + Second cell ID + + Returns + ------- + list or None + List of two (x, y) tuples representing the shared face endpoints, + or None if cells don't share a face + """ + if cellid1 == cellid2: + raise ValueError("cellid1 and cellid2 must be different") + + try: + if len(cellid1) == 3: + # Structured grid: (layer, row, col) + verts1 = mg.get_cell_vertices(cellid1[1], cellid1[2]) + verts2 = mg.get_cell_vertices(cellid2[1], cellid2[2]) + elif len(cellid1) == 2: + # Vertex grid: (layer, cell2d_id) + verts1 = mg.get_cell_vertices(cellid1[1]) + verts2 = mg.get_cell_vertices(cellid2[1]) + else: + # Unstructured grid: (node,) + verts1 = mg.get_cell_vertices(cellid1[0]) + verts2 = mg.get_cell_vertices(cellid2[0]) + except Exception: + return None + + tol = 1e-5 + shared_verts = [] + for v1 in verts1: + for v2 in verts2: + if np.allclose(v1, v2, rtol=tol): # reasonable tolerance? + if not any(np.allclose(v1, sv, rtol=tol) for sv in shared_verts): + shared_verts.append(v1) + break + + return shared_verts if len(shared_verts) >= 2 else None + + +def get_shared_face_3d(mg, cellid1, cellid2) -> list | None: + """ + Get the 3D coordinates of the shared face between two cells. + + Parameters + ---------- + mg : Grid + Model grid + cellid1 : tuple + First cell ID (DISU node) + cellid2 : tuple + Second cell ID (DISU node) + + Returns + ------- + list or None + List of (x, y, z) tuples representing the shared face vertices, + or None if cells don't share a face + """ + if cellid1 == cellid2: + raise ValueError("cellid1 and cellid2 must be different") + + # This method is specifically for unstructured grids + if not hasattr(mg, "vertices"): + return None + + try: + # Get 3D vertices for each cell + # For DISU, vertices are stored as (node_number, x, y, z) + node1 = cellid1[0] + node2 = cellid2[0] + + # Get cell vertex indices + if hasattr(mg, "iverts"): + # iverts is a jagged array of vertex indices for each cell + verts_idx1 = mg.iverts[node1] + verts_idx2 = mg.iverts[node2] + else: + return None + + # Get 3D coordinates for vertices + verts1 = [mg.vertices[idx] for idx in verts_idx1] + verts2 = [mg.vertices[idx] for idx in verts_idx2] + + # Find shared vertices in 3D + tol = 1e-5 + shared_verts = [] + for v1 in verts1: + for v2 in verts2: + if np.allclose(v1, v2, rtol=tol): + if not any(np.allclose(v1, sv, rtol=tol) for sv in shared_verts): + shared_verts.append(v1) + break + + return shared_verts if len(shared_verts) >= 2 else None + + except Exception: + return None + + +def is_vertical_barrier(mg, cellid1, cellid2) -> bool: + """ + Determine if a barrier is vertical (between vertically stacked cells). + + Parameters + ---------- + mg : Grid + Model grid + cellid1 : tuple + First cell ID + cellid2 : tuple + Second cell ID + + Returns + ------- + bool + True if barrier is vertical (cells differ only in layer), False otherwise + """ + if len(cellid1) == 3: + # Structured grid: (layer, row, col) + # Vertical if layers differ but row and col are the same + return ( + cellid1[0] != cellid2[0] + and cellid1[1] == cellid2[1] + and cellid1[2] == cellid2[2] + ) + elif len(cellid1) == 2: + # Vertex grid: (layer, cell2d) + # Vertical if layers differ but cell2d is the same + return cellid1[0] != cellid2[0] and cellid1[1] == cellid2[1] + else: + # Unstructured grid: (node,) + # Infer from geometry: check the orientation of the shared face + # If the face is horizontal (all z-coords equal), it's a vertical barrier + # If the face is vertical (z-coords differ), it's a horizontal barrier + shared_face_3d = get_shared_face_3d(mg, cellid1, cellid2) + if shared_face_3d is None: + # No shared face found, can't determine orientation + return False + + # Check if all z-coordinates are the same (horizontal face) + z_coords = [v[2] for v in shared_face_3d] + return np.allclose(z_coords, z_coords[0], rtol=1e-5)