diff --git a/HISTORY.rst b/HISTORY.rst
index d9f4fba2..e9a60f90 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -20,6 +20,9 @@ Unreleased Changes
   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)
 -------------------
diff --git a/src/pygeoprocessing/geoprocessing.py b/src/pygeoprocessing/geoprocessing.py
index 70847b3a..fbf5c0e5 100644
--- a/src/pygeoprocessing/geoprocessing.py
+++ b/src/pygeoprocessing/geoprocessing.py
@@ -3338,7 +3338,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
@@ -3429,7 +3429,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,
@@ -3477,7 +3477,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)
@@ -4074,7 +4074,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:
@@ -4093,7 +4093,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
@@ -4125,7 +4125,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 a2ef4d38..37678e36 100644
--- a/tests/test_geoprocessing.py
+++ b/tests/test_geoprocessing.py
@@ -3730,6 +3730,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