From adde5c16722bfe2384d6eb182e0509082b3e2faf Mon Sep 17 00:00:00 2001 From: David Hoese Date: Thu, 27 Jun 2024 09:32:49 -0500 Subject: [PATCH 1/8] Fix scipy 1.14 compatibility and trapz usage --- pyspectral/radiance_tb_conversion.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) 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, From 41265c79eceadc0dccf18ba0f4ea9a48442b1d10 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Thu, 27 Jun 2024 09:48:43 -0500 Subject: [PATCH 2/8] Replace legacy np.string_ usage with np.bytes_ --- pyspectral/tests/test_utils.py | 12 ++++++------ pyspectral/utils.py | 9 ++++++--- 2 files changed, 12 insertions(+), 9 deletions(-) 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 From 5c4508f155188cdb85648c329b7ce7e628487782 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Thu, 27 Jun 2024 10:18:51 -0500 Subject: [PATCH 3/8] Fix blackbody upcasting input radiances with float64 wavelength Numpy 2 used np.float64 wavelength as reason to upcast. --- pyspectral/blackbody.py | 3 +++ 1 file changed, 3 insertions(+) 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)) From 5b4e7556d184bdb3be36a1063b58491b57d5c6a9 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Thu, 27 Jun 2024 10:19:07 -0500 Subject: [PATCH 4/8] Fix dtype usage in numpy rayleigh test --- pyspectral/tests/test_rayleigh.py | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/pyspectral/tests/test_rayleigh.py b/pyspectral/tests/test_rayleigh.py index f11dbe5..8da5199 100644 --- a/pyspectral/tests/test_rayleigh.py +++ b/pyspectral/tests/test_rayleigh.py @@ -184,25 +184,26 @@ def test_get_reflectance_dask(self, fake_lut_hdf5, dtype): 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): + @pytest.mark.parametrize("dtype", [np.float32, np.float64]) + def test_get_reflectance_numpy(self, fake_lut_hdf5, dtype): """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.]) + sun_zenith = np.array([67., 32.], dtype=dtype) + sat_zenith = np.array([45., 18.], dtype=dtype) + azidiff = np.array([150., 110.], dtype=dtype) + redband_refl = np.array([14., 5.], 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) + np.testing.assert_allclose(refl_corr, TEST_RAYLEIGH_RESULT1.astype(dtype), atol=4.0e-06) assert isinstance(refl_corr, np.ndarray) - sun_zenith = np.array([60., 20.]) - sat_zenith = np.array([49., 26.]) - azidiff = np.array([140., 130.]) - redband_refl = np.array([12., 8.]) + sun_zenith = np.array([60., 20.], dtype=dtype) + sat_zenith = np.array([49., 26.], dtype=dtype) + azidiff = np.array([140., 130.], dtype=dtype) + redband_refl = np.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) + np.testing.assert_allclose(refl_corr, TEST_RAYLEIGH_RESULT2.astype(dtype), atol=4.0e-06) assert isinstance(refl_corr, np.ndarray) def test_get_reflectance_wvl_outside_range(self, fake_lut_hdf5): @@ -354,6 +355,7 @@ def test_get_reflectance(self, fake_lut_hdf5, sun_zenith, sat_zenith, azidiff, r refl_corr = rayl.get_reflectance( sun_zenith, sat_zenith, azidiff, 'ch3', redband_refl) assert isinstance(refl_corr, np.ndarray) + print(refl_corr.dtype) np.testing.assert_allclose(refl_corr, exp_result) @patch('pyspectral.rayleigh.da', None) From 8d32f34a63f255199a84902debd071ebb99a8a07 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Thu, 27 Jun 2024 10:26:27 -0500 Subject: [PATCH 5/8] Add better parametrize in rayleigh test --- pyspectral/tests/test_rayleigh.py | 62 +++++++++++++++---------------- 1 file changed, 29 insertions(+), 33 deletions(-) diff --git a/pyspectral/tests/test_rayleigh.py b/pyspectral/tests/test_rayleigh.py index 8da5199..2b96f0f 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,52 +163,44 @@ 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]) + def test_get_reflectance(self, fake_lut_hdf5, dtype, use_dask): """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([67., 32.], dtype=dtype) + sat_zenith = array_func([45., 18.], dtype=dtype) + azidiff = array_func([150., 110.], dtype=dtype) + redband_refl = array_func([14., 5.], dtype=dtype) rayl = _create_rayleigh() with mocked_rsr(): refl_corr = rayl.get_reflectance(sun_zenith, sat_zenith, azidiff, 'ch3', redband_refl) + + 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 + + assert isinstance(refl_corr, np.ndarray) 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) + sun_zenith = array_func([60., 20.], dtype=dtype) + sat_zenith = array_func([49., 26.], dtype=dtype) + azidiff = array_func([140., 130.], dtype=dtype) + redband_refl = array_func([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 - @pytest.mark.parametrize("dtype", [np.float32, np.float64]) - def test_get_reflectance_numpy(self, fake_lut_hdf5, dtype): - """Test getting the reflectance correction with dask inputs.""" - sun_zenith = np.array([67., 32.], dtype=dtype) - sat_zenith = np.array([45., 18.], dtype=dtype) - azidiff = np.array([150., 110.], dtype=dtype) - redband_refl = np.array([14., 5.], 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, 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.], dtype=dtype) - sat_zenith = np.array([49., 26.], dtype=dtype) - azidiff = np.array([140., 130.], dtype=dtype) - redband_refl = np.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, np.ndarray) + 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.""" From fedaed8411f048c8f6b3a38e083ee3b27a6bd3d7 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Thu, 27 Jun 2024 10:29:12 -0500 Subject: [PATCH 6/8] More parametrization in rayleigh tests --- pyspectral/tests/test_rayleigh.py | 36 +++++++++++-------------------- 1 file changed, 13 insertions(+), 23 deletions(-) diff --git a/pyspectral/tests/test_rayleigh.py b/pyspectral/tests/test_rayleigh.py index 2b96f0f..c6b97ea 100644 --- a/pyspectral/tests/test_rayleigh.py +++ b/pyspectral/tests/test_rayleigh.py @@ -164,13 +164,20 @@ def test_get_reflectance_dask_redband_outside_clip(self, fake_lut_hdf5): @pytest.mark.parametrize("dtype", [np.float32, np.float64]) @pytest.mark.parametrize("use_dask", [False, True]) - def test_get_reflectance(self, fake_lut_hdf5, dtype, use_dask): + @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.""" array_func = np.array if not use_dask else _create_dask_array - sun_zenith = array_func([67., 32.], dtype=dtype) - sat_zenith = array_func([45., 18.], dtype=dtype) - azidiff = array_func([150., 110.], dtype=dtype) - redband_refl = array_func([14., 5.], dtype=dtype) + 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) @@ -182,24 +189,7 @@ def test_get_reflectance(self, fake_lut_hdf5, dtype, use_dask): refl_corr = refl_corr_np assert isinstance(refl_corr, np.ndarray) - np.testing.assert_allclose(refl_corr, TEST_RAYLEIGH_RESULT1.astype(dtype), atol=4.0e-06) - assert refl_corr.dtype == dtype # check that the dask array's dtype is equal - - sun_zenith = array_func([60., 20.], dtype=dtype) - sat_zenith = array_func([49., 26.], dtype=dtype) - azidiff = array_func([140., 130.], dtype=dtype) - redband_refl = array_func([12., 8.], dtype=dtype) - with mocked_rsr(): - refl_corr = rayl.get_reflectance(sun_zenith, sat_zenith, azidiff, 'ch3', redband_refl) - - 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 - - np.testing.assert_allclose(refl_corr, TEST_RAYLEIGH_RESULT2.astype(dtype), atol=4.0e-06) - 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): From 3152387c4e54eabd51680cdd474c78433845469e Mon Sep 17 00:00:00 2001 From: David Hoese Date: Thu, 27 Jun 2024 10:43:35 -0500 Subject: [PATCH 7/8] Fix more dtype handling in rayleigh tests --- pyspectral/tests/test_rayleigh.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/pyspectral/tests/test_rayleigh.py b/pyspectral/tests/test_rayleigh.py index c6b97ea..ec8ca80 100644 --- a/pyspectral/tests/test_rayleigh.py +++ b/pyspectral/tests/test_rayleigh.py @@ -334,15 +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) - print(refl_corr.dtype) - 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): From 239c410622ea0472c5f19f2745bde25b72c5dbf4 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Thu, 27 Jun 2024 10:44:03 -0500 Subject: [PATCH 8/8] Fix solar_flux 64-bit float causing upcasting in NIR reflectance calculations --- pyspectral/near_infrared_reflectance.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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: