diff --git a/flopy/utils/geospatial_utils.py b/flopy/utils/geospatial_utils.py
index b201af9841..6d60971692 100644
--- a/flopy/utils/geospatial_utils.py
+++ b/flopy/utils/geospatial_utils.py
@@ -330,7 +330,7 @@ def __init__(self, obj, shapetype=None):
                     MultiPolygon,
                 ),
             ):
-                for geom in list(obj):
+                for geom in obj.geoms:
                     self.__collection.append(GeoSpatialUtil(geom))
 
     def __iter__(self):
diff --git a/flopy/utils/gridintersect.py b/flopy/utils/gridintersect.py
index 481e464298..843dd54c30 100644
--- a/flopy/utils/gridintersect.py
+++ b/flopy/utils/gridintersect.py
@@ -25,6 +25,33 @@
 except ImportError:
     shply = False
 
+import contextlib
+import shapely
+import warnings
+from distutils.version import LooseVersion
+
+SHAPELY_GE_20 = str(shapely.__version__) >= LooseVersion("2.0")
+
+try:
+    from shapely.errors import ShapelyDeprecationWarning as shapely_warning
+except ImportError:
+    shapely_warning = None
+
+if shapely_warning is not None and not SHAPELY_GE_20:
+
+    @contextlib.contextmanager
+    def ignore_shapely_warnings_for_object_array():
+        with warnings.catch_warnings():
+            warnings.filterwarnings("ignore", category=shapely_warning)
+            yield
+
+
+else:
+
+    @contextlib.contextmanager
+    def ignore_shapely_warnings_for_object_array():
+        yield
+
 
 def parse_shapely_ix_result(collection, ix_result, shptyps=None):
     """Recursive function for parsing shapely intersection results. Returns a
@@ -62,7 +89,7 @@ def parse_shapely_ix_result(collection, ix_result, shptyps=None):
         return collection
     # recursion for collections
     elif hasattr(ix_result, "geoms"):
-        for ishp in ix_result:
+        for ishp in ix_result.geoms:
             parse_shapely_ix_result(collection, ishp, shptyps=shptyps)
     # if collecting all types
     elif shptyps[0] is None:
@@ -465,7 +492,8 @@ def _intersect_point_shapely(self, shp, sort_by_cellid=True):
             names=["cellids", "vertices", "ixshapes"],
             formats=["O", "O", "O"],
         )
-        rec.ixshapes = isectshp
+        with ignore_shapely_warnings_for_object_array():
+            rec.ixshapes = isectshp
         rec.vertices = vertices
         rec.cellids = cellids
 
@@ -534,7 +562,8 @@ def _intersect_linestring_shapely(
             names=["cellids", "vertices", "lengths", "ixshapes"],
             formats=["O", "O", "f8", "O"],
         )
-        rec.ixshapes = isectshp
+        with ignore_shapely_warnings_for_object_array():
+            rec.ixshapes = isectshp
         rec.vertices = vertices
         rec.lengths = lengths
         rec.cellids = cellids
@@ -597,7 +626,8 @@ def _intersect_polygon_shapely(self, shp, sort_by_cellid=True):
             names=["cellids", "vertices", "areas", "ixshapes"],
             formats=["O", "O", "f8", "O"],
         )
-        rec.ixshapes = isectshp
+        with ignore_shapely_warnings_for_object_array():
+            rec.ixshapes = isectshp
         rec.vertices = vertices
         rec.areas = areas
         rec.cellids = cellids
@@ -649,10 +679,12 @@ def _intersect_point_structured(self, shp):
 
         Xe, Ye = self.mfgrid.xyedges
 
-        try:
-            iter(shp)
-        except TypeError:
+        if isinstance(shp, Point):
             shp = [shp]
+        elif isinstance(shp, MultiPoint):
+            shp = list(shp.geoms)
+        else:
+            raise ValueError("expected Point or MultiPoint")
 
         ixshapes = []
         for p in shp:
@@ -710,7 +742,8 @@ def _intersect_point_structured(self, shp):
             len(nodelist), names=["cellids", "ixshapes"], formats=["O", "O"]
         )
         rec.cellids = nodelist
-        rec.ixshapes = ixshapes
+        with ignore_shapely_warnings_for_object_array():
+            rec.ixshapes = ixshapes
         return rec
 
     def _intersect_linestring_structured(self, shp, keepzerolengths=False):
@@ -764,7 +797,7 @@ def _intersect_linestring_structured(self, shp, keepzerolengths=False):
         if lineclip.geom_type == "MultiLineString":  # there are multiple lines
             nodelist, lengths, vertices = [], [], []
             ixshapes = []
-            for ls in lineclip:
+            for ls in lineclip.geoms:
                 n, l, v, ixs = self._get_nodes_intersecting_linestring(ls)
                 nodelist += n
                 lengths += l
@@ -888,7 +921,8 @@ def _intersect_linestring_structured(self, shp, keepzerolengths=False):
         rec.vertices = vertices
         rec.lengths = lengths
         rec.cellids = nodelist
-        rec.ixshapes = ixshapes
+        with ignore_shapely_warnings_for_object_array():
+            rec.ixshapes = ixshapes
 
         return rec
 
@@ -1320,7 +1354,8 @@ def _intersect_polygon_structured(self, shp):
         rec.vertices = vertices
         rec.areas = areas
         rec.cellids = nodelist
-        rec.ixshapes = ixshapes
+        with ignore_shapely_warnings_for_object_array():
+            rec.ixshapes = ixshapes
 
         return rec