feat(Gridintersect): add shapetype kwarg (#1301)

* refactor(make-release): Software citation developed from version.py

* release: update mf6 v6.1.1 dfn and flopy mf6 classes (#1)

* docs(misc): clean up docstrings

* refactor(GridIntersect): fix error message
for unsupported methods

* bug(GridIntersect): fix #10
ensure one entry per grid cell in polygon intersection result
add test for this case

* bug(gridintersect): fixes #917

bug in structured mode where linestrings were not handled
correctly if they pass in and out of grid cell

* fix(GridIntersect): #1035
- fix vertices when intersecting linestring with offset grids

* feat(GridIntersect): add shapetype kwarg
allow intersection shape to be passed as list of vertices (modify test to test this)
specify keyword arguments explicitly in intersect calls
formatting

Co-authored-by: jdhughes-usgs <jdhughes@usgs.gov>
Co-authored-by: Christian Langevin <langevin@usgs.gov>

Closes #727

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()

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,23 +659,26 @@ 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
         ----------
         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
             a record array containing cell IDs of the gridcells
             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