Handle raster_map case with float32 array and float64 nodata values

HISTORY.rst

@@ -7,6 +7,9 @@ Release History
 ------------------
 * Updating for numpy 2.0 API changes.  Pygeoprocessing is now compatible with
   numpy 2.0 and later.  https://github.com/natcap/pygeoprocessing/issues/396
+* Handling a case in ``raster_map`` where an exception would be raised when a
+  float32 array was passed along with a float64 nodata value.
+  https://github.com/natcap/pygeoprocessing/issues/358
 
 2.4.4 (2024-05-21)
 ------------------

src/pygeoprocessing/geoprocessing.py

@@ -36,8 +36,9 @@
 from .geoprocessing_core import DEFAULT_CREATION_OPTIONS
 from .geoprocessing_core import DEFAULT_GTIFF_CREATION_TUPLE_OPTIONS
 from .geoprocessing_core import DEFAULT_OSR_AXIS_MAPPING_STRATEGY
+from .geoprocessing_core import gdal_use_exceptions
+from .geoprocessing_core import GDALUseExceptions
 from .geoprocessing_core import INT8_CREATION_OPTIONS
-from .geoprocessing_core import gdal_use_exceptions, GDALUseExceptions
 
 # This is used to efficiently pass data to the raster stats worker if available
 if sys.version_info >= (3, 8):
@@ -677,7 +678,7 @@ def raster_map(op, rasters, target_path, target_nodata=None, target_dtype=None,
         target_nodata = choose_nodata(target_dtype)
     else:
         if not numpy.can_cast(numpy.min_scalar_type(target_nodata),
-                              target_dtype):
+                              target_dtype, casting='same_kind'):
             raise ValueError(
                 f'Target nodata value {target_nodata} is incompatible with '
                 f'the target dtype {target_dtype}')
@@ -3825,7 +3826,9 @@ def get_gis_type(path):
         ``pygeoprocessing.RASTER_TYPE``, or ``pygeoprocessing.VECTOR_TYPE``.
 
     """
-    from pygeoprocessing import UNKNOWN_TYPE, RASTER_TYPE, VECTOR_TYPE
+    from pygeoprocessing import RASTER_TYPE
+    from pygeoprocessing import UNKNOWN_TYPE
+    from pygeoprocessing import VECTOR_TYPE
     gis_type = UNKNOWN_TYPE
     try:
         gis_raster = gdal.OpenEx(path, gdal.OF_RASTER)

tests/test_geoprocessing.py

@@ -32,10 +32,10 @@
 from pygeoprocessing.geoprocessing_core import DEFAULT_CREATION_OPTIONS
 from pygeoprocessing.geoprocessing_core import \
     DEFAULT_GTIFF_CREATION_TUPLE_OPTIONS
+from pygeoprocessing.geoprocessing_core import gdal_use_exceptions
 from pygeoprocessing.geoprocessing_core import INT8_CREATION_OPTIONS
 from pygeoprocessing.geoprocessing_core import \
     INT8_GTIFF_CREATION_TUPLE_OPTIONS
-from pygeoprocessing.geoprocessing_core import gdal_use_exceptions
 
 _DEFAULT_ORIGIN = (444720, 3751320)
 _DEFAULT_PIXEL_SIZE = (30, -30)
@@ -5371,6 +5371,21 @@ def test_raster_map_warn_on_multiband(self):
         self.assertEqual(log_messages[0].levelno, logging.WARNING)
         self.assertIn('has more than one band', log_messages[0].msg)
 
+    def test_raster_map_different_nodata_and_array_dtypes(self):
+        """PGP: raster_map can handle similar dtypes for nodata and arrays."""
+        band_1_array = numpy.array([[1, 1], [1, 1]], dtype=numpy.float32)
+        nodata = float(numpy.finfo(numpy.float32).min)  # this is a float64
+        source_path = os.path.join(self.workspace_dir, 'float32.tif')
+        _array_to_raster(band_1_array, nodata, source_path)
+
+        target_path = os.path.join(self.workspace_dir, 'target.tif')
+        pygeoprocessing.raster_map(lambda a: a, [source_path], target_path,
+                                   target_nodata=nodata)
+
+        expected_array = numpy.array(band_1_array, dtype=numpy.float32)
+        numpy.testing.assert_allclose(
+            expected_array, pygeoprocessing.raster_to_numpy_array(target_path))
+
     def test_choose_dtype(self):
         """PGP: choose_dtype picks smallest safe dtype for raster output."""
         uint8_raster = os.path.join(self.workspace_dir, 'uint8.tif')