Repository: natcap/invest (python)

Issue:
Forest Carbon test data: invalid geometries after reprojection
Enabling `gdal.UseExceptions` made this issue apparent. Two tests (`test_carbon_dup_output` and `test_carbon_full`) started failing because of invalid geometries in the vector resulting from reprojecting the original model data. For #1570 I'm going to just skip over any invalid geometries, since that's effectively what was happening before, but we should check why these are being produced and if there's any way to avoid it.

Pull request:
Forest Carbon: Global model vector clipping and reprojection
## Description
Fixes #1689 

### Background:
Enabling `gdal.UseExceptions` led to the discovery that invalid geometries were being created during reprojection of the Forest Carbon "Global Regression Models" vector (`forest_carbon_edge_regression_model_parameters.shp`). This was happening even with partial reprojection enabled. Digging into it, it seems partial reprojection doesn't necessarily remove entire geometries, just points that can't be reprojected; as such, you can end up with geometries that have insufficient points to be valid polygons.

Based on recommendations from @emlys and @phargogh, this PR does away with partial reprojection and instead clips the global vector prior to reprojection, so that only those polygons relevant to the model run need be reprojected. Clipping is based on the LULC raster's bounding box, plus a buffer (explained below). 

While this model does include an AIO input, it's optional and is only used to calculate aggregate statistics at the end of the model run. For consistency, I'm ignoring it for the purposes of clipping. 

### Bounding Box Buffering: 
The Global Regression Model vector data is used to construct a kd-tree. When querying the kd-tree, we set a `distance_upper_bound`. To ensure the clipped vector will contain all of the relevant data to the querying process, I'm buffering the LULC bounding box by the same `DISTANCE_UPPER_BOUND` value (currently set to 500km). 

I'm also keeping all polygons that intersect with the clipping polygon, rather than using `gdal.Layer.Clip()`, in order to preserve the shape of the polygons (this is important because we use the centroid when constructing the kd-tree). 

To summarize my clipping methodology: 
- Get the bounding box of the LULC raster; buffer it by `DISTANCE_UPPER_BOUND`
- Reproject the buffered bounding box to the coordinate system of the vector
- Copy all polygons from the global vector that intersect with the bounding box to a new vector
- Pass the clipped vector along to `_build_spatial_index`, to be reprojected and used to create the kd-tree

### Additional notes:
The User's Guide says, "Note that all spatial inputs must be in the same projected coordinate system and in linear meter units." However, we weren't validating the units within the model. I've updated the model spec to include `"projection_units": u.meter` for validation purposes.

Digging into this model, I also noticed a contradiction. The user's guide says, related to the "Number Of Points To Average," (in "How it works"): 
"The user can designate the number of local models used in the interpolation scheme which is defaulted to 10 but can range anywhere from 1 (only closest point) to 2635 (every regression model on the planet)."

However, within the code itself, we specify the aforementioned `DISTANCE_UPPER_BOUND` of 500km that we use when querying the kd-tree (there's a note that reads, "grid cells are 100km. Becky says 500km is a good upper bound to search"). We may want to update the User's Guide to reflect this, especially now that we're using that `DISTANCE_UPPER_BOUND` as a buffer when clipping the global vector to the LULC raster's bounding box.

## Checklist
- [ ] 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/src/natcap/invest/forest_carbon_edge_effect.py b/src/natcap/invest/forest_carbon_edge_effect.py
index 1b70a8eb8..711d104a8 100644
--- a/src/natcap/invest/forest_carbon_edge_effect.py
+++ b/src/natcap/invest/forest_carbon_edge_effect.py
@@ -12,6 +12,10 @@ import numpy
 import pandas
 import pygeoprocessing
 import scipy.spatial
+import shapely.errors
+import shapely.geometry
+import shapely.prepared
+import shapely.wkb
 import taskgraph
 from osgeo import gdal
 from osgeo import ogr
@@ -115,7 +119,8 @@ MODEL_SPEC = {
             "about": spec_utils.LULC['about'] + " " + gettext(
                 "All values in this raster must "
                 "have corresponding entries in the Biophysical Table."),
-            "projected": True
+            "projected": True,
+            "projection_units": u.meter
         },
         "pools_to_calculate": {
             "type": "option_string",
@@ -167,6 +172,8 @@ MODEL_SPEC = {
                         "Used only for the asymptotic model.")}
             },
             "geometries": spec_utils.POLYGONS,
+            "projected": True,
+            "projection_units": u.meter,
             "required": "compute_forest_edge_effects",
             "about": gettext(
                 "Map storing the optimal regression model for each tropical "
@@ -251,6 +258,13 @@ MODEL_SPEC = {
                     "bands": {1: {
                         "type": "number", "units": u.metric_ton/u.hectare
                     }}
+                },
+                "regression_model_params_clipped.shp": {
+                    "about": (
+                        "The Global Regression Models shapefile clipped "
+                        "to the study area."),
+                    "geometries": spec_utils.POLYGONS,
+                    "fields": {}
                 }
             }
         },
@@ -409,6 +423,9 @@ def execute(args):
             f'tropical_forest_edge_carbon_stocks{file_suffix}.tif')
         output_file_registry['non_forest_mask'] = os.path.join(
             intermediate_dir, f'non_forest_mask{file_suffix}.tif')
+        output_file_registry['tropical_forest_edge_clipped'] = os.path.join(
+            intermediate_dir,
+            f'regression_model_params_clipped{file_suffix}.shp')
 
     # Map non-forest landcover codes to carbon biomasses
     LOGGER.info('Calculating direct mapped carbon stocks')
@@ -445,15 +462,26 @@ def execute(args):
                               output_file_registry['non_forest_mask']],
             task_name='map_distance_from_forest_edge')
 
+        # Clip global regression model vector to LULC raster bounding box
+        LOGGER.info('Clipping global forest carbon edge regression models vector')
+        clip_forest_edge_carbon_vector_task = task_graph.add_task(
+            func=_clip_global_regression_models_vector,
+            args=(args['lulc_raster_path'],
+                  args['tropical_forest_edge_carbon_model_vector_path'],
+                  output_file_registry['tropical_forest_edge_clipped']),
+            target_path_list=[output_file_registry['tropical_forest_edge_clipped']],
+            task_name='clip_forest_edge_carbon_vector')
+
         # Build spatial index for gridded global model for closest 3 points
         LOGGER.info('Building spatial index for forest edge models.')
         build_spatial_index_task = task_graph.add_task(
             func=_build_spatial_index,
             args=(args['lulc_raster_path'], intermediate_dir,
-                  args['tropical_forest_edge_carbon_model_vector_path'],
+                  output_file_registry['tropical_forest_edge_clipped'],
                   output_file_registry['spatial_index_pickle']),
             target_path_list=[output_file_registry['spatial_index_pickle']],
-            task_name='build_spatial_index')
+            task_name='build_spatial_index',
+            dependent_task_list=[clip_forest_edge_carbon_vector_task])
 
         # calculate the carbon edge effect on forests
         LOGGER.info('Calculating forest edge carbon')
@@ -721,7 +749,8 @@ def _map_distance_from_tropical_forest_edge(
     edge_distance_raster = gdal.OpenEx(edge_distance_path, gdal.GA_Update)
     edge_distance_band = edge_distance_raster.GetRasterBand(1)
 
-    for offset_dict in pygeoprocessing.iterblocks((base_lulc_raster_path, 1), offset_only=True):
+    for offset_dict in pygeoprocessing.iterblocks((base_lulc_raster_path, 1),
+            offset_only=True):
         # where LULC has nodata, overwrite edge distance with nodata value
         lulc_block = lulc_band.ReadAsArray(**offset_dict)
         distance_block = edge_distance_band.ReadAsArray(**offset_dict)
@@ -733,6 +762,93 @@ def _map_distance_from_tropical_forest_edge(
             yoff=offset_dict['yoff'])
 
 
+def _clip_global_regression_models_vector(
+        lulc_raster_path, source_vector_path, target_vector_path):
+    """Clip the global carbon edge model shapefile
+
+    Clip the shapefile containing the global carbon edge model parameters
+    to the bounding box of the LULC raster (representing the target study area)
+    plus a buffer to account for the kd-tree lookup DISTANCE_UPPER_BOUND
+
+    Args:
+        lulc_raster_path (string): path to a raster that is used to define the
+            bounding box to use for clipping.
+        source_vector_path (string): a path to an OGR shapefile to be clipped.
+        target_vector_path (string): a path to an OGR shapefile to store the
+            clipped vector.
+
+    Returns:
+        None
+
+    """
+    raster_info = pygeoprocessing.get_raster_info(lulc_raster_path)
+    vector_info = pygeoprocessing.get_vector_info(source_vector_path)
+
+    bbox_minx, bbox_miny, bbox_maxx, bbox_maxy = raster_info['bounding_box']
+    buffered_bb = [
+        bbox_minx - DISTANCE_UPPER_BOUND,
+        bbox_miny - DISTANCE_UPPER_BOUND,
+        bbox_maxx + DISTANCE_UPPER_BOUND,
+        bbox_maxy + DISTANCE_UPPER_BOUND,
+    ]
+
+    # Reproject the LULC bounding box to the vector's projection for clipping
+    mask_bb = pygeoprocessing.transform_bounding_box(buffered_bb,
+        raster_info['projection_wkt'], vector_info['projection_wkt'])
+    shapely_mask = shapely.prepared.prep(shapely.geometry.box(*mask_bb))
+
+    base_vector = gdal.OpenEx(source_vector_path, gdal.OF_VECTOR)
+    base_layer = base_vector.GetLayer()
+    base_layer_defn = base_layer.GetLayerDefn()
+    base_geom_type = base_layer.GetGeomType()
+
+    target_driver = gdal.GetDriverByName('ESRI Shapefile')
+    target_vector = target_driver.Create(
+        target_vector_path, 0, 0, 0, gdal.GDT_Unknown)
+    target_layer = target_vector.CreateLayer(
+        base_layer_defn.GetName(), base_layer.GetSpatialRef(), base_geom_type)
+    target_layer.CreateFields(base_layer.schema)
+
+    target_layer.StartTransaction()
+    invalid_feature_count = 0
+    for feature in base_layer:
+        invalid = False
+        geometry = feature.GetGeometryRef()
+        try:
+            shapely_geom = shapely.wkb.loads(bytes(geometry.ExportToWkb()))
+        # Catch invalid geometries that cannot be loaded by Shapely;
+        # e.g. polygons with too few points for their type
+        except shapely.errors.ShapelyError:
+            invalid = True
+        else:
+            if shapely_geom.is_valid:
+                # Check for intersection rather than use gdal.Layer.Clip()
+                # to preserve the shape of the polygons (we use the centroid
+                # when constructing the kd-tree)
+                if shapely_mask.intersects(shapely_geom):
+                    new_feature = ogr.Feature(target_layer.GetLayerDefn())
+                    new_feature.SetGeometry(ogr.CreateGeometryFromWkb(
+                        shapely_geom.wkb))
+                    for field_name, field_value in feature.items().items():
+                        new_feature.SetField(field_name, field_value)
+                    target_layer.CreateFeature(new_feature)
+            else:
+                invalid = True
+        finally:
+            if invalid:
+                invalid_feature_count += 1
+                LOGGER.warning(
+                    f"The geometry at feature {feature.GetFID()} is invalid "
+                    "and will be skipped.")
+
+    target_layer.CommitTransaction()
+
+    if invalid_feature_count:
+        LOGGER.warning(
+            f"{invalid_feature_count} features in {source_vector_path} "
+            "were found to be invalid during clipping and were skipped.")
+
+
 def _build_spatial_index(
         base_raster_path, local_model_dir,
         tropical_forest_edge_carbon_model_vector_path,
@@ -770,10 +886,9 @@ def _build_spatial_index(
     lulc_projection_wkt = pygeoprocessing.get_raster_info(
         base_raster_path)['projection_wkt']
 
-    with utils._set_gdal_configuration('OGR_ENABLE_PARTIAL_REPROJECTION', 'TRUE'):
-        pygeoprocessing.reproject_vector(
-            tropical_forest_edge_carbon_model_vector_path, lulc_projection_wkt,
-            carbon_model_reproject_path)
+    pygeoprocessing.reproject_vector(
+        tropical_forest_edge_carbon_model_vector_path, lulc_projection_wkt,
+        carbon_model_reproject_path)
 
     model_vector = gdal.OpenEx(carbon_model_reproject_path)
     model_layer = model_vector.GetLayer()
@@ -785,17 +900,14 @@ def _build_spatial_index(
     # put all the polygons in the kd_tree because it's fast and simple
     for poly_feature in model_layer:
         poly_geom = poly_feature.GetGeometryRef()
-        if poly_geom.IsValid():
-            poly_centroid = poly_geom.Centroid()
-            # put in row/col order since rasters are row/col indexed
-            kd_points.append([poly_centroid.GetY(), poly_centroid.GetX()])
-
-            theta_model_parameters.append([
-                poly_feature.GetField(feature_id) for feature_id in
-                ['theta1', 'theta2', 'theta3']])
-            method_model_parameter.append(poly_feature.GetField('method'))
-        else:
-            LOGGER.warning(f'skipping invalid geometry {poly_geom}')
+        poly_centroid = poly_geom.Centroid()
+        # put in row/col order since rasters are row/col indexed
+        kd_points.append([poly_centroid.GetY(), poly_centroid.GetX()])
+
+        theta_model_parameters.append([
+            poly_feature.GetField(feature_id) for feature_id in
+            ['theta1', 'theta2', 'theta3']])
+        method_model_parameter.append(poly_feature.GetField('method'))
 
     method_model_parameter = numpy.array(
         method_model_parameter, dtype=numpy.int32)


Test patch:
diff --git a/tests/test_forest_carbon_edge.py b/tests/test_forest_carbon_edge.py
index 09435dace..93fd428ba 100644
--- a/tests/test_forest_carbon_edge.py
+++ b/tests/test_forest_carbon_edge.py
@@ -1,17 +1,18 @@
 """Module for Regression Testing the InVEST Forest Carbon Edge model."""
-import unittest
-import tempfile
-import shutil
 import os
+import pickle
+import shutil
+import tempfile
+import unittest
 
 import numpy
 import pandas
-from osgeo import gdal, osr, ogr
-from shapely.geometry import Polygon
-import pygeoprocessing
-import pickle
-
 import pygeoprocessing
+import shapely
+from osgeo import gdal
+from osgeo import ogr
+from osgeo import osr
+from shapely.geometry import Polygon
 
 gdal.UseExceptions()
 REGRESSION_DATA = os.path.join(
@@ -530,6 +531,146 @@ class ForestCarbonEdgeTests(unittest.TestCase):
         self.assertAlmostEqual(actual_output[3, 5], 274.47815)
         self.assertAlmostEqual(actual_output[0, 1], -1)
 
+    def test_invalid_geometries(self):
+        """Forest Carbon: Check handling of invalid vector geometries."""
+        from natcap.invest import forest_carbon_edge_effect
+        srs = osr.SpatialReference()
+        srs.ImportFromEPSG(26910)
+        projection_wkt = srs.ExportToWkt()
+
+        lulc_matrix = numpy.ones((5,5), dtype=numpy.int8)
+        lulc_raster_path = os.path.join(self.workspace_dir, 'lulc.tif')
+        pygeoprocessing.numpy_array_to_raster(
+            lulc_matrix, 255, (1, -1), (461262, 4923269), projection_wkt,
+            lulc_raster_path)
+
+        fields_polys = {'geom_id': ogr.OFTInteger}
+        attrs_polys = [
+            {'geom_id': 1},
+            {'geom_id': 2},
+            {'geom_id': 3}
+        ]
+
+        # valid polygon
+        valid_poly = Polygon([
+            (461265, 4923269), (461267, 4923269), (461267, 4923267),
+            (461265, 4923267), (461265, 4923269)])
+        self.assertTrue(shapely.is_valid(valid_poly))
+
+        # invalid bowtie polygon
+        invalid_bowtie_poly = Polygon([
+            (461262, 4923267), (461264, 4923267), (461262, 4923265),
+            (461264, 4923265), (461262, 4923267)])
+        self.assertFalse(shapely.is_valid(invalid_bowtie_poly))
+
+        # invalid bowtie polygon with vertex in the middle
+        invalid_alt_bowtie_poly = Polygon([
+            (461262, 4923267), (461264, 4923267), (461263, 4923266),
+            (461264, 4923265), (461262, 4923265), (461263, 4923266),
+            (461262, 4923267)])
+        self.assertFalse(shapely.is_valid(invalid_alt_bowtie_poly))
+
+        test_vector_path = os.path.join(self.workspace_dir, 'vector.gpkg')
+        pygeoprocessing.shapely_geometry_to_vector(
+            [valid_poly, invalid_bowtie_poly,
+             invalid_alt_bowtie_poly],
+            test_vector_path, projection_wkt, 'GPKG',
+            fields=fields_polys,
+            attribute_list=attrs_polys)
+
+        test_vector = gdal.OpenEx(
+            test_vector_path, gdal.OF_VECTOR | gdal.GA_ReadOnly)
+        test_layer = test_vector.GetLayer()
+        self.assertEqual(test_layer.GetFeatureCount(), 3)
+        test_layer = None
+        test_vector = None
+
+        target_vector_path = os.path.join(
+            self.workspace_dir, 'checked_geometries.gpkg')
+        forest_carbon_edge_effect._clip_global_regression_models_vector(
+            lulc_raster_path, test_vector_path, target_vector_path)
+
+        target_vector = gdal.OpenEx(target_vector_path, gdal.OF_VECTOR)
+        target_layer = target_vector.GetLayer()
+        self.assertEqual(
+            target_layer.GetFeatureCount(), 1)
+        # valid_polygon geom_id == 1
+        self.assertEqual(target_layer.GetFeature(0).GetField('geom_id'), 1)
+
+        target_layer = None
+        target_vector = None
+
+    def test_vector_clipping_buffer(self):
+        """Forest Carbon: Check DISTANCE_UPPER_BOUND buffer."""
+        from natcap.invest import forest_carbon_edge_effect
+        srs = osr.SpatialReference()
+        srs.ImportFromEPSG(26910)
+        projection_wkt = srs.ExportToWkt()
+
+        lulc_matrix = numpy.ones((5,5), dtype=numpy.int8)
+        lulc_raster_path = os.path.join(self.workspace_dir, 'lulc.tif')
+        pygeoprocessing.numpy_array_to_raster(
+            lulc_matrix, 255, (1000, -1000), (461261, 4923270), projection_wkt,
+            lulc_raster_path)
+
+        # get LULC bouding box; will be used to construct polygons
+        raster_info = pygeoprocessing.get_raster_info(lulc_raster_path)
+        bbox_minx, bbox_miny, bbox_maxx, bbox_maxy = raster_info['bounding_box']
+
+        fields_polys = {'geom_id': ogr.OFTInteger}
+        attrs_polys = [
+            {'geom_id': 1},
+            {'geom_id': 2},
+            {'geom_id': 3},
+            {'geom_id': 4}
+        ]
+
+        # polygon that intersects with the lulc
+        in_bbox_poly = shapely.geometry.box(
+            bbox_maxx - 150, bbox_maxy - 150,
+            bbox_maxx - 50, bbox_maxy - 50)
+        self.assertTrue(shapely.is_valid(in_bbox_poly))
+
+        # polygon outside of lulc bb, but within buffer distance
+        in_buffer_poly = shapely.geometry.box(
+            bbox_maxx + 50, bbox_maxy + 50,
+            bbox_maxx + 150, bbox_maxy + 150)
+        self.assertTrue(shapely.is_valid(in_buffer_poly))
+
+        buffer = forest_carbon_edge_effect.DISTANCE_UPPER_BOUND
+        # polygon on the edge of the buffer
+        buffer_edge_poly = shapely.geometry.box(
+            bbox_maxx + buffer + 50, bbox_maxy + buffer + 50,
+            bbox_maxx + buffer - 50, bbox_maxy + buffer - 50)
+        self.assertTrue(shapely.is_valid(buffer_edge_poly))
+
+        # polygon outside of buffer distance
+        out_of_buffer_poly = shapely.geometry.box(
+            bbox_maxx + buffer + 50, bbox_maxy + buffer + 50,
+            bbox_maxx + buffer + 150, bbox_maxy + buffer + 150)
+        self.assertTrue(shapely.is_valid(out_of_buffer_poly))
+
+        test_vector_path = os.path.join(self.workspace_dir, 'vector.gpkg')
+        pygeoprocessing.shapely_geometry_to_vector(
+            [in_bbox_poly, in_buffer_poly,
+             buffer_edge_poly, out_of_buffer_poly],
+            test_vector_path, projection_wkt, 'GPKG',
+            fields=fields_polys,
+            attribute_list=attrs_polys)
+
+        target_vector_path = os.path.join(
+            self.workspace_dir, 'checked_geometries.gpkg')
+        forest_carbon_edge_effect._clip_global_regression_models_vector(
+            lulc_raster_path, test_vector_path, target_vector_path)
+
+        vector = gdal.OpenEx(
+            target_vector_path, gdal.OF_VECTOR | gdal.GA_ReadOnly)
+        layer = vector.GetLayer()
+        self.assertEqual(layer.GetFeatureCount(), 3)
+
+        feat_ids = [feat.GetField('geom_id') for feat in layer]
+        self.assertEqual(feat_ids, [1, 2, 3])
+
     @staticmethod
     def _test_same_files(base_list_path, directory_path):
         """Assert files in `base_list_path` are in `directory_path`.