Repository: natcap/invest (python)

Issue:
Nodata casting error.
## Description
There is a potential mismatch between the nodata assigned to the numpy array and the numpy type at line 1187 (from task at line 672).

## To replicate

1. Use a LULC raster of type "Byte - Eight bit unsigned integer" with "n/a" nodata (as reported by QGIS).
2. Run NDR model in Python - I could not replicate this behavior with the Workbench

## Expected behavior
I expect the model to run! More specifically, I expect the LULC raster to be successfully masked by the watersheds.

## Observed behavior
OverflowError: Python integer -2147483648 out of bounds for uint8

## Logfile(s)
No log file output from python run.

## Screenshots
If applicable, add screenshots to help illustrate the problem.

## Environment
- Operating System: macOS Sequoia 15.4.1
- InVEST version number: 3.14.3 (although the code seems to be unchanged in 3.15.0)
- InVEST interface (Workbench, CLI, or Python API): Python API

## Additional context
I think this is a fairly straightforward fix. When passing a `source raster` of Int8 with undefined nodata and defining `default_nodata` as `numpy.iinfo(numpy.int32).min` in the `_mask_raster()` function (line 1160), the function tries to create a Byte array full of -2147483648 values, which breaks numpy.
This happens in particular with `mask_lulc_task`, line ~ 672.

<img width="1670" alt="Image" src="https://github.com/user-attachments/assets/0c878969-5f41-4cbe-a9e1-51b182d0b09e" />


Pull request:
Bugfix for NDR: OverflowError when filling nodata
`_mask_raster` was trying to fill a numpy array with a module-defined (as opposed to user-defined) default nodata value, but using the source raster's native dtype for the array.

This only occurred when the source raster didn't have it's own nodata value defined. Depending on the source raster's datatype, the default nodata value might overflow it.

Since `_mask_raster` is always specifying the target dtype for the raster, it didn't make sense to use the source raster's dtype to initialize the numpy array.

Also, since we're always specifying the target dtype, I didn't see a reason to try to re-use the source raster's nodata value. So I changed the parameter from `default_nodata` to `target_nodata`.

Fixes #1904 

## 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 d0fa4a467..63be9207d 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -86,6 +86,8 @@ NDR
 ====
 * ``stream.tif`` is now saved in the main output folder rather than the
   intermediate folder (`#1864 <https://github.com/natcap/invest/issues/1864>`_).
+* Fixed a bug where input rasters (e.g. LULC) without a defined nodata value could
+  cause an OverflowError. (`#1904 <https://github.com/natcap/invest/issues/1904>`_).
 
 Seasonal Water Yield
 ====================
diff --git a/src/natcap/invest/ndr/ndr.py b/src/natcap/invest/ndr/ndr.py
index 839db3f0c..11f5c875b 100644
--- a/src/natcap/invest/ndr/ndr.py
+++ b/src/natcap/invest/ndr/ndr.py
@@ -9,6 +9,7 @@ import pygeoprocessing
 import pygeoprocessing.routing
 import taskgraph
 from osgeo import gdal
+from osgeo import gdal_array
 from osgeo import ogr
 
 from .. import gettext
@@ -650,7 +651,7 @@ def execute(args):
             'mask_raster_path': f_reg['mask_path'],
             'target_masked_raster_path': f_reg['masked_runoff_proxy_path'],
             'target_dtype': gdal.GDT_Float32,
-            'default_nodata': _TARGET_NODATA,
+            'target_nodata': _TARGET_NODATA,
         },
         dependent_task_list=[mask_task, align_raster_task],
         target_path_list=[f_reg['masked_runoff_proxy_path']],
@@ -663,7 +664,7 @@ def execute(args):
             'mask_raster_path': f_reg['mask_path'],
             'target_masked_raster_path': f_reg['masked_dem_path'],
             'target_dtype': gdal.GDT_Float32,
-            'default_nodata': float(numpy.finfo(numpy.float32).min),
+            'target_nodata': float(numpy.finfo(numpy.float32).min),
         },
         dependent_task_list=[mask_task, align_raster_task],
         target_path_list=[f_reg['masked_dem_path']],
@@ -676,7 +677,7 @@ def execute(args):
             'mask_raster_path': f_reg['mask_path'],
             'target_masked_raster_path': f_reg['masked_lulc_path'],
             'target_dtype': gdal.GDT_Int32,
