Repository: natcap/invest (python)

Issue:
CV: TypeError: 'MultiLineString' object is not iterable
This is a somewhat rare edge case that comes up with complicated landmasses & AOIs. It's hard enough to contrive an example that we don't test for it, apparently. Our code to handle `MultLineString` iteration simply is not compatible with Shapely 2.0. https://shapely.readthedocs.io/en/stable/migration.html#multi-part-geometries-will-no-longer-be-sequences-length-iterable-indexable

I have some user's data so I can reproduce and fix it.

```python
Traceback (most recent call last):
  File "taskgraph\Task.py", line 674, in add_task
  File "taskgraph\Task.py", line 1093, in _call
  File "natcap\invest\coastal_vulnerability.py", line 1176, in prepare_landmass_line_index_and_interpolate_shore_points
TypeError: 'MultiLineString' object is not iterable
10/21/2024 17:29:44  natcap.invest.utils ERROR    Exception while executing natcap.invest.coastal_vulnerability
Traceback (most recent call last):
  File "natcap\invest\utils.py", line 165, in prepare_workspace
  File "invest\cli.py", line 470, in main
  File "natcap\invest\coastal_vulnerability.py", line 856, in execute
  File "taskgraph\Task.py", line 674, in add_task
  File "taskgraph\Task.py", line 1093, in _call
  File "natcap\invest\coastal_vulnerability.py", line 1176, in prepare_landmass_line_index_and_interpolate_shore_points
TypeError: 'MultiLineString' object is not iterable
```

Pull request description:
CV bugfix: handling multilinestrings
Shapely 2.0 had some breaking changes, including how to iterate over multipart geometries. We missed one case in this model when we updated to shapely 2.0.

We already had a helper function to help with these cases. And I modified a test to cover it.

Fixes #1657 -- see for more details.

I also alphabetized the HISTORY, after leading with "General" and "Workbench" sections.

## Checklist
- [x] Updated HISTORY.rst and link to any relevant issue (if these changes are user-facing)
- [ ] Updated the user's guide (if needed)
- [ ] Tested the Workbench UI (if relevant)


Patch:
diff --git a/HISTORY.rst b/HISTORY.rst
index a8bcf0514..71214591b 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -44,13 +44,6 @@ Unreleased Changes
       ``pygeoprocessing.get_raster_info`` and
       ``pygeoprocessing.get_vector_info``.
       https://github.com/natcap/invest/issues/1645
