Repository: natcap/invest (python)

Issue:
Error in "build binaries" job on macOS
"Build binaries" job is failing with this error when it tries to run `invest list` on the built binaries.
```
subprocess.CalledProcessError: Command '['/Users/runner/work/invest/invest/dist/invest/invest', '-sS', '-c', 'import platform; print(platform.mac_ver()[0])']' returned non-zero exit status 2.
```

Pull request description:
restrict scipy version to one that does not use `packaging.tags`
I'm not sure if we want to do this, but it does seem to allow our mac binary builds to work again. The change in scipy is probably this one in `_testutils.py`:

https://github.com/scipy/scipy/commit/96081aa10e929e807b199f8e441dd1693e304f03

Fixes #1508 

## Checklist
- [ ] 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 0853f86a3..8a76e41d4 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -52,15 +52,6 @@ Unreleased Changes
       administrative units vector now has the correct values for these fields,
       consistent with the user's guide chapter.
       https://github.com/natcap/invest/issues/1512
-    * Fixed an issue where certain nodata values were not being handled
-      correctly, leading to pixel values of +/- infinity in the urban nature
-      balance output raster.  https://github.com/natcap/invest/issues/1519
-* SDR
-    * Fixed an issue encountered in the sediment deposition function where
-      rasters with more than 2^32 pixels would raise a cryptic error relating
-      to negative dimensions. https://github.com/natcap/invest/issues/1431
-    * Optimized the creation of the summary vector by minimizing the number of
-      times the target vector needs to be rasterized.
 
 3.14.1 (2023-12-18)
 -------------------
diff --git a/requirements.txt b/requirements.txt
index 5b2583b0f..cad4fa0e7 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -16,7 +16,7 @@ pandas>=1.2.1
 numpy>=1.11.0,!=1.16.0
 Rtree>=0.8.2,!=0.9.1
 shapely>=2.0.0
-scipy>=1.9.0
+scipy>=1.9.0,!=1.12.*
 pygeoprocessing>=2.4.2  # pip-only
 taskgraph>=0.11.0
 psutil>=5.6.6
diff --git a/src/natcap/invest/sdr/sdr.py b/src/natcap/invest/sdr/sdr.py
index f1540982b..f23526b4a 100644
--- a/src/natcap/invest/sdr/sdr.py
+++ b/src/natcap/invest/sdr/sdr.py
@@ -19,7 +19,6 @@ from osgeo import ogr
 
 from .. import gettext
 from .. import spec_utils
-from .. import urban_nature_access
 from .. import utils
 from .. import validation
 from ..model_metadata import MODEL_METADATA
@@ -610,9 +609,7 @@ def execute(args):
     target_pixel_size = (min_pixel_size, -min_pixel_size)
 
     target_sr_wkt = dem_raster_info['projection_wkt']
