Repository: natcap/invest (python)

Issue:
SWY monthly QF outputs have 0 values instead of nodata
@newtpatrol was testing 3.14.2 and noticed that `QF.tif` had some 0 values around the watershed, where it looks like `nodata` values should be. I found that the monthly QF is where this is introduced. I'm not yet sure if it's part of the logic of the function `_calculate_monthly_quickflow_flow` or if it's a nodata masking bug.

![image](https://github.com/natcap/invest/assets/2659980/ba7eaadf-b9f6-4a6c-98c5-6b790a40b7de)


Pull request description:
Fix nodata handling in monthly QF output for SWY model
## Description
This PR fixes a bug in `_calculate_monthly_quick_flow` where pixels with nodata values in precipitation or number of rain events were incorrectly assigned a value of 0 in the QF output raster. Now, such pixels are correctly masked as nodata in the output. I also changed two variable names so that masking could be more easily interpreted and added a test to confirm that nodata values are treated correctly. 

Fixes #1592

## 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 Workbench UI (if relevant)


Patch:
diff --git a/HISTORY.rst b/HISTORY.rst
index 7fee4b929..699eb11dc 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -87,6 +87,12 @@ Crop Production
   of a corresponding regression yield table
   (`#1723 <https://github.com/natcap/invest/issues/1723>`_).
 
+Seasonal Water Yield
+====================
+* Fixed bug in quickflow calculation where `NoData` values in precipitation
+  were being incorrectly converted to 0's
+  (`#1592 <https://github.com/natcap/invest/issues/1592>`_).
+
 
 3.15.0 (2025-04-03)
 -------------------
diff --git a/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py b/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
index 8cc617e67..b283e3c6c 100644
--- a/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
+++ b/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
@@ -1097,8 +1097,8 @@ def _calculate_monthly_quick_flow(precip_path, n_events_path, stream_path,
         valid_p_mask = ~pygeoprocessing.array_equals_nodata(p_im, p_nodata)
         valid_n_mask = ~pygeoprocessing.array_equals_nodata(n_m, n_nodata)
         # precip mask: both p_im and n_m are defined and greater than 0
-        precip_mask = valid_p_mask & valid_n_mask & (p_im > 0) & (n_m > 0)
-        stream_mask = stream == 1
+        precip_ok = valid_p_mask & valid_n_mask & (p_im > 0) & (n_m > 0)
+        on_stream = stream == 1
         # stream_nodata is the only input that carries over nodata values from
         # the aligned DEM.
         valid_mask = (
@@ -1109,7 +1109,7 @@ def _calculate_monthly_quick_flow(precip_path, n_events_path, stream_path,
 
         # QF is defined in terms of three cases:
         #
-        # 1. Where there is no precipitation, QF = 0
+        # 1. Where precipitation <=0, QF = 0
         #    (even if stream or s_i is undefined)
         #
         # 2. Where there is precipitation and we're on a stream, QF = P
@@ -1163,14 +1163,15 @@ def _calculate_monthly_quick_flow(precip_path, n_events_path, stream_path,
         # qf_im is the quickflow at pixel i on month m
         qf_im = numpy.full(p_im.shape, TARGET_NODATA, dtype=numpy.float32)
 
-        # case 1: where there is no precipitation
-        qf_im[~precip_mask] = 0
+        # case 1: where precip or n are <=0 but not nodata
+        case_1_mask = ~precip_ok & valid_p_mask & valid_n_mask
+        qf_im[case_1_mask] = 0
 
         # case 2: where there is precipitation and we're on a stream
-        qf_im[precip_mask & stream_mask] = p_im[precip_mask & stream_mask]
+        qf_im[precip_ok & on_stream] = p_im[precip_ok & on_stream]
 
         # case 3: where there is precipitation and we're not on a stream
-        case_3_mask = valid_mask & precip_mask & ~stream_mask
+        case_3_mask = valid_mask & precip_ok & ~on_stream
 
         # for consistent indexing, make a_im the same shape as the other
         # arrays even though we only use it in case 3


Test patch:
diff --git a/tests/test_seasonal_water_yield_regression.py b/tests/test_seasonal_water_yield_regression.py
index 618f49cac..5a06384c5 100644
--- a/tests/test_seasonal_water_yield_regression.py
+++ b/tests/test_seasonal_water_yield_regression.py
@@ -61,7 +61,7 @@ def make_simple_shp(base_shp_path, origin):
     data_source = None
 
 
-def make_raster_from_array(base_array, base_raster_path):
+def make_raster_from_array(base_array, base_raster_path, nodata=-1):
     """Make a raster from an array on a designated path.
 
     Args:
@@ -78,7 +78,7 @@ def make_raster_from_array(base_array, base_raster_path):
 
     # Each pixel is 1x1 m
     pygeoprocessing.numpy_array_to_raster(
-        base_array, -1, (1, -1), (1180000, 690000), project_wkt,
+        base_array, nodata, (1, -1), (1180000, 690000), project_wkt,
         base_raster_path)
 
 
@@ -1263,6 +1263,52 @@ class SeasonalWaterYieldRegressionTests(unittest.TestCase):
             pygeoprocessing.raster_to_numpy_array(output_path),
             expected_quickflow_array, atol=1e-5)
 
+    def test_monthly_quickflow_nodata_propagation(self):
+        """Test correct nodata propagation in `_calculate_monthly_quick_flow`
+
+        This test checks that:
+        1. If n=nodata: output is nodata
+        2. If precip=nodata: output is nodata
+        3. If precip<0 and not nodata & n is valid: output is 0
+        4. If precip and n are valid & stream=1 & SI=nodata: output is valid
+        5. If precip and n are valid & stream=nodata: output is nodata
+        """
+        from natcap.invest.seasonal_water_yield import seasonal_water_yield
+
+        # Test a variety of valid/nodata combinations across the input layers
+        precip_array = numpy.array([[-1, -6, 32767, 32767],
+                                    [5, 6, 30, 8]], dtype=numpy.float32)
+        n_events_array = numpy.array([[-1, 1, -8, 8],
+                                      [-1, 6, 2, 9]], dtype=numpy.float32)
+        si_array = numpy.array([[1, -1, 3, 4],
+                                [5, -1, 7, 8]], dtype=numpy.float32)
+        stream_mask = numpy.array([[1, -1, 1, 1],
+                                   [1, 1, 0, -1]], dtype=numpy.float32)
+        expected_quickflow_array = numpy.array([[-1, 0, -1, -1],
+                                                [-1, 6, 0.382035, -1]])
+
+        precip_path = os.path.join(self.workspace_dir, 'precip.tif')
+        si_path = os.path.join(self.workspace_dir, 'si.tif')
+        n_events_path = os.path.join(self.workspace_dir, 'n_events.tif')
+        stream_path = os.path.join(self.workspace_dir, 'stream.tif')
+        output_path = os.path.join(self.workspace_dir, 'quickflow.tif')
+
+        # write all the test arrays to raster files
+        for array, path in [(n_events_array, n_events_path),
+                            (si_array, si_path),
+                            (stream_mask, stream_path)]:
+            # define a nodata value for intermediate outputs
+            make_raster_from_array(array, path)
+
+        # Ensure positive nodata value for precip is handled correctly
+        make_raster_from_array(precip_array, precip_path, nodata=32767)
+
+        seasonal_water_yield._calculate_monthly_quick_flow(
+            precip_path, n_events_path, stream_path, si_path, output_path)
+        numpy.testing.assert_allclose(
+            pygeoprocessing.raster_to_numpy_array(output_path),
+            expected_quickflow_array, atol=1e-6)
+
     def test_local_recharge_undefined_nodata(self):
         """Test `calculate_local_recharge` with undefined nodata values"""
         from natcap.invest.seasonal_water_yield import \