Handle GDAL exceptions raised in reproject_vector

HISTORY.rst

@@ -5,6 +5,8 @@ Unreleased Changes
 ------------------
 * Removing the ``numpy<2`` constraint for requirements.txt that should have
   been included in the 2.4.5 release. https://github.com/natcap/pygeoprocessing/issues/396
+* Handling GDAL-based ``RuntimeError``s raised during ``pygeoprocessing.reproject_vector``.
+  https://github.com/natcap/pygeoprocessing/issues/409
 
 2.4.5 (2024-10-08)
 ------------------

src/pygeoprocessing/geoprocessing.py

@@ -2242,21 +2242,29 @@ def reproject_vector(
             100.0 * float(feature_index+1) / (layer.GetFeatureCount()),
             os.path.basename(target_path))
 
-        geom = base_feature.GetGeometryRef()
-        if geom is None:
-            # we encountered this error occasionally when transforming clipped
-            # global polygons.  Not clear what is happening but perhaps a
-            # feature was retained that otherwise wouldn't have been included
-            # in the clip
-            error_count += 1
-            continue
+        try:
+            geom = base_feature.GetGeometryRef()
+            if geom is None:
+                # we encountered this error occasionally when transforming
+                # clipped global polygons.  Not clear what is happening but
+                # perhaps a feature was retained that otherwise wouldn't have
+                # been included in the clip
+                error_count += 1
+                continue
 
-        # Transform geometry into format desired for the new projection
-        error_code = geom.Transform(coord_trans)
-        if error_code != 0:  # error
-            # this could be caused by an out of range transformation
-            # whatever the case, don't put the transformed poly into the
-            # output set
+            # Transform geometry into format desired for the new projection
+            error_code = geom.Transform(coord_trans)
+            if error_code != 0:  # error
+                # this could be caused by an out of range transformation
+                # whatever the case, don't put the transformed poly into the
+                # output set
+                error_count += 1
+                continue
+        except RuntimeError as error:
+            # RuntimeError: GDAL's base error when geometries cannot be
+            # returned or transformed.
+            LOGGER.debug("Skipping feature %s due to %s",
+                         feature_index, str(error))
             error_count += 1
             continue
 

tests/test_geoprocessing.py

@@ -498,7 +498,6 @@ def test_reproject_vector(self):
             layer = None
             vector = None
 
-
     def test_reproject_vector_partial_fields(self):
         """PGP.geoprocessing: reproject vector with partial field copy."""
         # Create polygon shapefile to reproject
@@ -798,6 +797,49 @@ def test_reproject_vector_utm_to_utm(self):
             osr.SpatialReference(result_reference.ExportToWkt()).IsSame(
                 osr.SpatialReference(target_reference.ExportToWkt())))
 
+    def test_reproject_vector_handles_bad_data(self):
+        """PGP.geoprocessing: reproject vector with bad data."""
+        vector_path = os.path.join(self.workspace_dir, 'bad_vector.shp')
+        driver = ogr.GetDriverByName('ESRI Shapefile')
+        vector = driver.CreateDataSource(vector_path)
+        srs = osr.SpatialReference()
+        srs.ImportFromEPSG(26710)  # NAD27 / UTM zone 10N
+        layer = vector.CreateLayer(
+            'bad_layer', srs=srs, geom_type=ogr.wkbPoint)
+
+        # No/empty geometry
+        feature = ogr.Feature(layer.GetLayerDefn())
+        layer.CreateFeature(feature)
+
+        # Create a point at the centroid of the UTM zone 10N bounding box
+        feature = ogr.Feature(layer.GetLayerDefn())
+        feature.SetGeometry(
+            ogr.CreateGeometryFromWkt('POINT (2074757.82 7209331.79)'))
+
+        layer = None
+        vector = None
+
+        # Our target UTM zone is 59S (in NZ), so the points should not be
+        # usable.
+        target_srs = osr.SpatialReference()
+        target_srs.ImportFromEPSG(2134)  # NZGD2000 / UTM zone 59S
+        target_srs_wkt = target_srs.ExportToWkt()
+
+        target_path = os.path.join(self.workspace_dir, 'target_vector.shp')
+        pygeoprocessing.reproject_vector(
+            base_vector_path=vector_path,
+            target_projection_wkt=target_srs_wkt,
+            target_path=target_path,
+        )
+
+        # verify that both fields were skipped.
+        try:
+            target_vector = gdal.OpenEx(target_path)
+            target_layer = target_vector.GetLayer()
+            self.assertEqual(target_layer.GetFeatureCount(), 0)
+        finally:
+            target_vector = target_layer = None
+
     def test_calculate_disjoint_polygon_set(self):
         """PGP.geoprocessing: test calc_disjoint_poly no/intersection."""
         gpkg_driver = ogr.GetDriverByName('GPKG')