Repository: natcap/invest (python)

Issue:
NDR runoff proxy index depends on watershed extent
@phawthorne raised this issue:
The runoff proxy index is calculated as $RPI_i = RP_i / RP_{av}$ (dividing runoff proxy by average runoff proxy within the watersheds). This is problematic because average runoff proxy $RP_{av}$ changes depending on which pixels/watersheds are included in the raster. @phawthorne said they have even considered inserting fake pixel values into the RP raster in order to change the $RP_{av}$. We could expose the $RP_{av}$ parameter and allow the user to specify it instead.

Pull request description:
Expose runoff proxy average as a user-specified parameter in NDR RPI calculation
## Description
Allows a user to specify their own average runoff proxy value. More specifically, this PR modifies the calculation of the runoff proxy index (RPI) in the NDR model by exposing the average runoff proxy as an optional user-defined parameter. Previously, RP average was dynamically computed based on the included watershed pixels, leading to inconsistencies when model inputs changed (e.g., different watershed extents or climate scenarios). By allowing users to specify a runoff proxy average, this change ensures consistent RPI values across scenarios, enabling better comparisons between baseline and future conditions.

Changes:

- Introduces new optional model parameter `Runoff Proxy Average`
- Updates RPI calculation to use the user-defined average runoff proxy (if specified) or the automatically, dynamically calculated average if not specified by the user
- Adds a test and updates documentation

Fixes #1741

UG PR: https://github.com/natcap/invest.users-guide/pull/162

## Checklist
- [x] Updated HISTORY.rst and link to any relevant issue (if these changes are user-facing)
- [x] Updated the user's guide (if needed)
- [x] Tested the Workbench UI (if relevant)


Patch:
diff --git a/HISTORY.rst b/HISTORY.rst
index 2a1639d56..1b5e4765a 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -100,6 +100,10 @@ Unreleased Changes
     * Raster outputs that previously contained per-pixel values (e.g., kg/pixel)
       now contain per-hectare values (e.g., kg/ha). (`InVEST #1270
       <https://github.com/natcap/invest/issues/1270>`_).
