diff --git a/pyspectral/blackbody.py b/pyspectral/blackbody.py index a684684..901a835 100644 --- a/pyspectral/blackbody.py +++ b/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)) diff --git a/pyspectral/near_infrared_reflectance.py b/pyspectral/near_infrared_reflectance.py index 1c86ce2..807b2a9 100644 --- a/pyspectral/near_infrared_reflectance.py +++ b/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: diff --git a/pyspectral/radiance_tb_conversion.py b/pyspectral/radiance_tb_conversion.py index ca6a2f5..8521b8e 100644 --- a/pyspectral/radiance_tb_conversion.py +++ b/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, diff --git a/pyspectral/tests/test_rayleigh.py b/pyspectral/tests/test_rayleigh.py index f11dbe5..ec8ca80 100644 --- a/pyspectral/tests/test_rayleigh.py +++ b/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): diff --git a/pyspectral/tests/test_utils.py b/pyspectral/tests/test_utils.py index 3969bc0..c54bd8d 100644 --- a/pyspectral/tests/test_utils.py +++ b/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): diff --git a/pyspectral/utils.py b/pyspectral/utils.py index 425db6d..2c1b002 100644 --- a/pyspectral/utils.py +++ b/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