-* Forest Carbon Edge Effects
-    * Updating vector reprojection to allow partial reprojection.  Related to
-      https://github.com/natcap/invest/issues/1645
-* Urban Nature Access
-    * The model now works as expected when the user provides an LULC raster
-      that does not have a nodata value defined.
-      https://github.com/natcap/invest/issues/1293
 * Workbench
     * Several small updates to the model input form UI to improve usability
       and visual consistency (https://github.com/natcap/invest/issues/912).
@@ -64,6 +57,13 @@ Unreleased Changes
       (https://github.com/natcap/invest/issues/1609).
     * Improved error handling when a datastack cannot be saved with relative
       paths across drives (https://github.com/natcap/invest/issues/1608).
+* Coastal Vulnerability
+    * Fixed a regression where an AOI with multiple features could raise a
+      TypeError after intersecting with the landmass polygon.
+      https://github.com/natcap/invest/issues/1657
+* Forest Carbon Edge Effects
+    * Updating vector reprojection to allow partial reprojection.  Related to
+      https://github.com/natcap/invest/issues/1645
 * Habitat Quality
     * Access raster is now generated from the reprojected access vector
       (https://github.com/natcap/invest/issues/1615).
@@ -72,6 +72,10 @@ Unreleased Changes
 * Urban Flood Risk
     * Fields present on the input AOI vector are now retained in the output.
       (https://github.com/natcap/invest/issues/1600)
+* Urban Nature Access
+    * The model now works as expected when the user provides an LULC raster
+      that does not have a nodata value defined.
+      https://github.com/natcap/invest/issues/1293
 
 3.14.2 (2024-05-29)
 -------------------
diff --git a/src/natcap/invest/coastal_vulnerability.py b/src/natcap/invest/coastal_vulnerability.py
index ff30d3a24..a1637451a 100644
--- a/src/natcap/invest/coastal_vulnerability.py
+++ b/src/natcap/invest/coastal_vulnerability.py
@@ -21,6 +21,7 @@ import taskgraph
 from osgeo import gdal
 from osgeo import ogr
 from osgeo import osr
+from shapely.geometry import LineString, MultiLineString
 from shapely.geometry.base import BaseMultipartGeometry
 from shapely.strtree import STRtree
 
@@ -1104,7 +1105,6 @@ def prepare_landmass_line_index_and_interpolate_shore_points(
     # Get shapely geometries from landmass
     landmass_polygon_shapely_list = _ogr_to_geometry_list(landmass_vector_path)
     landmass_shapely = shapely.ops.unary_union(landmass_polygon_shapely_list)
-
     landmass_polygon_shapely_list = None
 
     # store polygon geom for point-in-poly check later in ray-casting
@@ -1170,19 +1170,13 @@ def prepare_landmass_line_index_and_interpolate_shore_points(
             if aoi_shapely_prepped.intersects(landmass_line):
                 intersected_shapely_geom = aoi_shapely.intersection(
                     landmass_line)
-                if intersected_shapely_geom.geom_type == 'LineString':
-                    lines_in_aoi_list.append(intersected_shapely_geom)
-                elif intersected_shapely_geom.geom_type == 'MultiLineString':
-                    shapely_geom_explode = [
-                        shapely.geometry.LineString(x)
-                        for x in intersected_shapely_geom]
-
-                    lines_in_aoi_list.extend(shapely_geom_explode)
-                else:
-                    # intersection could generate a point geom
-                    # or if somehow the intersection is empty,
-                    # type will be GeometryCollection.
-                    continue
+                # intersection could generate a point geom,
+                # or if somehow the intersection is empty,
+                # type will be GeometryCollection.
+                if isinstance(intersected_shapely_geom,
+                              (LineString, MultiLineString)):
+                    lines_in_aoi_list.extend(
+                        _list_geometry(intersected_shapely_geom))
 
         # if none of the lines were disjoint before this linemerge,
         # unioned_line will now be a LineString.


Test patch:
diff --git a/tests/test_coastal_vulnerability.py b/tests/test_coastal_vulnerability.py
index 0b3be853a..a6e6f6088 100644
--- a/tests/test_coastal_vulnerability.py
+++ b/tests/test_coastal_vulnerability.py
@@ -1190,11 +1190,15 @@ class CoastalVulnerabilityTests(unittest.TestCase):
         srs = osr.SpatialReference()
         srs.ImportFromEPSG(26910)  # UTM Zone 10N
         wkt = srs.ExportToWkt()
+        # These two disjoint AOI polygons intersect the same landmass line
+        # segment. This tests an edge case where a MultiLineString
+        # geometry is created when landmass lines are clipped by the AOI.
         poly_a = Polygon([
             (-200, -200), (-100, -200), (-100, -100), (-200, -100),
             (-200, -200)])
         poly_b = Polygon([
-            (100, 100), (200, 100), (200, 200), (100, 200), (100, 100)])
+            (100, -200), (200, -200), (200, -100), (100, -100),
+            (100, -200)])
 
         pygeoprocessing.shapely_geometry_to_vector(
             [poly_a, poly_b], aoi_path, wkt, 'GeoJSON')
@@ -1358,6 +1362,48 @@ class CoastalVulnerabilityTests(unittest.TestCase):
                 aoi_path, landmass_path, model_resolution, target_vector_path,
                 polygon_pickle, lines_pickle, lines_rtree)
 
+    def test_prepare_landmass_invalid_geometry(self):
+        """CV: test handling invalid geometries in landmass vector."""
+        from natcap.invest import coastal_vulnerability
+        aoi_path = os.path.join(self.workspace_dir, 'aoi.geojson')
+        srs = osr.SpatialReference()
+        srs.ImportFromEPSG(26910)  # UTM Zone 10N
+        wkt = srs.ExportToWkt()
+
+        aoi_geometries = [Polygon([
+            (-200, -200), (200, -200), (200, 200), (-200, 200), (-200, -200)])]
+        pygeoprocessing.shapely_geometry_to_vector(
+            aoi_geometries, aoi_path, wkt, 'GeoJSON')
+
+        landmass_vector_path = os.path.join(self.workspace_dir, 'vector.gpkg')
+        n_features = make_vector_of_invalid_geoms(landmass_vector_path)
+
+        target_polygon_pickle_path = os.path.join(
+            self.workspace_dir, 'polygon.pickle')
+        target_lines_pickle_path = os.path.join(
+            self.workspace_dir, 'lines.pickle')
+        target_rtree_path = os.path.join(self.workspace_dir, 'rtree.dat')
+        # Create rtree files to exercise the function's logic of removing
+        # pre-exisiting files
+        target_rtree_path_base = os.path.splitext(target_rtree_path)[0]
+        open(target_rtree_path, 'a').close()
+        open(target_rtree_path_base + '.idx', 'a').close()
+
+        model_resolution = 100
+        target_vector_path = os.path.join(
+            self.workspace_dir, 'temp-shore-pts.gpkg')
+        coastal_vulnerability.prepare_landmass_line_index_and_interpolate_shore_points(
+            aoi_path, landmass_vector_path, model_resolution,
+            target_vector_path, target_polygon_pickle_path,
+            target_lines_pickle_path, target_rtree_path)
+
+        with open(target_polygon_pickle_path, 'rb') as polygon_file:
+            shapely_geom = pickle.load(polygon_file)
+
+        # Expect 1 input geometry to be skipped, and the rest to be in
+        # shapely_geom_list.
+        self.assertTrue(len(shapely_geom.geoms) == n_features - 1)
+
     def test_no_wwiii_coverage(self):
         """CV: test exception when shore points are outside max wwiii dist."""
         from natcap.invest import coastal_vulnerability
@@ -1434,48 +1480,6 @@ class CoastalVulnerabilityTests(unittest.TestCase):
         layer = None
         vector = None
 
-    def test_prepare_landmass_invalid_geometry(self):
-        """CV: test handling invalid geometries in landmass vector."""
-        from natcap.invest import coastal_vulnerability
-        aoi_path = os.path.join(self.workspace_dir, 'aoi.geojson')
-        srs = osr.SpatialReference()
-        srs.ImportFromEPSG(26910)  # UTM Zone 10N
-        wkt = srs.ExportToWkt()
-
-        aoi_geometries = [Polygon([
-            (-200, -200), (200, -200), (200, 200), (-200, 200), (-200, -200)])]
-        pygeoprocessing.shapely_geometry_to_vector(
-            aoi_geometries, aoi_path, wkt, 'GeoJSON')
-
-        landmass_vector_path = os.path.join(self.workspace_dir, 'vector.gpkg')
-        n_features = make_vector_of_invalid_geoms(landmass_vector_path)
-
-        target_polygon_pickle_path = os.path.join(
-            self.workspace_dir, 'polygon.pickle')
-        target_lines_pickle_path = os.path.join(
-            self.workspace_dir, 'lines.pickle')
-        target_rtree_path = os.path.join(self.workspace_dir, 'rtree.dat')
-        # Create rtree files to exercise the function's logic of removing
-        # pre-exisiting files
-        target_rtree_path_base = os.path.splitext(target_rtree_path)[0]
-        open(target_rtree_path, 'a').close()
-        open(target_rtree_path_base + '.idx', 'a').close()
-
-        model_resolution = 100
-        target_vector_path = os.path.join(
-            self.workspace_dir, 'temp-shore-pts.gpkg')
-        coastal_vulnerability.prepare_landmass_line_index_and_interpolate_shore_points(
-            aoi_path, landmass_vector_path, model_resolution,
-            target_vector_path, target_polygon_pickle_path,
-            target_lines_pickle_path, target_rtree_path)
-
-        with open(target_polygon_pickle_path, 'rb') as polygon_file:
-            shapely_geom = pickle.load(polygon_file)
-
-        # Expect 1 input geometry to be skipped, and the rest to be in
-        # shapely_geom_list.
-        self.assertTrue(len(shapely_geom.geoms) == n_features - 1)
-
     def test_clip_project_already_projected_raster(self):
         """CV: test clip_and_project_raster on an already projected raster."""
         from natcap.invest import coastal_vulnerability