Repository: natcap/invest (python)

Issue:
Shapely 1.8 deprecation warnings
It looks like shapely `1.8.0` is out and it's a transitional release in prep for `2.0`. 
https://shapely.readthedocs.io/en/stable/migration.html#migrating-to-shapely-1-8-2-0

There are some new deprecation warnings we can deal with while we're on `1.8.0`. And we should probably restrict our requirements to `<2.0`. I'm not sure exactly when 2.0 will drop.

I observed some warnings when running the CV sample data.


Pull request:
Resolve shapely 1.8.0 deprecation warnings
# Description
Fixes #763 
* `cascaded_union` -> `unary_union`
* multipart geometries are no longer iterable
* can no longer assign attributes to geometries

# Checklist
- [ ] ~Updated HISTORY.rst (if these changes are user-facing)~
- [ ] ~Updated the user's guide (if needed)~


Patch:
diff --git a/src/natcap/invest/coastal_vulnerability.py b/src/natcap/invest/coastal_vulnerability.py
index d7e1c2ac7..062b3c24a 100644
--- a/src/natcap/invest/coastal_vulnerability.py
+++ b/src/natcap/invest/coastal_vulnerability.py
@@ -17,6 +17,7 @@ import shapely.ops
 import shapely.speedups
 import shapely.errors
 from shapely.strtree import STRtree
+from shapely.geometry.base import BaseMultipartGeometry
 import pygeoprocessing
 import taskgraph
 
@@ -675,6 +676,21 @@ def execute(args):
         target_exposure_csv_path)
 
 