-    vector_mask_options = {
-        'mask_vector_path': args['watersheds_path'],
-    }
+    vector_mask_options = {'mask_vector_path': args['watersheds_path']}
     align_task = task_graph.add_task(
         func=pygeoprocessing.align_and_resize_raster_stack,
         args=(
@@ -1494,26 +1491,18 @@ def _generate_report(
     target_layer = target_vector.GetLayer()
     target_layer.SyncToDisk()
 
-    # It's worth it to check if the geometries don't significantly overlap.
-    # On large rasters, this can save a TON of time rasterizing even a
-    # relatively simple vector.
-    geometries_might_overlap = urban_nature_access._geometries_overlap(
-        watershed_results_sdr_path)
-    fields_and_rasters = [
-        ('usle_tot', usle_path), ('sed_export', sed_export_path),
-        ('sed_dep', sed_deposition_path), ('avoid_exp', avoided_export_path),
-        ('avoid_eros', avoided_erosion_path)]
-
-    # Using the list option for raster path bands so that we can reduce
-    # rasterizations, which are costly on large datasets.
-    zonal_stats_results = pygeoprocessing.zonal_statistics(
-        [(raster_path, 1) for (_, raster_path) in fields_and_rasters],
-        watershed_results_sdr_path,
-        polygons_might_overlap=geometries_might_overlap)
-
     field_summaries = {
-        field: stats for ((field, _), stats) in
-        zip(fields_and_rasters, zonal_stats_results)}
+        'usle_tot': pygeoprocessing.zonal_statistics(
+            (usle_path, 1), watershed_results_sdr_path),
+        'sed_export': pygeoprocessing.zonal_statistics(
+            (sed_export_path, 1), watershed_results_sdr_path),
+        'sed_dep': pygeoprocessing.zonal_statistics(
+            (sed_deposition_path, 1), watershed_results_sdr_path),
+        'avoid_exp': pygeoprocessing.zonal_statistics(
+            (avoided_export_path, 1), watershed_results_sdr_path),
+        'avoid_eros': pygeoprocessing.zonal_statistics(
+            (avoided_erosion_path, 1), watershed_results_sdr_path),
+    }
 
     for field_name in field_summaries:
         field_def = ogr.FieldDefn(field_name, ogr.OFTReal)
diff --git a/src/natcap/invest/sdr/sdr_core.pyx b/src/natcap/invest/sdr/sdr_core.pyx
index 45b71b25b..c6848bfcf 100644
--- a/src/natcap/invest/sdr/sdr_core.pyx
+++ b/src/natcap/invest/sdr/sdr_core.pyx
@@ -436,17 +436,17 @@ def calculate_sediment_deposition(
     flow_dir_info = pygeoprocessing.get_raster_info(mfd_flow_direction_path)
     n_cols, n_rows = flow_dir_info['raster_size']
     cdef int mfd_nodata = 0
-    cdef stack[long] processing_stack
+    cdef stack[int] processing_stack
     cdef float sdr_nodata = pygeoprocessing.get_raster_info(
         sdr_path)['nodata'][0]
     cdef float e_prime_nodata = pygeoprocessing.get_raster_info(
         e_prime_path)['nodata'][0]
     cdef long col_index, row_index, win_xsize, win_ysize, xoff, yoff
     cdef long global_col, global_row, j, k
-    cdef long flat_index
+    cdef unsigned long flat_index
     cdef long seed_col = 0
     cdef long seed_row = 0
-    cdef long neighbor_row, neighbor_col, ds_neighbor_row, ds_neighbor_col
+    cdef long neighbor_row, neighbor_col
     cdef int flow_val, neighbor_flow_val, ds_neighbor_flow_val
     cdef int flow_weight, neighbor_flow_weight
     cdef float flow_sum, neighbor_flow_sum
diff --git a/src/natcap/invest/urban_nature_access.py b/src/natcap/invest/urban_nature_access.py
index b2b8a67ca..da083bbbb 100644
--- a/src/natcap/invest/urban_nature_access.py
+++ b/src/natcap/invest/urban_nature_access.py
@@ -1328,15 +1328,20 @@ def execute(args):
                 output_dir,
                 f'urban_nature_balance_percapita_{pop_group}{suffix}.tif')
             per_cap_urban_nature_balance_pop_group_task = graph.add_task(
-                _calculate_urban_nature_balance_percapita,
+                pygeoprocessing.raster_calculator,
                 kwargs={
-                    'urban_nature_supply_path':
-                        urban_nature_supply_percapita_to_group_path,
-                    'urban_nature_demand': float(args['urban_nature_demand']),
-                    'target_path':
-                        per_cap_urban_nature_balance_pop_group_path},
+                    'base_raster_path_band_const_list': [
+                        (urban_nature_supply_percapita_to_group_path, 1),
+                        (float(args['urban_nature_demand']), 'raw')
+                    ],
+                    'local_op': _urban_nature_balance_percapita_op,
+                    'target_raster_path':
+                        per_cap_urban_nature_balance_pop_group_path,
+                    'datatype_target': gdal.GDT_Float32,
+                    'nodata_target': FLOAT32_NODATA
+                },
                 task_name=(
-                    f'Calculate per-capita urban nature balance-{pop_group}'),
+                    f'Calculate per-capita urban nature balance - {pop_group}'),
                 target_path_list=[
                     per_cap_urban_nature_balance_pop_group_path],
                 dependent_task_list=[
@@ -1414,17 +1419,20 @@ def execute(args):
             ])
 
         per_capita_urban_nature_balance_task = graph.add_task(
-            _calculate_urban_nature_balance_percapita,
+            pygeoprocessing.raster_calculator,
             kwargs={
-                'urban_nature_supply_path':
-                    file_registry['urban_nature_supply_percapita'],
-                'urban_nature_demand': float(args['urban_nature_demand']),
-                'target_path':
-                    file_registry['urban_nature_balance_percapita']},
-            task_name=(
-                'Calculate per-capita urban nature balance}'),
-            target_path_list=[
-                file_registry['urban_nature_balance_percapita']],
+                'base_raster_path_band_const_list': [
+                    (file_registry['urban_nature_supply_percapita'], 1),
+                    (float(args['urban_nature_demand']), 'raw')
+                ],
+                'local_op': _urban_nature_balance_percapita_op,
+                'target_raster_path':
+                    file_registry['urban_nature_balance_percapita'],
+                'datatype_target': gdal.GDT_Float32,
+                'nodata_target': FLOAT32_NODATA
+            },
+            task_name='Calculate per-capita urban nature balance',
+            target_path_list=[file_registry['urban_nature_balance_percapita']],
             dependent_task_list=[
                 urban_nature_supply_percapita_task,
             ])
@@ -1471,17 +1479,20 @@ def execute(args):
                                       RADIUS_OPT_URBAN_NATURE):
         # This is "SUP_DEMi_cap" from the user's guide
         per_capita_urban_nature_balance_task = graph.add_task(
-            _calculate_urban_nature_balance_percapita,
+            pygeoprocessing.raster_calculator,
             kwargs={
-                'urban_nature_supply_path':
-                    file_registry['urban_nature_supply_percapita'],
-                'urban_nature_demand': float(args['urban_nature_demand']),
-                'target_path':
-                    file_registry['urban_nature_balance_percapita']},
-            task_name=(
-                'Calculate per-capita urban nature balance'),
-            target_path_list=[
-                file_registry['urban_nature_balance_percapita']],
+                'base_raster_path_band_const_list': [
+                    (file_registry['urban_nature_supply_percapita'], 1),
+                    (float(args['urban_nature_demand']), 'raw')
+                ],
+                'local_op': _urban_nature_balance_percapita_op,
+                'target_raster_path':
+                    file_registry['urban_nature_balance_percapita'],
+                'datatype_target': gdal.GDT_Float32,
+                'nodata_target': FLOAT32_NODATA
+            },
+            task_name='Calculate per-capita urban nature balance',
+            target_path_list=[file_registry['urban_nature_balance_percapita']],
             dependent_task_list=[
                 urban_nature_supply_percapita_task,
             ])
