fix(gridintersect): fix multiple issues (#2343)

Fix #2342:

    * fix typo in shapely.multilinestrings
    * fix issue with np.apply_along_axis
    * consider geometry collections when computing overlaps
    * add test for vertex mode
    * add inactive test for structured mode (not yet working)

Notes:

    * np.apply_along_axis was reducing result from multiple GeometryCollections to a single MultiLineString causing duplication of shapes in the intersection result.
    * structured support will be dropped in future releases (see #2344)

autotest/test_gridintersect.py

@@ -775,6 +775,26 @@ def test_tri_grid_linestring_cell_boundary_return_all_ix_shapely(rtree):
     assert len(r) == 3
 
 
+@requires_pkg("shapely")
+def test_rect_vertex_grid_linestring_geomcollection():
+    gr = get_rect_vertex_grid()
+    ix = GridIntersect(gr, method="vertex")
+    ls = LineString(
+        [
+            (20.0, 0.0),
+            (5.0, 5.0),
+            (15.0, 7.5),
+            (10.0, 10.0),
+            (5.0, 15.0),
+            (10.0, 19.0),
+            (10.0, 20.0),
+        ]
+    )
+    result = ix.intersect(ls)
+    assert len(result) == 3
+    assert np.allclose(result.lengths.sum(), ls.length)
+
+
 # %% test polygon structured
 
 

flopy/utils/gridintersect.py

@@ -547,7 +547,7 @@ def parse_linestrings_in_geom_collection(gc):
                 parts = shapely.get_parts(gc)
                 parts = parts[np.isin(shapely.get_type_id(parts), line_ids)]
                 if len(parts) > 1:
-                    p = shapely.multilinestring(parts)
+                    p = shapely.multilinestrings(parts)
                 elif len(parts) == 0:
                     p = shapely.LineString()
                 else:
@@ -558,19 +558,25 @@ def parse_linestrings_in_geom_collection(gc):
                 geomtype_ids[~mask_empty & mask_type]
                 == shapely.GeometryType.GEOMETRYCOLLECTION
             )
-            ixresult[mask_gc] = np.apply_along_axis(
-                parse_linestrings_in_geom_collection,
-                axis=0,
-                arr=ixresult[mask_gc],
-            )
+            # NOTE: not working for multiple geometry collections, result is reduced
+            # to a single multilinestring, which causes doubles in the result
+            # ixresult[mask_gc] = np.apply_along_axis(
+            #     parse_linestrings_in_geom_collection,
+            #     axis=0,
+            #     arr=ixresult[mask_gc],
+            # )
+            ixresult[mask_gc] = [
+                parse_linestrings_in_geom_collection(gc)
+                for gc in ixresult[mask_gc]
+            ]
 
         if not return_all_intersections:
             # intersection with grid cell boundaries
             ixbounds = shapely.intersection(
                 shp, shapely.get_exterior_ring(self.geoms[qcellids])
             )
             mask_bnds_empty = shapely.is_empty(ixbounds)
-            mask_bnds_type = np.isin(shapely.get_type_id(ixbounds), line_ids)
+            mask_bnds_type = np.isin(shapely.get_type_id(ixbounds), all_ids)
             # get ids of boundary intersections
             idxs = np.nonzero(~mask_bnds_empty & mask_bnds_type)[0]
 
@@ -586,7 +592,7 @@ def parse_linestrings_in_geom_collection(gc):
                 mask_bnds_empty = shapely.is_empty(
                     isect
                 )  # select boundary ix result
-                mask_overlap = np.isin(shapely.get_type_id(isect), line_ids)
+                mask_overlap = np.isin(shapely.get_type_id(isect), all_ids)
 
                 # calculate difference between self and overlapping result
                 diff = shapely.difference(