Merge pull request #1520 from phargogh/bugfix/1519-inf-values-in-una-…

…output

Fix nodata value masking issue in UNA's urban nature balance

HISTORY.rst

@@ -52,6 +52,9 @@ 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
 
 3.14.1 (2023-12-18)
 -------------------

src/natcap/invest/urban_nature_access.py

@@ -1328,20 +1328,15 @@ def decay_func(dist_array):
                 output_dir,
                 f'urban_nature_balance_percapita_{pop_group}{suffix}.tif')
             per_cap_urban_nature_balance_pop_group_task = graph.add_task(
-                pygeoprocessing.raster_calculator,
+                _calculate_urban_nature_balance_percapita,
                 kwargs={
-                    '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
-                },
+                    '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},
                 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=[
@@ -1419,20 +1414,17 @@ def decay_func(dist_array):
             ])
 
         per_capita_urban_nature_balance_task = graph.add_task(
-            pygeoprocessing.raster_calculator,
+            _calculate_urban_nature_balance_percapita,
             kwargs={
-                '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']],
+                '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']],
             dependent_task_list=[
                 urban_nature_supply_percapita_task,
             ])
@@ -1479,20 +1471,17 @@ def decay_func(dist_array):
                                       RADIUS_OPT_URBAN_NATURE):
         # This is "SUP_DEMi_cap" from the user's guide
         per_capita_urban_nature_balance_task = graph.add_task(
-            pygeoprocessing.raster_calculator,
+            _calculate_urban_nature_balance_percapita,
             kwargs={
-                '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']],
+                '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']],
             dependent_task_list=[
                 urban_nature_supply_percapita_task,
             ])
@@ -2136,28 +2125,44 @@ def _write_supply_demand_vector(source_aoi_vector_path, feature_attrs,
     target_vector = None
 
 
-def _urban_nature_balance_percapita_op(urban_nature_supply, urban_nature_demand):
-    """Calculate the per-capita urban nature balance.
+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]
 
-    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.
+    def _urban_nature_balance_percapita_op(urban_nature_supply,
+                                           urban_nature_demand):
+        """Calculate the per-capita urban nature 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.
+        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 = ~numpy.isclose(urban_nature_supply, FLOAT32_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.
+            urban_nature_demand (float): The policy-defined urban nature
+            requirement, in square meters per person.
+
+        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
+
+    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)
 
 
 def _urban_nature_balance_totalpop_op(urban_nature_balance, population):

tests/test_urban_nature_access.py

@@ -217,14 +217,19 @@ def test_urban_nature_balance(self):
             [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')
 
-        population = numpy.array([
-            [50, 100],
-            [40.75, nodata]], dtype=numpy.float32)
+        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)
 
-        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)