From 4a114be5152a2e63ae6263dabfd3279aa27e29e3 Mon Sep 17 00:00:00 2001 From: Claire Simpson Date: Wed, 11 Mar 2026 16:07:26 -0600 Subject: [PATCH 1/3] Fix nan nodata handling in convolve2d #473 --- src/pygeoprocessing/geoprocessing.py | 12 ++++++------ tests/test_geoprocessing.py | 26 ++++++++++++++++++++++++++ 2 files changed, 32 insertions(+), 6 deletions(-) diff --git a/src/pygeoprocessing/geoprocessing.py b/src/pygeoprocessing/geoprocessing.py index 13d35edc..cd808c51 100644 --- a/src/pygeoprocessing/geoprocessing.py +++ b/src/pygeoprocessing/geoprocessing.py @@ -3335,7 +3335,7 @@ def convolve_2d( kernel_sum = 0.0 for _, kernel_block in iterblocks(kernel_path_band): if kernel_nodata is not None and ignore_nodata_and_edges: - kernel_block[numpy.isclose(kernel_block, kernel_nodata)] = 0.0 + kernel_block[array_equals_nodata(kernel_block, kernel_nodata)] = 0.0 kernel_sum += numpy.sum(kernel_block) # limit the size of the work queue since a large kernel / signal with small @@ -3426,7 +3426,7 @@ def _fill_work_queue(): # guard against a None nodata value if s_nodata is not None and mask_nodata: valid_mask[:] = ( - ~numpy.isclose(potential_nodata_signal_array, s_nodata)) + ~array_equals_nodata(potential_nodata_signal_array, s_nodata)) output_array[:] = target_nodata output_array[valid_mask] = ( (result[top_index_result:bottom_index_result, @@ -3474,7 +3474,7 @@ def _fill_work_queue(): signal_block = signal_band.ReadAsArray(**target_offset_data) mask_block = mask_band.ReadAsArray(**target_offset_data) if mask_nodata and signal_nodata is not None: - valid_mask = ~numpy.isclose(signal_block, signal_nodata) + valid_mask = ~array_equals_nodata(signal_block, signal_nodata) else: valid_mask = numpy.ones(target_block.shape, dtype=bool) valid_mask &= (mask_block > 0) @@ -4071,7 +4071,7 @@ def _convolve_2d_worker( kernel_sum = 0.0 for _, kernel_block in iterblocks(kernel_path_band): if kernel_nodata is not None and ignore_nodata: - kernel_block[numpy.isclose(kernel_block, kernel_nodata)] = 0.0 + kernel_block[array_equals_nodata(kernel_block, kernel_nodata)] = 0.0 kernel_sum += numpy.sum(kernel_block) while True: @@ -4090,7 +4090,7 @@ def _convolve_2d_worker( # don't ever convolve the nodata value if signal_nodata is not None: - signal_nodata_mask = numpy.isclose(signal_block, signal_nodata) + signal_nodata_mask = array_equals_nodata(signal_block, signal_nodata) signal_block[signal_nodata_mask] = 0.0 if not ignore_nodata: signal_nodata_mask[:] = 0 @@ -4122,7 +4122,7 @@ def _convolve_2d_worker( continue if kernel_nodata is not None: - kernel_block[numpy.isclose(kernel_block, kernel_nodata)] = 0.0 + kernel_block[array_equals_nodata(kernel_block, kernel_nodata)] = 0.0 if normalize_kernel: kernel_block /= kernel_sum diff --git a/tests/test_geoprocessing.py b/tests/test_geoprocessing.py index 0567ef12..f984f09b 100644 --- a/tests/test_geoprocessing.py +++ b/tests/test_geoprocessing.py @@ -3708,6 +3708,32 @@ def test_convolve_2d_ignore_undefined_nodata(self): numpy.isclose(signal_nodata_array, signal_nodata_none_array).all(), 'signal with nodata should be the same as signal with none') + def test_convolve_2d_nan_nodata(self): + """PGP.geoprocessing: test convolve_2d on raster w nan nodata value.""" + # Create signal array with one nan pixel + signal_array = numpy.arange(25).reshape(5, 5).astype(numpy.float32) + signal_array = numpy.ones((5, 5), numpy.float32) + signal_array[1, 1] = numpy.nan + signal_path = os.path.join(self.workspace_dir, 'signal.tif') + _array_to_raster(signal_array, numpy.nan, signal_path) + + kernel_array = numpy.array(((0, 1, 0), + (1, 1, 1), + (0, 1, 0)), numpy.float32) + kernel_path = os.path.join(self.workspace_dir, 'kernel.tif') + _array_to_raster(kernel_array, numpy.nan, kernel_path) + + target_path = os.path.join(self.workspace_dir, 'target.tif') + pygeoprocessing.convolve_2d( + (signal_path, 1), (kernel_path, 1), target_path, + ignore_nodata_and_edges=True, mask_nodata=True, + normalize_kernel=True, target_datatype=6, + target_nodata=numpy.nan) + + expected_output = pygeoprocessing.raster_to_numpy_array(target_path) + numpy.testing.assert_allclose(signal_array, expected_output) + self.assertEqual(1, numpy.count_nonzero(numpy.isnan(expected_output))) + def test_convolve_2d_error_on_worker_timeout(self): """PGP.geoprocessing: test convolve 2d error when worker times out.""" n_pixels = 10000 From 08e3a2e5e5bba5f7893c83a752c102709649c73a Mon Sep 17 00:00:00 2001 From: Claire Simpson Date: Wed, 11 Mar 2026 16:07:45 -0600 Subject: [PATCH 2/3] Add history note #473 --- HISTORY.rst | 3 +++ 1 file changed, 3 insertions(+) diff --git a/HISTORY.rst b/HISTORY.rst index 2fc61c9f..59160976 100644 --- a/HISTORY.rst +++ b/HISTORY.rst @@ -16,6 +16,9 @@ Unreleased Changes * Updating pyproject.toml to use the standard ``license-files`` key and replacing the license-related Trove classifier with the approved SPDX string. https://github.com/natcap/pygeoprocessing/issues/466 +* Fixed a bug where ``convolve_2d`` was not handling ``nan`` NoData + values correctly. + https://github.com/natcap/pygeoprocessing/issues/473 2.4.10 (2026-01-13) ------------------- From f961779c59b4fb68454a7797f0fb271530f44411 Mon Sep 17 00:00:00 2001 From: Claire Simpson <112011324+claire-simpson@users.noreply.github.com> Date: Wed, 11 Mar 2026 16:18:42 -0600 Subject: [PATCH 3/3] Fix order of history post-merge --- HISTORY.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/HISTORY.rst b/HISTORY.rst index bb96c2ff..e9a60f90 100644 --- a/HISTORY.rst +++ b/HISTORY.rst @@ -16,13 +16,13 @@ Unreleased Changes * Updating pyproject.toml to use the standard ``license-files`` key and replacing the license-related Trove classifier with the approved SPDX string. https://github.com/natcap/pygeoprocessing/issues/466 -* Fixed a bug where ``convolve_2d`` was not handling ``nan`` NoData - values correctly. - https://github.com/natcap/pygeoprocessing/issues/473 * Fixing a side effect in ``pygeoprocessing.align_and_resize_raster_stack`` where aligning a raster stack using the bounding box mode would result in the input bounding box being modified in-place. https://github.com/natcap/pygeoprocessing/issues/471 +* Fixed a bug where ``convolve_2d`` was not handling ``nan`` NoData + values correctly. + https://github.com/natcap/pygeoprocessing/issues/473 2.4.10 (2026-01-13) -------------------