fixed a shapely 2.0 incompatibility; modified a test to cover it. natcap#1657

Author: davemfish

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

src/natcap/invest/coastal_vulnerability.py

@@ -21,6 +21,7 @@
 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.

tests/test_coastal_vulnerability.py

@@ -1184,17 +1184,22 @@ def test_zero_shorepoints_created(self):
     def test_aoi_multiple_features(self):
         """CV: test shore point creation in AOI with multiple features."""
         from natcap.invest import coastal_vulnerability
-        workspace_dir = self.workspace_dir
+        # workspace_dir = self.workspace_dir
+        workspace_dir = 'scratch/cv_test'
 
         aoi_path = os.path.join(workspace_dir, 'aoi.geojson')
         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')