+def _list_geometry(geom):
+    """Iterate over a possibly multipart geometry.
+
+
+    Args:
+        geom (shapely.Geometry): geometry to iterate over
+
+    Returns:
+        list of shapely single-part geometries
+    """
+    if isinstance(geom, BaseMultipartGeometry):
+        return list(geom.geoms)
+    return [geom]
+
+
 def prepare_landmass_line_index_and_interpolate_shore_points(
         aoi_vector_path, landmass_vector_path, model_resolution,
         target_points_path, target_polygon_pickle_path,
@@ -706,8 +722,7 @@ def prepare_landmass_line_index_and_interpolate_shore_points(
     LOGGER.info("preparing landmass geometry and interpolating shore points")
     # Get shapely geometries from landmass
     landmass_polygon_shapely_list = _ogr_to_geometry_list(landmass_vector_path)
-    landmass_shapely = shapely.ops.cascaded_union(
-        landmass_polygon_shapely_list)
+    landmass_shapely = shapely.ops.unary_union(landmass_polygon_shapely_list)
 
     landmass_polygon_shapely_list = None
 
@@ -733,7 +748,7 @@ def prepare_landmass_line_index_and_interpolate_shore_points(
     aoi_spatial_reference.ImportFromWkt(aoi_vector_info['projection_wkt'])
 
     aoi_shapely_list = _ogr_to_geometry_list(aoi_vector_path)
-    aoi_shapely = shapely.ops.cascaded_union(aoi_shapely_list)
+    aoi_shapely = shapely.ops.unary_union(aoi_shapely_list)
     aoi_shapely_prepped = shapely.prepared.prep(aoi_shapely)
 
     gpkg_driver = ogr.GetDriverByName("GPKG")
@@ -753,14 +768,8 @@ def prepare_landmass_line_index_and_interpolate_shore_points(
     # unique ID to iterate for shore points target SHORE_ID_FIELD field
     point_idx_tracker = 0
 
-    # if landmass is a single polygon, convert to list for iteration
-    try:
-        iter(landmass_shapely)
-    except(TypeError):
-        landmass_shapely = [landmass_shapely]
-
     LOGGER.info("indexing geometry of landmass and creating shore points")
-    for landmass_polygon in landmass_shapely:
+    for landmass_polygon in _list_geometry(landmass_shapely):
         lines_in_aoi_list = []
         for landmass_line in geometry_to_lines(
                 landmass_polygon, include_interiors=False):
@@ -790,19 +799,13 @@ def prepare_landmass_line_index_and_interpolate_shore_points(
                     # type will be GeometryCollection.
                     continue
 
-        unioned_line = shapely.ops.linemerge(lines_in_aoi_list)
-
         # if none of the lines were disjoint before this linemerge,
-        # unioned_line will now be a non-iterable LineString.
+        # unioned_line will now be a LineString.
         # If some lines were disjoint before the linemerge, unioned_line
-        # will now be an iterable MultiLineString.
-        try:
-            iter(unioned_line)
-        except(TypeError):
-            unioned_line = [unioned_line]
-
+        # will now be a MultiLineString.
+        unioned_line = shapely.ops.linemerge(lines_in_aoi_list)
         point_list = []
-        for line in unioned_line:
+        for line in _list_geometry(unioned_line):
             n_pts = int(math.ceil(line.length / model_resolution))
             points_along = [
                 line.interpolate(float(i) / n_pts, normalized=True)
@@ -1662,14 +1665,8 @@ def calculate_surge_exposure(
                     "Could not transform feature from %s to spatial reference"
                     "system %s", shelf_contour_path, target_srs_wkt)
                 continue
-            shapely_geom = shapely.wkb.loads(
-                bytes(clipped_geometry.ExportToWkb()))
-            if shapely_geom.type == 'MultiLineString':
-                shapely_geom_explode = [
-                    shapely.geometry.LineString(x) for x in shapely_geom]
-                shapely_geometry_list.extend(shapely_geom_explode)
-            else:
-                shapely_geometry_list.append(shapely_geom)
+            shapely_geom = shapely.wkb.loads(bytes(clipped_geometry.ExportToWkb()))
+            shapely_geometry_list.extend(_list_geometry(shapely_geom))
 
     if not shapely_geometry_list:
         raise ValueError(
@@ -2106,12 +2103,7 @@ def search_for_vector_habitat(
                         continue
                 shapely_geom = shapely.wkb.loads(
                     bytes(clipped_geometry.ExportToWkb()))
-                if shapely_geom.type == 'MultiPolygon':
-                    shapely_geom_explode = [
-                        shapely.geometry.Polygon(x) for x in shapely_geom]
-                    shapely_geometry_list.extend(shapely_geom_explode)
-                else:
-                    shapely_geometry_list.append(shapely_geom)
+                shapely_geometry_list.extend(_list_geometry(shapely_geom))
         else:
             LOGGER.warning(
                 "FID %d in %s has invalid geometry and will be excluded"
@@ -2267,15 +2259,16 @@ def calculate_geomorphology_exposure(
         geomorphology_vector_path, gdal.OF_VECTOR | gdal.GA_Update)
     geomorph_layer = geomorph_vector.GetLayer()
     LOGGER.info("Build spatial index of geomorphology segments")
+    geom_id_to_rank = {}
     for feature in geomorph_layer:
         geometry = feature.GetGeometryRef()
         shapely_geometry = shapely.wkb.loads(
             bytes(geometry.ExportToWkb()))
-        # Store the rank field value as an attribute of a shapely geometry so
-        # it's accessible when the geometry is returned from a STRtree query.
-        shapely_geometry.rank = feature.GetField('RANK')
         if shapely_geometry.is_valid:
             geomorph_line_shapely_list.append(shapely_geometry)
+            # use id because geometries aren't hashable. see note:
+            # https://shapely.readthedocs.io/en/1.8.0/manual.html#strtree.STRtree.strtree.query
+            geom_id_to_rank[id(shapely_geometry)] = feature.GetField('RANK')
         else:
             LOGGER.warning(
                 'geomorphology FID:%d is excluded due to invalid geometry',
@@ -2284,6 +2277,7 @@ def calculate_geomorphology_exposure(
         feature = None
     geomorph_layer = None
     geomorph_vector = None
+
     tree = STRtree(geomorph_line_shapely_list)
 
     # Vector to track which points ended up with the geomorphology_fill_value:
@@ -2317,8 +2311,8 @@ def calculate_geomorphology_exposure(
         found_geoms = [g for g in tree.query(poly_shapely)]
         # the tree.query returns geometries found in envelope of poly_shapely,
         # so follow-up with a check for intersection.
-        found_ranks = set(
-            [g.rank for g in found_geoms if poly_shapely.intersects(g)])
+        found_ranks = set(geom_id_to_rank[id(g)] for g in found_geoms
+                          if poly_shapely.intersects(g))
         if found_ranks:
             mean_geomorph = sum(found_ranks) / float(len(found_ranks))
         else:
diff --git a/src/natcap/invest/recreation/recmodel_client.py b/src/natcap/invest/recreation/recmodel_client.py
index d4ff713f5..8062ca379 100644
--- a/src/natcap/invest/recreation/recmodel_client.py
+++ b/src/natcap/invest/recreation/recmodel_client.py
@@ -602,7 +602,7 @@ def _grid_vector(vector_path, grid_type, cell_size, out_grid_vector_path):
         original_vector_shapes.append(wkt_feat)
     vector_layer.ResetReading()
     original_polygon = shapely.prepared.prep(
-        shapely.ops.cascaded_union(original_vector_shapes))
+        shapely.ops.unary_union(original_vector_shapes))
 
     out_grid_vector = driver.Create(
         out_grid_vector_path, 0, 0, 0, gdal.GDT_Unknown)
diff --git a/src/natcap/invest/utils.py b/src/natcap/invest/utils.py
index 52b4c2e86..b5fecfea6 100644
--- a/src/natcap/invest/utils.py
+++ b/src/natcap/invest/utils.py
@@ -727,6 +727,12 @@ def _assert_vectors_equal(
     Assert spatial reference, feature count, geometries, field names, and
     values are equal with no respect to order of field names or geometries.
 
+    Note: this assertion may fail incorrectly. The comparison algorithm
+    allows for float imprecision if geometry coords are in the same order. It
+    also allows coords to be in a different order if the numbers match exactly.
+    If the coords are in a different order and have float imprecision, the
+    assertion will fail.
+
     Args:
         actual_vector_path (string): path on disk to a gdal Vector dataset.
         expected_vector_path (string): path on disk to a gdal Vector dataset
@@ -821,18 +827,18 @@ def _assert_vectors_equal(
             actual_geom_wkt = actual_geom.ExportToWkt()
             expected_geom_shapely = loads(expected_geom_wkt)
             actual_geom_shapely = loads(actual_geom_wkt)
-            # Try comparing geoms exactly equal allowing for different
-            # geometry ordering
-            geoms_equal = expected_geom_shapely.equals(actual_geom_shapely)
+            # confusingly named, `equals` compares numbers exactly but
+            # allows coords to be in a different order.
+            # `equals_exact` compares numbers with a tolerance but
+            # coords must be in the same order.
+            geoms_equal = (
+                expected_geom_shapely.equals(actual_geom_shapely) or
+                expected_geom_shapely.equals_exact(actual_geom_shapely, 1e-6))
             if not geoms_equal:
-                # Try almost_equal allowing for precision differences
-                geoms_almost_eq = expected_geom_shapely.almost_equals(
-                    actual_geom_shapely)
-                if not geoms_almost_eq:
-                    raise AssertionError(
-                        "Vector geometry assertion fail. \n"
-                        f"Expected geometry: {expected_geom_wkt}. \n"
-                        f"Actual geometry: {actual_geom_wkt}. ")
+                raise AssertionError(
+                    "Vector geometry assertion fail. \n"
+                    f"Expected geometry: {expected_geom_wkt}. \n"
+                    f"Actual geometry: {actual_geom_wkt}. ")
 
             expected_feature = None
             actual_feature = None


Test patch:
diff --git a/tests/test_coastal_vulnerability.py b/tests/test_coastal_vulnerability.py
index 07a118036..12850b35a 100644
--- a/tests/test_coastal_vulnerability.py
+++ b/tests/test_coastal_vulnerability.py
@@ -1355,11 +1355,11 @@ class CoastalVulnerabilityTests(unittest.TestCase):
             target_lines_pickle_path, target_rtree_path)
 
         with open(target_polygon_pickle_path, 'rb') as polygon_file:
-            shapely_geom_list = pickle.load(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_list) == n_features - 1)
+        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."""