@@ -2125,44 +2136,28 @@ def _write_supply_demand_vector(source_aoi_vector_path, feature_attrs,
     target_vector = None
 
 
-def _calculate_urban_nature_balance_percapita(
-        urban_nature_supply_path, urban_nature_demand, target_path):
-    supply_nodata = pygeoprocessing.get_raster_info(
-        urban_nature_supply_path)['nodata'][0]
+def _urban_nature_balance_percapita_op(urban_nature_supply, urban_nature_demand):
+    """Calculate the per-capita urban nature balance.
 
-    def _urban_nature_balance_percapita_op(urban_nature_supply,
-                                           urban_nature_demand):
-        """Calculate the per-capita urban nature balance.
-
-        This is the amount of urban nature that each pixel has above (positive
-        values) or below (negative values) the user-defined
-        ``urban_nature_demand`` value.
-
-        Args:
-            urban_nature_supply (numpy.array): The supply of urban nature
-                available to each person in the population.  This is ``Ai`` in
-                the User's Guide.
-            urban_nature_demand (float): The policy-defined urban nature
-            requirement, in square meters per person.
+    This is the amount of urban nature that each pixel has above (positive
+    values) or below (negative values) the user-defined ``urban_nature_demand``
+    value.
 
-        Returns:
-            A ``numpy.array`` of the calculated urban nature budget.
-        """
-        balance = numpy.full(
-            urban_nature_supply.shape, FLOAT32_NODATA, dtype=numpy.float32)
-        valid_pixels = ~pygeoprocessing.array_equals_nodata(
-            urban_nature_supply, supply_nodata)
-        balance[valid_pixels] = (
-            urban_nature_supply[valid_pixels] - urban_nature_demand)
-        return balance
+    Args:
+        urban_nature_supply (numpy.array): The supply of urban nature available to
+            each person in the population.  This is ``Ai`` in the User's Guide.
+            This matrix must have ``FLOAT32_NODATA`` as its nodata value.
+        urban_nature_demand (float): The policy-defined urban nature requirement,
+            in square meters per person.
 
-    pygeoprocessing.raster_calculator(
-        base_raster_path_band_const_list=[
-            (urban_nature_supply_path, 1), (urban_nature_demand, 'raw')],
-        local_op=_urban_nature_balance_percapita_op,
-        target_raster_path=target_path,
-        datatype_target=gdal.GDT_Float32,
-        nodata_target=FLOAT32_NODATA)
+    Returns:
+        A ``numpy.array`` of the calculated urban nature budget.
+    """
+    balance = numpy.full(
+        urban_nature_supply.shape, FLOAT32_NODATA, dtype=numpy.float32)
+    valid_pixels = ~numpy.isclose(urban_nature_supply, FLOAT32_NODATA)
+    balance[valid_pixels] = urban_nature_supply[valid_pixels] - urban_nature_demand
+    return balance
 
 
 def _urban_nature_balance_totalpop_op(urban_nature_balance, population):


Test patch:
diff --git a/tests/test_urban_nature_access.py b/tests/test_urban_nature_access.py
index e09a158f7..e0791022a 100644
--- a/tests/test_urban_nature_access.py
+++ b/tests/test_urban_nature_access.py
@@ -217,19 +217,14 @@ class UNATests(unittest.TestCase):
             [nodata, 100.5],
             [75, 100]], dtype=numpy.float32)
         urban_nature_demand = 50
-        supply_path = os.path.join(self.workspace_dir, 'supply.path')
-        target_path = os.path.join(self.workspace_dir, 'target.path')
 
-        pygeoprocessing.numpy_array_to_raster(
-            urban_nature_supply_percapita, nodata, _DEFAULT_PIXEL_SIZE,
-            _DEFAULT_ORIGIN, _DEFAULT_SRS.ExportToWkt(), supply_path)
-
-        urban_nature_access._calculate_urban_nature_balance_percapita(
-            supply_path, urban_nature_demand, target_path)
-
-        urban_nature_budget = pygeoprocessing.raster_to_numpy_array(
-            target_path)
+        population = numpy.array([
+            [50, 100],
+            [40.75, nodata]], dtype=numpy.float32)
 
+        urban_nature_budget = (
+            urban_nature_access._urban_nature_balance_percapita_op(
+                urban_nature_supply_percapita, urban_nature_demand))
         expected_urban_nature_budget = numpy.array([
             [nodata, 50.5],
             [25, 50]], dtype=numpy.float32)