diff --git a/HISTORY.rst b/HISTORY.rst index 8a76e41d4..85d524f62 100644 --- a/HISTORY.rst +++ b/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) ------------------- diff --git a/src/natcap/invest/urban_nature_access.py b/src/natcap/invest/urban_nature_access.py index da083bbbb..b2b8a67ca 100644 --- a/src/natcap/invest/urban_nature_access.py +++ b/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): diff --git a/tests/test_urban_nature_access.py b/tests/test_urban_nature_access.py index e0791022a..e09a158f7 100644 --- a/tests/test_urban_nature_access.py +++ b/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)