Repository: natcap/invest
Pull request #1064 (2022-08-25 22:56:49)

Handling a ZeroDivisionError case in HRA.
This case occurs when a subregion might not overlap a particular habitat
or stressor.  RE:#819

I didn't update HISTORY because this HRA rewrite hasn't yet been released.  Tests have been updated to cover this case, and have been confirmed to trigger the issue in the absence of the fix.



Patch:
diff --git a/src/natcap/invest/hra.py b/src/natcap/invest/hra.py
index 3534782c5..31aed8da8 100644
--- a/src/natcap/invest/hra.py
+++ b/src/natcap/invest/hra.py
@@ -1197,8 +1197,11 @@ def _create_summary_statistics_file(
             for prefix in ('E', 'C', 'R'):
                 record[f'{prefix}_MIN'] = float(stats[f'{prefix}_MIN'])
                 record[f'{prefix}_MAX'] = float(stats[f'{prefix}_MAX'])
-                record[f'{prefix}_MEAN'] = float(
-                    stats[f'{prefix}_SUM'] / stats[f'{prefix}_N_PIXELS'])
+                mean = 0
+                if stats[f'{prefix}_N_PIXELS'] > 0:  # avoid dividing by 0
+                    mean = float(
+                        stats[f'{prefix}_SUM'] / stats[f'{prefix}_N_PIXELS'])
+                record[f'{prefix}_MEAN'] = mean
             records.append(record)
             all_subregion_stats[(habitat, stressor, subregion_id)] = stats
 
@@ -1241,8 +1244,11 @@ def _create_summary_statistics_file(
                 sums_across_stressors['R_N_LOW'] / reclassified_count) * 100
 
             for prefix in ('E', 'C', 'R'):
-                record[f'{prefix}_MEAN'] = (
-                    record[f'{prefix}_SUM'] / record[f'{prefix}_N_PIXELS'])
+                mean = 0
+                if record[f'{prefix}_N_PIXELS'] > 0:
+                    mean = (
+                        record[f'{prefix}_SUM'] / record[f'{prefix}_N_PIXELS'])
+                record[f'{prefix}_MEAN'] = mean
                 del record[f'{prefix}_SUM']
                 del record[f'{prefix}_N_PIXELS']
 


Test patch:
diff --git a/tests/test_hra.py b/tests/test_hra.py
index d6952a884..17fd7f845 100644
--- a/tests/test_hra.py
+++ b/tests/test_hra.py
@@ -1006,10 +1006,17 @@ class HRAModelTests(unittest.TestCase):
             'n_overlapping_stressors': 2,
             'visualize_outputs': False,
         }
-        aoi_geoms = [shapely.geometry.box(
-            *shapely.geometry.Point(ORIGIN).buffer(1000).bounds)]
+        aoi_geoms = [
+            shapely.geometry.box(  # This geometry covers all areas
+                *shapely.geometry.Point(ORIGIN).buffer(100).bounds),
+            shapely.geometry.box(  # Geometry covers only 1 stressor
+                *shapely.geometry.Point(
+                    (ORIGIN[0], ORIGIN[1]-200)).buffer(100).bounds)
+        ]
         pygeoprocessing.shapely_geometry_to_vector(
-            aoi_geoms, args['aoi_vector_path'], SRS_WKT, 'ESRI Shapefile')
+            aoi_geoms, args['aoi_vector_path'], SRS_WKT, 'ESRI Shapefile',
+            fields={'name': ogr.OFTString}, attribute_list=[
+                {'name': 'wholearea'}, {'name': 'noarea'}])
 
         with open(args['criteria_table_path'], 'w') as criteria_table:
             criteria_table.write(textwrap.dedent(
@@ -1169,7 +1176,6 @@ class HRAModelTests(unittest.TestCase):
                 pygeoprocessing.geoprocessing.raster_to_numpy_array(
                     rasterized_path))
 
-
     def test_model_habitat_mismatch(self):
         """HRA: check errors when habitats are mismatched."""
         from natcap.invest import hra