Repository: natcap/invest (python)

Issue:
HRA outputs always high-risk when E/C pairs are entirely 0 for a habitat
A user on the forums ([link to private convo](https://community.naturalcapitalproject.org/t/hra-outputs/3495) with @phargogh and @davemfish ) has an interesting setup for their run where all of the criteria for all stressors affecting a single habitat are being ignored by setting those table values to 0.  This is resulting in some invalid values (numpy `WARNING`s from divide-by-zero), which I suspect is then throwing off the results.

Pull request description:
HRA fix for when the user skips all criteria in E or C
This PR handles a case where the user may request to skip all habitat/stressor overlap pairs when calculating exposure or consequence.  The current behavior is that you'll get a bunch of division-by-zero warnings, which then percolates up to the risk classification where `inf` values are being interpreted as valid, which they are not.

This PR fixes the issue by handling this case and using 0 as the calculated score, which works within the existing model structure.

Holding in draft until the user confirms the fix, and (realistically) until after the 3.13.0 release.

Fixes #1250 

## Checklist
- [x] Updated HISTORY.rst and link to any relevant issue (if these changes are user-facing)
- [ ] Updated the user's guide (if needed)
- [ ] Tested the affected models' UIs (if relevant)


Patch:
diff --git a/HISTORY.rst b/HISTORY.rst
index cdb084aa9..5a3d1c572 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -38,11 +38,16 @@
 Unreleased Changes
 ------------------
 
+* HRA
+    * Fixed a bug in HRA where the model would error when all exposure and
+      consequence criteria were skipped for a single habitat. The model now
+      correctly handles this case. https://github.com/natcap/invest/issues/1250
 * Scenic Quality
     * The Scenic Quality model will now raise an error when it encounters a
       geometry that is not a simple Point.  This is in line with the user's
       guide chapter.  https://github.com/natcap/invest/issues/1245
 
+
 3.13.0 (2023-03-17)
 -------------------
 * General
diff --git a/src/natcap/invest/hra.py b/src/natcap/invest/hra.py
index 91e3376f6..dee30fddf 100644
--- a/src/natcap/invest/hra.py
+++ b/src/natcap/invest/hra.py
@@ -2275,8 +2275,20 @@ def _calc_criteria(attributes_list, habitat_mask_raster_path,
         if decayed_edt_raster_path:
             numerator[valid_mask] *= decayed_edt_band.ReadAsArray(
                 **block_info)[valid_mask]
-        criterion_score[valid_mask] = (
-            numerator[valid_mask] / denominator[valid_mask])
+
+        # It is possible for the user to skip all pairwise calculations by
+        # setting the rating to 0 for every habitat/stressor combination.
+        # Doing so will leave the denominator at 0, resulting in a numpy
+        # divide-by-zero warning and +/- inf pixel values.
+        # JD is making the call to instead write out a raster filled with 0
+        # everywhere it is valid.
+        if numpy.sum(denominator[valid_mask]) == 0:
+            criterion_score[valid_mask] = 0
+        else:
+            # This is the normal calculation when the user has defined at least
+            # one rating for a habitat/stressor pair.
+            criterion_score[valid_mask] = (
+                numerator[valid_mask] / denominator[valid_mask])
 
         target_criterion_band.WriteArray(criterion_score,
                                          xoff=block_info['xoff'],


Test patch:
diff --git a/tests/test_hra.py b/tests/test_hra.py
index 7ce7b3e9d..d7f096d2e 100644
--- a/tests/test_hra.py
+++ b/tests/test_hra.py
@@ -77,6 +77,44 @@ class HRAUnitTests(unittest.TestCase):
         numpy.testing.assert_allclose(
             exposure_array, expected_exposure_array)
 
+    def test_calc_criteria_skip_all_criteria(self):
+        """HRA: handle user skipping all criteria."""
+        from natcap.invest import hra
+
+        habitat_mask = numpy.array([
+            [0, 1, 1]], dtype=numpy.uint8)
+        habitat_mask_path = os.path.join(self.workspace_dir,
+                                         'habitat_mask.tif')
+        pygeoprocessing.numpy_array_to_raster(
+            habitat_mask, 255, (30, -30), ORIGIN, SRS_WKT, habitat_mask_path)
+
+        decayed_distance_array = numpy.array([
+            [0, 1, 1]], dtype=numpy.float32)
+        decayed_distance_raster_path = os.path.join(self.workspace_dir,
+                                                    'decayed_dist.tif')
+        pygeoprocessing.numpy_array_to_raster(
+            decayed_distance_array, -1, (30, -30), ORIGIN, SRS_WKT,
+            decayed_distance_raster_path)
+
+        attributes_list = [
+            {'rating': 0, 'dq': 3, 'weight': 3},
+            {'rating': 0, 'dq': 2, 'weight': 1},
+            {'rating': 0, 'dq': 3, 'weight': 3},
+            {'rating': 0, 'dq': 3, 'weight': 3},
+        ]
+        target_exposure_path = os.path.join(self.workspace_dir, 'exposure.tif')
+        hra._calc_criteria(attributes_list, habitat_mask_path,
+                           target_exposure_path, decayed_distance_raster_path)
+
+        exposure_array = pygeoprocessing.raster_to_numpy_array(
+            target_exposure_path)
+        nodata = hra._TARGET_NODATA_FLOAT32
+        # These expected values were calculated by hand.
+        expected_exposure_array = numpy.array([
+            [nodata, 0, 0]], dtype=numpy.float32)
+        numpy.testing.assert_allclose(
+            exposure_array, expected_exposure_array)
+
     def test_decayed_distance_linear(self):
         """HRA: linear decay over a distance."""
         from natcap.invest import hra