diff --git a/autotest/test_binaryfile.py b/autotest/test_binaryfile.py
index e3308ef9c..101f6a1cc 100644
--- a/autotest/test_binaryfile.py
+++ b/autotest/test_binaryfile.py
@@ -10,7 +10,7 @@
 import pytest
 from matplotlib import pyplot as plt
 from matplotlib.axes import Axes
-from modflow_devtools.markers import requires_exe
+from modflow_devtools.markers import requires_exe, requires_pkg
 
 import flopy
 from flopy.utils import (
@@ -679,6 +679,7 @@ def dis_sim(function_tmpdir):
 
 
 @pytest.mark.requires_exe("mf6")
+@requires_pkg("shapely")
 def test_headfile_get_ts_disv_grid(dis_sim, function_tmpdir):
     """Test HeadFile.get_ts() with DISV grid using both new and old index formats."""
     from flopy.mf6 import ModflowGwfchd, ModflowGwfdisv
@@ -744,6 +745,7 @@ def test_headfile_get_ts_disv_grid(dis_sim, function_tmpdir):
 
 
 @pytest.mark.requires_exe("mf6")
+@requires_pkg("shapely")
 def test_headfile_get_ts_disu_grid(dis_sim, function_tmpdir):
     """Test HeadFile.get_ts() with DISU grid using both new and old index formats."""
     from flopy.mf6 import ModflowGwfchd, ModflowGwfdisu
diff --git a/autotest/test_cellbudgetfile.py b/autotest/test_cellbudgetfile.py
index f99382be2..e01d598f0 100644
--- a/autotest/test_cellbudgetfile.py
+++ b/autotest/test_cellbudgetfile.py
@@ -3,6 +3,7 @@
 import numpy as np
 import pandas as pd
 import pytest
+from modflow_devtools.markers import requires_pkg
 
 from autotest.conftest import get_example_data_path
 from flopy.mf6.modflow.mfsimulation import MFSimulation
@@ -861,6 +862,7 @@ def test_cellbudgetfile_get_ts_aux_vars_mf6_dis(dis_sim):
 
 
 @pytest.mark.requires_exe("mf6")
+@requires_pkg("shapely")
 def test_cellbudgetfile_get_ts_aux_vars_mf6_disv(dis_sim):
     from flopy.mf6 import ModflowGwfchd, ModflowGwfdisv
 
@@ -938,6 +940,7 @@ def test_cellbudgetfile_get_ts_aux_vars_mf6_disv(dis_sim):
 
 
 @pytest.mark.requires_exe("mf6")
+@requires_pkg("shapely")
 def test_cellbudgetfile_get_ts_aux_vars_mf6_disu(dis_sim):
     from flopy.mf6 import ModflowGwfchd, ModflowGwfdisu
 
@@ -1013,6 +1016,7 @@ def test_cellbudgetfile_get_ts_aux_vars_mf6_disu(dis_sim):
 
 
 @pytest.mark.requires_exe("mf6")
+@requires_pkg("shapely")
 def test_cellbudgetfile_get_ts_imeth1_disv_grid(dis_sim, function_tmpdir):
     """Test that IMETH=1 budget terms (like STO-SS) work with DISV grids.
 
@@ -1074,6 +1078,7 @@ def test_cellbudgetfile_get_ts_imeth1_disv_grid(dis_sim, function_tmpdir):
 
 
 @pytest.mark.requires_exe("mf6")
+@requires_pkg("shapely")
 def test_cellbudgetfile_get_ts_backwards_compatible_idx_format(
     dis_sim, function_tmpdir
 ):
diff --git a/autotest/test_export.py b/autotest/test_export.py
index 0c107cb2f..1bd9a9514 100644
--- a/autotest/test_export.py
+++ b/autotest/test_export.py
@@ -7,7 +7,6 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import pytest
-import shapefile
 from flaky import flaky
 from modflow_devtools.markers import excludes_platform, requires_exe, requires_pkg
 from modflow_devtools.misc import has_pkg
@@ -2132,6 +2131,8 @@ def test_to_shapefile_raises_attributeerror():
 @pytest.mark.parametrize("use_pandas", [True])  # TODO: test non-pandas
 @pytest.mark.parametrize("sparse", [True, False])
 def test_mf6_chd_shapefile_export_structured(function_tmpdir, use_pandas, sparse):
+    import shapefile
+
     from flopy.mf6 import (
         MFSimulation,
         ModflowGwf,
@@ -2183,6 +2184,8 @@ def test_mf6_chd_shapefile_export_structured(function_tmpdir, use_pandas, sparse
 @pytest.mark.parametrize("sparse", [True])  # TODO: test non-sparse
 def test_mf6_chd_shapefile_export_unstructured(function_tmpdir, use_pandas, sparse):
     """Test CHD package shapefile export for DISU (unstructured) grids"""
+    import shapefile
+
     from flopy.mf6 import (
         MFSimulation,
         ModflowGwf,
diff --git a/autotest/test_grid.py b/autotest/test_grid.py
index 9ff8a1608..704324da7 100644
--- a/autotest/test_grid.py
+++ b/autotest/test_grid.py
@@ -12,7 +12,6 @@
 from matplotlib import pyplot as plt
 from modflow_devtools.markers import requires_exe, requires_pkg
 from modflow_devtools.misc import has_pkg
-from scipy.spatial import Delaunay
 
 from autotest.test_dis_cases import case_dis, case_disv
 from autotest.test_grid_cases import GridCases
@@ -223,6 +222,7 @@ def test_structured_grid_get_cell_vertices():
     assert v2 == v5, "Named i,j should match"
 
 
+@requires_pkg("shapely")
 def test_vertex_grid_get_cell_vertices():
     """Test VertexGrid.get_cell_vertices() with various input forms"""
     disv_props = get_disv_kwargs(2, 10, 10, 10.0, 10.0, 100.0, [50.0, 0.0])
@@ -250,6 +250,7 @@ def test_vertex_grid_get_cell_vertices():
     assert v4 == v5, "Node and (layer, cell2d) should match"
 
 
+@requires_pkg("shapely")
 def test_unstructured_grid_get_cell_vertices():
     """Test UnstructuredGrid.get_cell_vertices() with various input forms"""
     disu_props = get_disu_kwargs(
@@ -548,8 +549,11 @@ def test_structured_grid_intersect_array(simple_structured_grid):
     assert np.isnan(rows_mixed[1])  # Second point out of bounds
 
 
+@requires_pkg("scipy")
 def test_vertex_grid_intersect_array():
     """Test VertexGrid.intersect() with array inputs."""
+    from scipy.spatial import Delaunay
+
     # Create a simple vertex grid using Delaunay triangulation
     np.random.seed(42)
     n_points = 50
@@ -589,8 +593,11 @@ def test_vertex_grid_intersect_array():
     assert np.isnan(results_mixed[1])  # Second point out of bounds
 
 
+@requires_pkg("scipy")
 def test_unstructured_grid_intersect_array():
     """Test UnstructuredGrid.intersect() with array inputs."""
+    from scipy.spatial import Delaunay
+
     # Create a simple unstructured grid using Delaunay triangulation
     np.random.seed(42)
     n_points = 50
@@ -1477,9 +1484,7 @@ def test_voronoi_vertex_grid(function_tmpdir):
         GridCases.voronoi_nested_circles,
         GridCases.voronoi_polygons,
         GridCases.voronoi_many_polygons,
-    ]
-    if (has_pkg("shapely", True) and has_pkg("scipy", True))
-    else [],
+    ],
     ids=[
         "voronoi_polygon",
         "voronoi_rectangle",
diff --git a/autotest/test_gridintersect.py b/autotest/test_gridintersect.py
index 2a1a876fd..e925929fa 100644
--- a/autotest/test_gridintersect.py
+++ b/autotest/test_gridintersect.py
@@ -1804,6 +1804,7 @@ def test_tri_grid_intersect_multiple_polygon_array(rtree):
         ix.intersect(p, geo_dataframe=df_toggle)
 
 
+@requires_pkg("shapely")
 def test_rtree_false_raises_in_points_to_cellids():
     """rtree=False raises error in points_to_cellids."""
     gr = get_rect_grid()
@@ -1816,6 +1817,7 @@ def test_rtree_false_raises_in_points_to_cellids():
         ix.points_to_cellids(pts)
 
 
+@requires_pkg("shapely")
 def test_rtree_false_raises_with_arrays_in_intersects():
     """rtree=False raises error in points_to_cellids."""
     gr = get_rect_grid()
diff --git a/autotest/test_hfb_util.py b/autotest/test_hfb_util.py
index 31f9958df..9f2a194b2 100644
--- a/autotest/test_hfb_util.py
+++ b/autotest/test_hfb_util.py
@@ -1,11 +1,14 @@
 import numpy as np
-from modflow_devtools.markers import requires_exe
+from modflow_devtools.markers import requires_exe, requires_pkg
 
 import flopy
 from flopy.utils.hfb_util import make_hfb_array
 from flopy.utils.triangle import Triangle
 from flopy.utils.voronoi import VoronoiGrid
 
+# every test builds the hfb array, which is done with shapely
+pytestmark = requires_pkg("shapely")
+
 
 def structured_sim():
     lx = 100
diff --git a/autotest/test_lgrutil.py b/autotest/test_lgrutil.py
index 9ff50c419..4c5b35673 100644
--- a/autotest/test_lgrutil.py
+++ b/autotest/test_lgrutil.py
@@ -1,5 +1,6 @@
 import numpy as np
 import pytest
+from modflow_devtools.markers import requires_pkg
 
 from flopy.discretization import StructuredGrid
 from flopy.utils.cvfdutil import get_disv_gridprops, gridlist_to_verts
@@ -148,6 +149,7 @@ def test_lgr_variable_rc_spacing():
     assert np.allclose(lgr.delc, answer), f"{lgr.delc} /= {answer}"
 
 
+@requires_pkg("shapely")
 def test_lgr_hanging_vertices():
     # Define parent grid information
     xoffp = 0.0
diff --git a/autotest/test_plot_cross_section.py b/autotest/test_plot_cross_section.py
index 460babdf0..e5829b45f 100644
--- a/autotest/test_plot_cross_section.py
+++ b/autotest/test_plot_cross_section.py
@@ -298,6 +298,9 @@ def hfb_xc_model(request):
     from flopy.utils.gridutil import get_disu_kwargs, get_disv_kwargs
 
     grid_type = request.param
+    if grid_type != "dis":
+        # the vertex and unstructured grids are built with shapely
+        pytest.importorskip("shapely")
 
     # Create simulation
     sim = flopy.mf6.MFSimulation(sim_name=f"test_hfb_xc_{grid_type}")
@@ -414,6 +417,9 @@ def vertical_hfb_xc_model(request):
     from flopy.utils.gridutil import get_disv_kwargs
 
     grid_type = request.param
+    if grid_type != "dis":
+        # the vertex and unstructured grids are built with shapely
+        pytest.importorskip("shapely")
 
     # Create simulation
     sim = flopy.mf6.MFSimulation(sim_name=f"test_vhfb_xc_{grid_type}")
diff --git a/autotest/test_plot_map_view.py b/autotest/test_plot_map_view.py
index 7eed743ca..d42ff156d 100644
--- a/autotest/test_plot_map_view.py
+++ b/autotest/test_plot_map_view.py
@@ -421,6 +421,9 @@ def hfb_model(request):
     from flopy.utils.gridutil import get_disu_kwargs, get_disv_kwargs
 
     grid_type = request.param
+    if grid_type != "dis":
+        # the vertex and unstructured grids are built with shapely
+        pytest.importorskip("shapely")
 
     # Create simulation
     sim = flopy.mf6.MFSimulation(sim_name=f"test_hfb_{grid_type}")
@@ -554,6 +557,9 @@ def vertical_hfb_model(request):
     from flopy.utils.gridutil import get_disv_kwargs
 
     grid_type = request.param
+    if grid_type != "dis":
+        # the vertex and unstructured grids are built with shapely
+        pytest.importorskip("shapely")
 
     # Create simulation
     sim = flopy.mf6.MFSimulation(sim_name=f"test_vhfb_{grid_type}")
diff --git a/autotest/test_postprocessing.py b/autotest/test_postprocessing.py
index d3fbbd2f6..fcebe5fbf 100644
--- a/autotest/test_postprocessing.py
+++ b/autotest/test_postprocessing.py
@@ -2,7 +2,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import pytest
-from modflow_devtools.markers import requires_exe
+from modflow_devtools.markers import requires_exe, requires_pkg
 
 import flopy
 from flopy.mf6 import (
@@ -556,6 +556,7 @@ def test_get_transmissivities_mf6_structured(function_tmpdir):
     assert np.array_equal(Tcoords, Tcellids)
 
 
+@requires_pkg("shapely")
 def test_get_transmissivities_mf6_vertex(function_tmpdir):
     nl = 1
     nr = 8
@@ -593,6 +594,7 @@ def test_get_transmissivities_mf6_vertex(function_tmpdir):
         assert "r, c parameters only valid for structured grids" in str(e.value)
 
 
+@requires_pkg("shapely")
 def test_get_transmissivities_mf6_unstructured(function_tmpdir):
     nl = 1
     nr = 8