From 4f354dda05a5048786df3577d8fa82d9ac7b416e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dav=C3=ADd=20Brakenhoff?= <d.brakenhoff@artesia-water.nl>
Date: Tue, 22 Oct 2024 16:09:39 +0200
Subject: [PATCH] fix bug: typo shapely function call within gridintersect
 #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)

---
 autotest/test_gridintersect.py | 53 +++++++++++++++++++++++++++-------
 flopy/utils/gridintersect.py   | 23 ++++++++++-----
 2 files changed, 57 insertions(+), 19 deletions(-)

diff --git a/autotest/test_gridintersect.py b/autotest/test_gridintersect.py
index 268de09f90..4773dbe4da 100644
--- a/autotest/test_gridintersect.py
+++ b/autotest/test_gridintersect.py
@@ -484,6 +484,28 @@ def test_rect_grid_linestring_starting_on_vertex():
     assert result.cellids[0] == (1, 1)
 
 
+# @requires_pkg("shapely")
+# def test_rect_grid_linestring_geomcolletion():
+#     gr = get_rect_grid()
+#     ix = GridIntersect(gr, method="structured")
+#     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  # TODO: currently returns 4
+#     assert np.allclose(
+#         result.lengths.sum(), ls.length
+#     )  # TODO: length is 1m longer than ls
+
+
 # %% test linestring shapely
 
 
@@ -737,6 +759,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
 
 
@@ -1555,14 +1597,3 @@ def test_create_raster_from_array_transform(example_data_path):
         or not ((ymax - ymin) / (rymax - rymin)) == 2
     ):
         raise AssertionError("Transform based raster not working properly")
-
-
-if __name__ == "__main__":
-    sgr = get_rect_grid(angrot=45.0, xyoffset=10.0)
-    ls = LineString([(5, 10.0 + np.sqrt(200.0)), (15, 10.0 + np.sqrt(200.0))])
-    ix = GridIntersect(sgr, method="structured")
-    result = ix.intersect(ls)
-    assert len(result) == 2
-    ix = GridIntersect(sgr, method="structured", local=True)
-    result = ix.intersect(ls)
-    assert len(result) == 0
diff --git a/flopy/utils/gridintersect.py b/flopy/utils/gridintersect.py
index eed25d3102..b54631b460 100644
--- a/flopy/utils/gridintersect.py
+++ b/flopy/utils/gridintersect.py
@@ -966,7 +966,7 @@ def parse_linestrings_in_geom_collection(gc):
                 parts = shapely.get_parts(gc)
                 parts = parts[np.isin(shapely.get_type_id(parts), [1, 5])]
                 if len(parts) > 1:
-                    p = shapely.multilinestring(parts)
+                    p = shapely.multilinestrings(parts)
                 elif len(parts) == 0:
                     p = shapely.LineString()
                 else:
@@ -974,11 +974,17 @@ def parse_linestrings_in_geom_collection(gc):
                 return p
 
             mask_gc = geomtype_ids[~mask_empty & mask_type] == 7
-            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
@@ -986,7 +992,8 @@ def parse_linestrings_in_geom_collection(gc):
                 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), [1, 5])
+            mask_bnds_type = np.isin(shapely.get_type_id(ixbounds), [1, 5, 7])
+
             # get ids of boundary intersections
             idxs = np.nonzero(~mask_bnds_empty & mask_bnds_type)[0]
 
@@ -1002,7 +1009,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), [1, 5])
+                mask_overlap = np.isin(shapely.get_type_id(isect), [1, 5, 7])
 
                 # calculate difference between self and overlapping result
                 diff = shapely.difference(