-            'default_nodata': numpy.iinfo(numpy.int32).min,
+            'target_nodata': numpy.iinfo(numpy.int32).min,
         },
         dependent_task_list=[mask_task, align_raster_task],
         target_path_list=[f_reg['masked_lulc_path']],
@@ -1158,7 +1159,7 @@ def _create_mask_raster(source_raster_path, source_vector_path,
 
 
 def _mask_raster(source_raster_path, mask_raster_path,
-                 target_masked_raster_path, default_nodata, target_dtype):
+                 target_masked_raster_path, target_nodata, target_dtype):
     """Using a raster of 1s and 0s, determine which pixels remain in output.
 
     Args:
@@ -1170,8 +1171,8 @@ def _mask_raster(source_raster_path, mask_raster_path,
             target raster.
         target_masked_raster_path (str): The path to where the target raster
             should be written.
-        default_nodata (int, float, None): The nodata value that should be used
-            if ``source_raster_path`` does not have a defined nodata value.
+        target_nodata (int, float): The target nodata value that should match
+            ``target_dtype``.
         target_dtype (int): The ``gdal.GDT_*`` datatype of the target raster.
 
     Returns:
@@ -1179,22 +1180,20 @@ def _mask_raster(source_raster_path, mask_raster_path,
     """
     source_raster_info = pygeoprocessing.get_raster_info(source_raster_path)
     source_nodata = source_raster_info['nodata'][0]
-    nodata = source_nodata
-    if nodata is None:
-        nodata = default_nodata
+    target_numpy_dtype = gdal_array.GDALTypeCodeToNumericTypeCode(target_dtype)
 
     def _mask_op(mask, raster):
-        result = numpy.full(mask.shape, nodata,
-                            dtype=source_raster_info['numpy_type'])
+        result = numpy.full(mask.shape, target_nodata,
+                            dtype=target_numpy_dtype)
         valid_pixels = (
-            ~pygeoprocessing.array_equals_nodata(raster, nodata) &
+            ~pygeoprocessing.array_equals_nodata(raster, source_nodata) &
             (mask == 1))
         result[valid_pixels] = raster[valid_pixels]
         return result
 
     pygeoprocessing.raster_calculator(
         [(mask_raster_path, 1), (source_raster_path, 1)], _mask_op,
-        target_masked_raster_path, target_dtype, nodata)
+        target_masked_raster_path, target_dtype, target_nodata)
 
 
 def _add_fields_to_shapefile(field_pickle_map, target_vector_path):


Test patch:
diff --git a/tests/test_ndr.py b/tests/test_ndr.py
index bc73d58aa..f120aa5a5 100644
--- a/tests/test_ndr.py
+++ b/tests/test_ndr.py
@@ -376,6 +376,39 @@ class NDRTests(unittest.TestCase):
         if mismatch_list:
             raise RuntimeError("results not expected: %s" % mismatch_list)
 
+    def test_mask_raster_nodata_overflow(self):
+        """NDR test when target nodata value overflows source dtype."""
+        from natcap.invest.ndr import ndr
+
+        source_raster_path = os.path.join(self.workspace_dir, 'source.tif')
+        target_raster_path = os.path.join(
+            self.workspace_dir, 'target.tif')
+        source_dtype = numpy.int8
+        target_dtype = gdal.GDT_Int32
+        target_nodata = numpy.iinfo(numpy.int32).min
+
+        pygeoprocessing.numpy_array_to_raster(
+            base_array=numpy.full((4, 4), 1, dtype=source_dtype),
+            target_nodata=None,
+            pixel_size=(1, -1),
+            origin=(0, 0),
+            projection_wkt=None,
+            target_path=source_raster_path)
+
+        ndr._mask_raster(
+            source_raster_path=source_raster_path,
+            mask_raster_path=source_raster_path,  # mask=source for convenience
+            target_masked_raster_path=target_raster_path,
+            target_nodata=target_nodata,
+            target_dtype=target_dtype)
+
+        # Mostly we're testing that _mask_raster did not raise an OverflowError,
+        # but we can assert the results anyway.
+        array = pygeoprocessing.raster_to_numpy_array(target_raster_path)
+        numpy.testing.assert_array_equal(
+            array,
+            numpy.full((4, 4), 1, dtype=numpy.int32))  # matches target_dtype
+
     def test_validation(self):
         """NDR test argument validation."""
         from natcap.invest import validation