Merge pull request natcap#1658 from davemfish/bugfix/CV-1657-multilin…

…estring

CV bugfix: handling multilinestrings

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

@@ -1190,11 +1190,15 @@ def test_aoi_multiple_features(self):
         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 @@ def test_aoi_invalid_geometry(self):
                 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 @@ def test_projected_wwiii_input(self):
         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