+    * Made the runoff proxy index calculation more robust by allowing users to
+      specify the average runoff proxy, preventing normalization issues across
+      different climate scenarios and watershed selections.
+      https://github.com/natcap/invest/issues/1741
 * SDR
     * Raster outputs that previously contained per-pixel values (e.g., t/pixel)
       now contain per-hectare values (e.g., t/ha). (`InVEST #1270
diff --git a/src/natcap/invest/ndr/ndr.py b/src/natcap/invest/ndr/ndr.py
index f5e406675..802f882ee 100644
--- a/src/natcap/invest/ndr/ndr.py
+++ b/src/natcap/invest/ndr/ndr.py
@@ -140,6 +140,21 @@ MODEL_SPEC = {
                 "actually reaches the stream)."),
             "name": gettext("Borselli k parameter"),
         },
+        "runoff_proxy_av": {
+            "type": "number",
+            "units": u.none,
+            'expression': 'value > 0',
+            "required": False,
+            "name": gettext("average runoff proxy"),
+            "about": gettext(
+                "This parameter allows the user to specify a predefined "
+                "average value for the runoff proxy. This value is used "
+                "to normalize the Runoff Proxy raster when calculating "
+                "the Runoff Proxy Index (RPI). If a user does not specify "
+                "the runoff proxy average, this value will be automatically "
+                "calculated from the Runoff Proxy raster. The units will "
+                "be the same as those in the Runoff Proxy raster."),
+        },
         "subsurface_critical_length_n": {
             "type": "number",
             "units": u.meter,
@@ -522,6 +537,9 @@ def execute(args):
         args['k_param'] (number): The Borselli k parameter. This is a
             calibration parameter that determines the shape of the
             relationship between hydrologic connectivity.
+        args['runoff_proxy_av'] (number): (optional) The average runoff proxy.
+            Used to calculate the runoff proxy index. If not specified,
+            it will be automatically calculated.
         args['subsurface_critical_length_n'] (number): The distance (traveled
             subsurface and downslope) after which it is assumed that soil
             retains nutrient at its maximum capacity, given in meters. If
@@ -576,6 +594,11 @@ def execute(args):
         args['biophysical_table_path'],
         **MODEL_SPEC['args']['biophysical_table_path'])
 
+    # Ensure that if user doesn't explicitly assign a value,
+    # runoff_proxy_av = None
+    runoff_proxy_av = args.get("runoff_proxy_av")
+    runoff_proxy_av = float(runoff_proxy_av) if runoff_proxy_av else None
+
     # these are used for aggregation in the last step
     field_pickle_map = {}
 
@@ -718,6 +741,7 @@ def execute(args):
         func=_normalize_raster,
         args=((f_reg['masked_runoff_proxy_path'], 1),
               f_reg['runoff_proxy_index_path']),
+        kwargs={'user_provided_mean': runoff_proxy_av},
         target_path_list=[f_reg['runoff_proxy_index_path']],
         dependent_task_list=[align_raster_task, mask_runoff_proxy_task],
         task_name='runoff proxy mean')
@@ -1198,7 +1222,8 @@ def validate(args, limit_to=None):
     return validation_warnings
 
 
-def _normalize_raster(base_raster_path_band, target_normalized_raster_path):
+def _normalize_raster(base_raster_path_band, target_normalized_raster_path,
+                      user_provided_mean=None):
     """Calculate normalize raster by dividing by the mean value.
 
     Args:
@@ -1206,20 +1231,27 @@ def _normalize_raster(base_raster_path_band, target_normalized_raster_path):
             mean.
         target_normalized_raster_path (string): path to target normalized
             raster from base_raster_path_band.
+        user_provided_mean (float, optional): user-provided average.
+            If provided, this value will be used instead of computing
+            the mean from the raster.
 
     Returns:
         None.
 
     """
-    value_sum, value_count = pygeoprocessing.raster_reduce(
-        function=lambda sum_count, block: (  # calculate both in one pass
-            sum_count[0] + numpy.sum(block), sum_count[1] + block.size),
-        raster_path_band=base_raster_path_band,
-        initializer=(0, 0))
-
-    value_mean = value_sum
-    if value_count > 0:
-        value_mean /= value_count
+    if user_provided_mean is None:
+        value_sum, value_count = pygeoprocessing.raster_reduce(
+            function=lambda sum_count, block: (  # calculate both in one pass
+                sum_count[0] + numpy.sum(block), sum_count[1] + block.size),
+            raster_path_band=base_raster_path_band,
+            initializer=(0, 0))
+        value_mean = value_sum
+        if value_count > 0:
+            value_mean /= value_count
+        LOGGER.info(f"Normalizing raster ({base_raster_path_band[0]}) using "
+                    f"auto-calculated mean: {value_mean}")
+    else:
+        value_mean = user_provided_mean
 
     pygeoprocessing.raster_map(
         op=lambda array: array if value_mean == 0 else array / value_mean,
diff --git a/workbench/src/renderer/ui_config.js b/workbench/src/renderer/ui_config.js
index b9a7af3d6..9ae3a6970 100644
--- a/workbench/src/renderer/ui_config.js
+++ b/workbench/src/renderer/ui_config.js
@@ -202,7 +202,7 @@ const UI_SPEC = {
       ['dem_path', 'lulc_path', 'runoff_proxy_path', 'watersheds_path', 'biophysical_table_path'],
       ['calc_p'],
       ['calc_n', 'subsurface_critical_length_n', 'subsurface_eff_n'],
-      ['threshold_flow_accumulation', 'k_param'],
+      ['threshold_flow_accumulation', 'k_param', 'runoff_proxy_av'],
     ],
     enabledFunctions: {
       subsurface_critical_length_n: isSufficient.bind(null, 'calc_n'),


Test patch:
diff --git a/tests/test_ndr.py b/tests/test_ndr.py
index f70bdb9ca..e99d900a4 100644
--- a/tests/test_ndr.py
+++ b/tests/test_ndr.py
@@ -5,6 +5,7 @@ import tempfile
 import unittest
 
 import numpy
+import pandas
 import pygeoprocessing
 import shapely.geometry
 from osgeo import gdal
@@ -435,3 +436,48 @@ class NDRTests(unittest.TestCase):
         numpy.testing.assert_array_equal(
             expected_array,
             pygeoprocessing.raster_to_numpy_array(target_raster_path))
+
+    def test_synthetic_runoff_proxy_av(self):
+        """
+        Test RPI given user-entered or auto-calculated runoff proxy average.
+
+        Test that the runoff proxy index (RPI) is calculated correctly if
+        (1) the user specifies a runoff proxy average value,
+        (2) the user does not specify a value so the runoff proxy average
+            is auto-calculated.
+        """
+        from natcap.invest.ndr import ndr
+
+        # make simple raster
+        runoff_proxy_path = os.path.join(self.workspace_dir, "ppt.tif")
+        runoff_proxy_array = numpy.array(
+            [[800, 799, 567, 234], [765, 867, 765, 654]], dtype=numpy.float32)
+        srs = osr.SpatialReference()
+        srs.ImportFromEPSG(26910)
+        projection_wkt = srs.ExportToWkt()
+        origin = (461251, 4923445)
+        pixel_size = (30, -30)
+        no_data = -1
+        pygeoprocessing.numpy_array_to_raster(
+            runoff_proxy_array, no_data, pixel_size, origin, projection_wkt,
+            runoff_proxy_path)
+        target_rpi_path = os.path.join(self.workspace_dir, "out_raster.tif")
+
+        # Calculate RPI with user-specified runoff proxy average
+        runoff_proxy_av = 2
+        ndr._normalize_raster((runoff_proxy_path, 1), target_rpi_path,
+                              user_provided_mean=runoff_proxy_av)
+
+        actual_rpi = pygeoprocessing.raster_to_numpy_array(target_rpi_path)
+        expected_rpi = runoff_proxy_array/runoff_proxy_av
+
+        numpy.testing.assert_allclose(actual_rpi, expected_rpi)
+
+        # Now calculate RPI with auto-calculated RP average
+        ndr._normalize_raster((runoff_proxy_path, 1), target_rpi_path,
+                              user_provided_mean=None)
+
+        actual_rpi = pygeoprocessing.raster_to_numpy_array(target_rpi_path)
+        expected_rpi = runoff_proxy_array/numpy.mean(runoff_proxy_array)
+
+        numpy.testing.assert_allclose(actual_rpi, expected_rpi)