Repository: natcap/invest (python)
Issue: CV does not handle `nan` in bathymetry
A raster with `nan` as `nodata` will yield `nan` values from `extract_bathymetry_along_ray` if any of the fetch ray overlaps nodata.

This looks like a case where we're using `math.isclose` where we should be using `numpy.isclose(x, y, equal_nan=True)`.
Pull request description: Bugfix CV: mask out nan nodata
We were not properly masking out `nan` values from the bathymetry raster.

I updated an `isclose` check and updated a test to confirm that `nan` values can be properly masked like any other nodata value. Other raster ops in this model were already doing the proper checks via `pygeoprocessing.array_equals_nodata`.

Fixes #1528 

## Checklist
- [x] 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 cf88c717b..3a0a55700 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -41,6 +41,11 @@ Unreleased Changes
 * Annual Water Yield
     * Added the results_suffix to a few intermediate files where it was
       missing. https://github.com/natcap/invest/issues/1517
+* Coastal Vulnerability
+    * Fixed a bug in handling ``nan`` as the nodata value of the bathymetry
+      raster. ``nan`` pixels will now be propertly ignored before calculating
+      mean depths along fetch rays.
+      https://github.com/natcap/invest/issues/1528
 * Urban Nature Access
     * Fixed a ``NameError`` that occurred when running the model using
       search radii defined per population group with an exponential search
diff --git a/src/natcap/invest/coastal_vulnerability.py b/src/natcap/invest/coastal_vulnerability.py
index ced65f9d7..8ebf03298 100644
--- a/src/natcap/invest/coastal_vulnerability.py
+++ b/src/natcap/invest/coastal_vulnerability.py
@@ -1754,7 +1754,8 @@ def extract_bathymetry_along_ray(
             raise ValueError(
                 f'got a {value} when trying to read bathymetry at {location}. '
                 'Does the bathymetry input fully cover the fetch ray area?')
-        if bathy_nodata is None or not math.isclose(value[0][0], bathy_nodata):
+        if bathy_nodata is None or not numpy.isclose(
+                value[0][0], bathy_nodata, equal_nan=True):
             bathy_values.append(value)
 
     # Gaps between shoreline and bathymetry input datasets could result in no
@@ -3160,7 +3161,7 @@ def _aggregate_raster_values_in_radius(
         max_distance=pixel_dist,
         normalize=False)
     radial_kernel_mask = pygeoprocessing.raster_to_numpy_array(
-    kernel_path).astype(bool)
+        kernel_path).astype(bool)
     shutil.rmtree(temp_dir, ignore_errors=True)
 
     result = {}


Test patch:
diff --git a/tests/test_coastal_vulnerability.py b/tests/test_coastal_vulnerability.py
index 9f3f7b51d..8d8a6a37d 100644
--- a/tests/test_coastal_vulnerability.py
+++ b/tests/test_coastal_vulnerability.py
@@ -224,27 +224,29 @@ class CoastalVulnerabilityTests(unittest.TestCase):
         projection_wkt = srs.ExportToWkt()
         geotransform = [0, 1.0, 0.0, 0, 0.0, -1.0]
         n = 5
-        nodata_val = 9999
+        band1_nodata = numpy.nan
         gtiff_driver = gdal.GetDriverByName('GTiff')
         new_raster = gtiff_driver.Create(
-            raster_path, n, n, 2, gdal.GDT_Int32, options=[
+            raster_path, n, n, 2, gdal.GDT_Float32, options=[
                 'TILED=YES', 'BIGTIFF=YES', 'COMPRESS=LZW',
                 'BLOCKXSIZE=16', 'BLOCKYSIZE=16'])
         new_raster.SetProjection(projection_wkt)
         new_raster.SetGeoTransform(geotransform)
-        array = numpy.array([-1]*n*n).reshape((n, n))
+        valid_value = -10.0
+        array = numpy.array([valid_value]*n*n).reshape((n, n))
 
         # nodata across the top row for Band 1
         new_band = new_raster.GetRasterBand(1)
-        array[:1] = nodata_val
+        array[:1] = band1_nodata
         new_band.WriteArray(array)
-        new_band.SetNoDataValue(nodata_val)
+        new_band.SetNoDataValue(band1_nodata)
 
         # all nodata for Band 2
+        band2_nodata = 999.99
         nodata_band = new_raster.GetRasterBand(2)
-        array[:] = nodata_val
+        array[:] = band2_nodata
         nodata_band.WriteArray(array)
-        nodata_band.SetNoDataValue(nodata_val)
+        nodata_band.SetNoDataValue(band2_nodata)
 
         new_raster.FlushCache()
         new_band = None
@@ -276,25 +278,25 @@ class CoastalVulnerabilityTests(unittest.TestCase):
         band = raster.GetRasterBand(1)  # nodata across top row
 
         values = coastal_vulnerability.extract_bathymetry_along_ray(
-            all_valid_ray, geotransform, nodata_val, band)
-        self.assertTrue(numpy.mean(values) == -1)
+            all_valid_ray, geotransform, band1_nodata, band)
+        self.assertEqual(numpy.mean(values), valid_value)
 
         values = coastal_vulnerability.extract_bathymetry_along_ray(
-            some_nodata_ray, geotransform, nodata_val, band)
-        self.assertTrue(numpy.mean(values) == -1)
+            some_nodata_ray, geotransform, band1_nodata, band)
+        self.assertEqual(numpy.mean(values), valid_value)
 
         values = coastal_vulnerability.extract_bathymetry_along_ray(
-            all_nodata_ray, geotransform, nodata_val, band)
-        self.assertTrue(numpy.mean(values) == -1)
+            all_nodata_ray, geotransform, band1_nodata, band)
+        self.assertEqual(numpy.mean(values), valid_value)
 
         with self.assertRaises(ValueError):
             values = coastal_vulnerability.extract_bathymetry_along_ray(
-                out_of_bounds_ray, geotransform, nodata_val, band)
+                out_of_bounds_ray, geotransform, band1_nodata, band)
 
         nodata_band = raster.GetRasterBand(2)  # all nodata band
         with self.assertRaises(ValueError):
             values = coastal_vulnerability.extract_bathymetry_along_ray(
-                all_valid_ray, geotransform, nodata_val, nodata_band)
+                all_valid_ray, geotransform, band2_nodata, nodata_band)
 
         raster = None
         band = None