diff --git a/autotest/t065_test_gridintersect.py b/autotest/t065_test_gridintersect.py
index 1ec78ed68..c1bc084e4 100644
--- a/autotest/t065_test_gridintersect.py
+++ b/autotest/t065_test_gridintersect.py
@@ -205,7 +205,8 @@ def test_rect_grid_point_outside():
         return
     gr = get_rect_grid()
     ix = GridIntersect(gr, method="structured")
-    result = ix.intersect(Point(25.0, 25.0))
+    # use GeoSpatialUtil to convert to shapely geometry
+    result = ix.intersect((25.0, 25.0), shapetype="point")
     assert len(result) == 0
     return result
 
@@ -1289,9 +1290,9 @@ def test_all_intersections_shapely_no_strtree():
     # offset and rotated grids
     test_point_offset_rot_structured_grid_shapely(rtree=False)
     test_linestring_offset_rot_structured_grid_shapely(rtree=False)
-    ix = test_polygon_offset_rot_structured_grid_shapely(rtree=False)
+    test_polygon_offset_rot_structured_grid_shapely(rtree=False)
 
-    return ix
+    return
 
 
 # %% test rasters
@@ -1416,8 +1417,3 @@ def test_raster_sampling_methods():
             raise AssertionError(
                 f"{method} resampling returning incorrect values"
             )
-
-
-if __name__ == "__main__":
-    test_rasters()
-    test_raster_sampling_methods()
\ No newline at end of file
diff --git a/flopy/utils/gridintersect.py b/flopy/utils/gridintersect.py
index 6e3d9141a..18d418314 100644
--- a/flopy/utils/gridintersect.py
+++ b/flopy/utils/gridintersect.py
@@ -186,29 +186,34 @@ def __init__(self, mfgrid, method=None, rtree=True):
                 )
             )
 
-    def intersect(self, shp, **kwargs):
-        """
-        Method to intersect a shape with a model grid
+    def intersect(
+        self, shp, shapetype=None, sort_by_cellid=True, keepzerolengths=False
+    ):
+        """Method to intersect a shape with a model grid.
 
         Parameters
         ----------
         shp : shapely.geometry, geojson object, shapefile.Shape,
-              or flopy geomerty object
+              or flopy geometry object
+        shapetype : str, optional
+            type of shape (i.e. "point", "linestring", "polygon" or
+            their multi-variants), used by GeoSpatialUtil if shp is
+            passed as a list of vertices, default is None
         sort_by_cellid : bool
-            Sort results by cellid
+            sort results by cellid, ensures cell with lowest cellid is
+            returned for boundary cases when using vertex methods, default
+            is True
         keepzerolengths : bool
             boolean method to keep zero length intersections for
-            linestring intersection
+            linestring intersection, used when shp is of type "linestring"
 
         Returns
         -------
         numpy.recarray
             a record array containing information about the intersection
         """
-        gu = GeoSpatialUtil(shp)
+        gu = GeoSpatialUtil(shp, shapetype=shapetype)
         shp = gu.shapely
-        sort_by_cellid = kwargs.pop("sort_by_cellid", True)
-        keepzerolengths = kwargs.pop("keepzerolengths", False)
 
         if gu.shapetype in ("Point", "MultiPoint"):
             if (
@@ -654,7 +659,7 @@ def _intersect_polygon_shapely(self, shp, sort_by_cellid=True):
 
         return rec
 
-    def intersects(self, shp):
+    def intersects(self, shp, shapetype=None):
         """Return cellIDs for shapes that intersect with shape.
 
         Parameters
@@ -662,7 +667,10 @@ def intersects(self, shp):
         shp : shapely.geometry, geojson geometry, shapefile.shape,
               or flopy geometry object
             shape to intersect with the grid
-
+        shapetype : str, optional
+            type of shape (i.e. "point", "linestring", "polygon" or
+            their multi-variants), used by GeoSpatialUtil if shp is
+            passed as a list of vertices, default is None
         Returns
         -------
         rec : numpy.recarray
@@ -670,7 +678,7 @@ def intersects(self, shp):
             the shape intersects with
         """
         # query grid
-        shp = GeoSpatialUtil(shp).shapely
+        shp = GeoSpatialUtil(shp, shapetype=shapetype).shapely
 
         qresult = self.query_grid(shp)
         # filter result further if possible (only strtree and filter methods)
@@ -754,7 +762,6 @@ def _intersect_point_structured(self, shp):
                 tempnodes.append(node)
                 tempshapes.append(ixs)
             else:
-                # TODO: not sure if this is correct
                 tempshapes[-1] = shapely_geo.MultiPoint([tempshapes[-1], ixs])
 
         ixshapes = tempshapes