Merge pull request #228 from djhoese/bugfix-scipy114-compat

Fix scipy 1.14 compatibility and trapz usage

pyspectral/blackbody.py

@@ -54,6 +54,9 @@ def blackbody_rad2temp(wavelength, radiance):
         The derived temperature in Kelvin.
 
     """
+    if getattr(wavelength, "dtype", None) != radiance.dtype:
+        # avoid a wavelength numpy scalar upcasting radiances (ex. 32-bit to 64-bit float)
+        wavelength = radiance.dtype.type(wavelength)
     with np.errstate(invalid='ignore'):
         return PLANCK_C1 / (wavelength * np.log(PLANCK_C2 / (radiance * wavelength**5) + 1.0))
 

pyspectral/near_infrared_reflectance.py

@@ -259,7 +259,7 @@ def reflectance_from_tbs(self, sun_zenith, tb_near_ir, tb_thermal, **kwargs):
         mu0 = np.cos(np.deg2rad(sunz))
 
         # mu0 = np.where(np.less(mu0, 0.1), 0.1, mu0)
-        self._solar_radiance = self.solar_flux * mu0 / np.pi
+        self._solar_radiance = (self.solar_flux * mu0 / np.pi).astype(tb_nir.dtype)
 
         # CO2 correction to the 3.9 radiance, only if tbs of a co2 band around
         # 13.4 micron is provided:

pyspectral/radiance_tb_conversion.py

@@ -28,12 +28,16 @@
 from numbers import Number
 
 import numpy as np
-from scipy import integrate
 
 from pyspectral.blackbody import C_SPEED, H_PLANCK, K_BOLTZMANN, blackbody, blackbody_wn
 from pyspectral.rsr_reader import RelativeSpectralResponse
 from pyspectral.utils import BANDNAMES, WAVE_LENGTH, WAVE_NUMBER, convert2wavenumber, get_bandname_from_wavelength
 
+try:
+    from scipy.integrate import trapezoid
+except ImportError:
+    from scipy.integrate import trapz as trapezoid
+
 LOG = logging.getLogger(__name__)
 
 BLACKBODY_FUNC = {WAVE_LENGTH: blackbody,
@@ -221,9 +225,9 @@ def tb2radiance(self, tb_, **kwargs):
 
         planck = self.blackbody_function(self.wavelength_or_wavenumber, tb_) * self.response
         if normalized:
-            radiance = integrate.trapz(planck, self.wavelength_or_wavenumber) / self.rsr_integral
+            radiance = trapezoid(planck, self.wavelength_or_wavenumber) / self.rsr_integral
         else:
-            radiance = integrate.trapz(planck, self.wavelength_or_wavenumber)
+            radiance = trapezoid(planck, self.wavelength_or_wavenumber)
 
         return {'radiance': radiance,
                 'unit': unit,

pyspectral/tests/test_rayleigh.py

@@ -136,6 +136,10 @@ def mocked_rsr():
         yield mymock
 
 
+def _create_dask_array(input_data, dtype):
+    return da.from_array(np.array(input_data, dtype=dtype))
+
+
 class TestRayleighDask:
     """Class for testing pyspectral.rayleigh - with dask-arrays as input."""
 
@@ -159,51 +163,34 @@ def test_get_reflectance_dask_redband_outside_clip(self, fake_lut_hdf5):
         assert isinstance(refl_corr, da.Array)
 
     @pytest.mark.parametrize("dtype", [np.float32, np.float64])
-    def test_get_reflectance_dask(self, fake_lut_hdf5, dtype):
+    @pytest.mark.parametrize("use_dask", [False, True])
+    @pytest.mark.parametrize(
+        ("input_data", "exp_result"),
+        [
+            (([67.0, 32.0], [45.0, 18.0], [150.0, 110.0], [14.0, 5.0]), TEST_RAYLEIGH_RESULT1),
+            (([60.0, 20.0], [49.0, 26.0], [140.0, 130.0], [12.0, 8.0]), TEST_RAYLEIGH_RESULT2),
+        ]
+    )
+    def test_get_reflectance(self, fake_lut_hdf5, dtype, use_dask, input_data, exp_result):
         """Test getting the reflectance correction with dask inputs."""
-        sun_zenith = da.array([67., 32.], dtype=dtype)
-        sat_zenith = da.array([45., 18.], dtype=dtype)
-        azidiff = da.array([150., 110.], dtype=dtype)
-        redband_refl = da.array([14., 5.], dtype=dtype)
+        array_func = np.array if not use_dask else _create_dask_array
+        sun_zenith = array_func(input_data[0], dtype=dtype)
+        sat_zenith = array_func(input_data[1], dtype=dtype)
+        azidiff = array_func(input_data[2], dtype=dtype)
+        redband_refl = array_func(input_data[3], dtype=dtype)
         rayl = _create_rayleigh()
         with mocked_rsr():
             refl_corr = rayl.get_reflectance(sun_zenith, sat_zenith, azidiff, 'ch3', redband_refl)
-        np.testing.assert_allclose(refl_corr, TEST_RAYLEIGH_RESULT1.astype(dtype), atol=4.0e-06)
-        assert isinstance(refl_corr, da.Array)
-        assert refl_corr.dtype == dtype  # check that the dask array's dtype is equal
-        assert refl_corr.compute().dtype == dtype  # check that the final numpy array's dtype is equal
-
-        sun_zenith = da.array([60., 20.], dtype=dtype)
-        sat_zenith = da.array([49., 26.], dtype=dtype)
-        azidiff = da.array([140., 130.], dtype=dtype)
-        redband_refl = da.array([12., 8.], dtype=dtype)
-        with mocked_rsr():
-            refl_corr = rayl.get_reflectance(sun_zenith, sat_zenith, azidiff, 'ch3', redband_refl)
-        np.testing.assert_allclose(refl_corr, TEST_RAYLEIGH_RESULT2.astype(dtype), atol=4.0e-06)
-        assert isinstance(refl_corr, da.Array)
-        assert refl_corr.dtype == dtype  # check that the dask array's dtype is equal
-        assert refl_corr.compute().dtype == dtype  # check that the final numpy array's dtype is equal
 
-    def test_get_reflectance_numpy_dask(self, fake_lut_hdf5):
-        """Test getting the reflectance correction with dask inputs."""
-        sun_zenith = np.array([67., 32.])
-        sat_zenith = np.array([45., 18.])
-        azidiff = np.array([150., 110.])
-        redband_refl = np.array([14., 5.])
-        rayl = _create_rayleigh()
-        with mocked_rsr():
-            refl_corr = rayl.get_reflectance(sun_zenith, sat_zenith, azidiff, 'ch3', redband_refl)
-        np.testing.assert_allclose(refl_corr, TEST_RAYLEIGH_RESULT1)
-        assert isinstance(refl_corr, np.ndarray)
+        if use_dask:
+            assert isinstance(refl_corr, da.Array)
+            refl_corr_np = refl_corr.compute()
+            assert refl_corr_np.dtype == refl_corr.dtype  # check that the final numpy array's dtype is equal
+            refl_corr = refl_corr_np
 
-        sun_zenith = np.array([60., 20.])
-        sat_zenith = np.array([49., 26.])
-        azidiff = np.array([140., 130.])
-        redband_refl = np.array([12., 8.])
-        with mocked_rsr():
-            refl_corr = rayl.get_reflectance(sun_zenith, sat_zenith, azidiff, 'ch3', redband_refl)
-        np.testing.assert_allclose(refl_corr, TEST_RAYLEIGH_RESULT2)
         assert isinstance(refl_corr, np.ndarray)
+        np.testing.assert_allclose(refl_corr, exp_result.astype(dtype), atol=4.0e-06)
+        assert refl_corr.dtype == dtype  # check that the dask array's dtype is equal
 
     def test_get_reflectance_wvl_outside_range(self, fake_lut_hdf5):
         """Test getting the reflectance correction with wavelength outside correction range."""
@@ -347,14 +334,19 @@ def test_get_reflectance_redband_outside_clip(self, fake_lut_hdf5):
              TEST_RAYLEIGH_RESULT5),
         ]
     )
-    def test_get_reflectance(self, fake_lut_hdf5, sun_zenith, sat_zenith, azidiff, redband_refl, exp_result):
+    @pytest.mark.parametrize("dtype", [np.float32, np.float64])
+    def test_get_reflectance(self, fake_lut_hdf5, sun_zenith, sat_zenith, azidiff, redband_refl, exp_result, dtype):
         """Test getting the reflectance correction."""
         rayl = _create_rayleigh()
         with mocked_rsr():
             refl_corr = rayl.get_reflectance(
-                sun_zenith, sat_zenith, azidiff, 'ch3', redband_refl)
+                sun_zenith.astype(dtype),
+                sat_zenith.astype(dtype),
+                azidiff.astype(dtype),
+                'ch3',
+                redband_refl.astype(dtype))
         assert isinstance(refl_corr, np.ndarray)
-        np.testing.assert_allclose(refl_corr, exp_result)
+        np.testing.assert_allclose(refl_corr, exp_result.astype(dtype), atol=4.0e-06)
 
     @patch('pyspectral.rayleigh.da', None)
     def test_get_reflectance_no_rsr(self, fake_lut_hdf5):

pyspectral/tests/test_utils.py

@@ -225,11 +225,11 @@ def test_get_wave_range(self):
 @pytest.mark.parametrize(
     ("input_value", "exp_except"),
     [
-        (np.string_("hey"), False),
-        (np.array([np.string_("hey")]), False),
-        (np.array(np.string_("hey")), False),
+        (np.bytes_("hey"), False),
+        (np.array([np.bytes_("hey")]), False),
+        (np.array(np.bytes_("hey")), False),
         ("hey", False),
-        (np.array([np.string_("hey"), np.string_("hey")]), True),
+        (np.array([np.bytes_("hey"), np.bytes_("hey")]), True),
         (5, True),
     ],
 )
@@ -247,8 +247,8 @@ def test_np2str(input_value, exp_except):
     [
         (b"Hello", "Hello"),
         ("Hello", "Hello"),
-        (np.string_("Hello"), "Hello"),
-        (np.array(np.string_("Hello")), np.array(np.string_("Hello"))),
+        (np.bytes_("Hello"), "Hello"),
+        (np.array(np.bytes_("Hello")), np.array(np.bytes_("Hello"))),
     ]
 )
 def test_bytes2string(input_value, exp_result):

pyspectral/utils.py

@@ -506,20 +506,23 @@ def convert2str(value):
 
 
 def np2str(value):
-    """Convert an `numpy.string_` to str.
+    """Convert an ``numpy.bytes_`` to str.
+
+    Note: ``numpy.string_`` was deprecated in numpy 2.0 in favor of
+    ``numpy.bytes_``.
 
     Args:
         value (ndarray): scalar or 1-element numpy array to convert
     Raises:
         ValueError: if value is array larger than 1-element or it is not of
-                    type `numpy.string_` or it is not a numpy array
+                    type `numpy.bytes_` or it is not a numpy array
 
     """
     if isinstance(value, str):
         return value
 
     if hasattr(value, 'dtype') and \
-            issubclass(value.dtype.type, (np.str_, np.string_, np.object_)) \
+            issubclass(value.dtype.type, (np.str_, np.bytes_, np.object_)) \
             and value.size == 1:
         value = value.item()
         # python 3 - was scalar numpy array of bytes