diff --git a/pyspectral/blackbody.py b/pyspectral/blackbody.py index ed9faf1..a684684 100644 --- a/pyspectral/blackbody.py +++ b/pyspectral/blackbody.py @@ -99,24 +99,12 @@ def planck(wave, temperature, wavelength=True): LOG.debug("Using {0} when calculating the Blackbody radiance".format( units[(wavelength is True) - 1])) - if np.isscalar(temperature): - temperature = np.array([temperature, ], dtype='float64') - elif isinstance(temperature, (list, tuple)): - temperature = np.array(temperature, dtype='float64') - - shape = temperature.shape - if np.isscalar(wave): - wln = np.array([wave, ], dtype='float64') - else: - wln = np.array(wave, dtype='float64') - if wavelength: - const = 2 * H_PLANCK * C_SPEED ** 2 - nom = const / wln ** 5 - arg1 = H_PLANCK * C_SPEED / (K_BOLTZMANN * wln) + nom = PLANCK_C2 / wave ** 5 + arg1 = PLANCK_C1 / wave else: - nom = 2 * H_PLANCK * (C_SPEED ** 2) * (wln ** 3) - arg1 = H_PLANCK * C_SPEED * wln / K_BOLTZMANN + nom = PLANCK_C2 * wave ** 3 + arg1 = PLANCK_C1 * wave with np.errstate(divide='ignore', invalid='ignore'): # use dask functions when needed @@ -132,7 +120,7 @@ def planck(wave, temperature, wavelength=True): str(np.nanmax(arg2)), str(np.nanmin(arg2))) try: - exp_arg = np.multiply(arg1.astype('float64'), arg2.astype('float64')) + exp_arg = np.multiply(arg1, arg2) except MemoryError: LOG.warning(("Dimensions used in numpy.multiply probably reached " "limit!\n" @@ -151,21 +139,7 @@ def planck(wave, temperature, wavelength=True): with np.errstate(over='ignore'): denom = np.exp(exp_arg) - 1 - rad = nom / denom - radshape = rad.shape - if wln.shape[0] == 1: - if temperature.shape[0] == 1: - return rad[0, 0] - else: - return rad[:, 0].reshape(shape) - else: - if temperature.shape[0] == 1: - return rad[0, :] - else: - if len(shape) == 1: - return rad.reshape((shape[0], radshape[1])) - else: - return rad.reshape((shape[0], shape[1], radshape[1])) + return nom / denom def blackbody_wn(wavenumber, temp): diff --git a/pyspectral/near_infrared_reflectance.py b/pyspectral/near_infrared_reflectance.py index 080e22b..920c0ff 100644 --- a/pyspectral/near_infrared_reflectance.py +++ b/pyspectral/near_infrared_reflectance.py @@ -243,11 +243,11 @@ def reflectance_from_tbs(self, sun_zenith, tb_near_ir, tb_thermal, **kwargs): # Assume rsr is in microns!!! # FIXME! self._rad3x_t11 = self.tb2radiance(tb_therm, lut=lut)['radiance'] - thermal_emiss_one = self._rad3x_t11 * self.rsr_integral - + rsr_integral = tb_therm.dtype.type(self.rsr_integral) + thermal_emiss_one = self._rad3x_t11 * rsr_integral l_nir = self.tb2radiance(tb_nir, lut=lut)['radiance'] self._rad3x = l_nir.copy() - l_nir *= self.rsr_integral + l_nir *= rsr_integral if thermal_emiss_one.ravel().shape[0] < 10: LOG.info('thermal_emiss_one = %s', str(thermal_emiss_one)) @@ -268,7 +268,8 @@ def reflectance_from_tbs(self, sun_zenith, tb_near_ir, tb_thermal, **kwargs): self.derive_rad39_corr(tb_therm, tbco2) LOG.info("CO2 correction applied...") else: - self._rad3x_correction = 1.0 + self._rad3x_correction = np.float64(1.0) + self._rad3x_correction = self._rad3x_correction.astype(tb_nir.dtype) corrected_thermal_emiss_one = thermal_emiss_one * self._rad3x_correction nomin = l_nir - corrected_thermal_emiss_one diff --git a/pyspectral/radiance_tb_conversion.py b/pyspectral/radiance_tb_conversion.py index d246065..ca6a2f5 100644 --- a/pyspectral/radiance_tb_conversion.py +++ b/pyspectral/radiance_tb_conversion.py @@ -201,15 +201,16 @@ def tb2radiance(self, tb_, **kwargs): scale = 1.0 if lut: + lut_radiance = lut['radiance'].astype(tb_.dtype) ntb = (tb_ * self.tb_scale).astype('int16') start = int(lut['tb'][0] * self.tb_scale) retv = {} - bounds = 0, lut['radiance'].shape[0] - 1 + bounds = 0, lut_radiance.shape[0] - 1 index = (ntb - start).clip(bounds[0], bounds[1]) try: - retv['radiance'] = index.map_blocks(self._getitem, lut['radiance'], dtype=lut['radiance'].dtype) + retv['radiance'] = index.map_blocks(self._getitem, lut_radiance, dtype=lut_radiance.dtype) except AttributeError: - retv['radiance'] = lut['radiance'][index] + retv['radiance'] = lut_radiance[index] try: retv['radiance'] = retv['radiance'].item() except (ValueError, AttributeError): diff --git a/pyspectral/tests/test_blackbody.py b/pyspectral/tests/test_blackbody.py index a40e2d5..a3c1cc5 100644 --- a/pyspectral/tests/test_blackbody.py +++ b/pyspectral/tests/test_blackbody.py @@ -19,12 +19,13 @@ """Unit testing the Blackbody/Plack radiation derivation.""" -import unittest - +import dask +import dask.array as da import numpy as np +import pytest from pyspectral.blackbody import blackbody, blackbody_rad2temp, blackbody_wn, blackbody_wn_rad2temp -from pyspectral.tests.unittest_helpers import assertNumpyArraysEqual +from pyspectral.tests.unittest_helpers import ComputeCountingScheduler RAD_11MICRON_300KELVIN = 9573176.935507433 RAD_11MICRON_301KELVIN = 9714686.576498277 @@ -36,92 +37,65 @@ __unittest = True -class CustomScheduler(object): - """Custom dask scheduler that raises an exception if dask is computed too many times.""" - - def __init__(self, max_computes=1): - """Set starting and maximum compute counts.""" - self.max_computes = max_computes - self.total_computes = 0 - - def __call__(self, dsk, keys, **kwargs): - """Compute dask task and keep track of number of times we do so.""" - import dask - self.total_computes += 1 - if self.total_computes > self.max_computes: - raise RuntimeError("Too many dask computations were scheduled: {}".format(self.total_computes)) - return dask.get(dsk, keys, **kwargs) - - -class TestBlackbody(unittest.TestCase): +class TestBlackbody: """Unit testing the blackbody function.""" def test_blackbody(self): """Calculate the blackbody radiation from wavelengths and temperatures.""" wavel = 11. * 1E-6 - black = blackbody((wavel, ), [300., 301]) - self.assertEqual(black.shape[0], 2) - self.assertAlmostEqual(black[0], RAD_11MICRON_300KELVIN) - self.assertAlmostEqual(black[1], RAD_11MICRON_301KELVIN) - - temp1 = blackbody_rad2temp(wavel, black[0]) - self.assertAlmostEqual(temp1, 300.0, 4) - temp2 = blackbody_rad2temp(wavel, black[1]) - self.assertAlmostEqual(temp2, 301.0, 4) + black = blackbody(wavel, 300.) + assert black == pytest.approx(RAD_11MICRON_300KELVIN) + temp1 = blackbody_rad2temp(wavel, black) + assert temp1 == pytest.approx(300.0, 4) - black = blackbody(13. * 1E-6, 200.) - self.assertTrue(np.isscalar(black)) + black = blackbody(wavel, 301.) + assert black == pytest.approx(RAD_11MICRON_301KELVIN) + temp2 = blackbody_rad2temp(wavel, black) + assert temp2 == pytest.approx(301.0, 4) tb_therm = np.array([[300., 301], [299, 298], [279, 286]]) - black = blackbody((10. * 1E-6, 11.e-6), tb_therm) - self.assertIsInstance(black, np.ndarray) + black = blackbody(10. * 1E-6, tb_therm) + assert isinstance(black, np.ndarray) - tb_therm = np.array([[300., 301], [0., 298], [279, 286]]) - black = blackbody((10. * 1E-6, 11.e-6), tb_therm) - self.assertIsInstance(black, np.ndarray) - - def test_blackbody_dask(self): + def test_blackbody_dask_wave_tuple(self): """Calculate the blackbody radiation from wavelengths and temperatures with dask arrays.""" - import dask - import dask.array as da - tb_therm = da.from_array([[300., 301], [299, 298], [279, 286]], chunks=2) - with dask.config.set(scheduler=CustomScheduler(0)): - black = blackbody((10. * 1E-6, 11.e-6), tb_therm) - self.assertIsInstance(black, da.Array) - - tb_therm = da.from_array([[300., 301], [0., 298], [279, 286]], chunks=2) - with dask.config.set(scheduler=CustomScheduler(0)): - black = blackbody((10. * 1E-6, 11.e-6), tb_therm) - self.assertIsInstance(black, da.Array) + tb_therm = da.array([[300., 301], [299, 298], [279, 286]]) + with dask.config.set(scheduler=ComputeCountingScheduler(0)): + black = blackbody(10. * 1E-6, tb_therm) + assert isinstance(black, da.Array) + assert black.dtype == np.float64 + + @pytest.mark.parametrize("dtype", (np.float32, np.float64, float)) + def test_blackbody_dask_wave_array(self, dtype): + """Test blackbody calculations with dask arrays as inputs.""" + tb_therm = da.array([[300., 301], [0., 298], [279, 286]], dtype=dtype) + with dask.config.set(scheduler=ComputeCountingScheduler(0)): + black = blackbody(da.array([10. * 1E-6, 11.e-6], dtype=dtype), tb_therm) + assert isinstance(black, da.Array) + assert black.dtype == dtype def test_blackbody_wn(self): """Calculate the blackbody radiation from wavenumbers and temperatures.""" wavenumber = 90909.1 # 11 micron band - black = blackbody_wn((wavenumber, ), [300., 301]) - self.assertEqual(black.shape[0], 2) - self.assertAlmostEqual(black[0], WN_RAD_11MICRON_300KELVIN) - self.assertAlmostEqual(black[1], WN_RAD_11MICRON_301KELVIN) + black = blackbody_wn(wavenumber, 300.) + assert black == pytest.approx(WN_RAD_11MICRON_300KELVIN) + temp1 = blackbody_wn_rad2temp(wavenumber, black) + assert temp1 == pytest.approx(300.0, 4) - temp1 = blackbody_wn_rad2temp(wavenumber, black[0]) - self.assertAlmostEqual(temp1, 300.0, 4) - temp2 = blackbody_wn_rad2temp(wavenumber, black[1]) - self.assertAlmostEqual(temp2, 301.0, 4) + black = blackbody_wn(wavenumber, 301.) + assert black == pytest.approx(WN_RAD_11MICRON_301KELVIN) + temp2 = blackbody_wn_rad2temp(wavenumber, black) + assert temp2 == pytest.approx(301.0, 4) t__ = blackbody_wn_rad2temp(wavenumber, np.array([0.001, 0.0009])) expected = [290.3276916, 283.76115441] - self.assertAlmostEqual(t__[0], expected[0]) - self.assertAlmostEqual(t__[1], expected[1]) + np.testing.assert_allclose(t__, expected) radiances = np.array([0.001, 0.0009, 0.0012, 0.0018]).reshape(2, 2) t__ = blackbody_wn_rad2temp(wavenumber, radiances) expected = np.array([290.3276916, 283.76115441, 302.4181330, 333.1414164]).reshape(2, 2) - self.assertAlmostEqual(t__[1, 1], expected[1, 1], 5) - self.assertAlmostEqual(t__[0, 0], expected[0, 0], 5) - self.assertAlmostEqual(t__[0, 1], expected[0, 1], 5) - self.assertAlmostEqual(t__[1, 0], expected[1, 0], 5) - - assertNumpyArraysEqual(t__, expected) + np.testing.assert_allclose(t__, expected) def test_blackbody_wn_dask(self): """Test that blackbody rad2temp preserves dask arrays.""" @@ -129,15 +103,10 @@ def test_blackbody_wn_dask(self): import dask.array as da wavenumber = 90909.1 # 11 micron band radiances = da.from_array([0.001, 0.0009, 0.0012, 0.0018], chunks=2).reshape(2, 2) - with dask.config.set(scheduler=CustomScheduler(0)): + with dask.config.set(scheduler=ComputeCountingScheduler(0)): t__ = blackbody_wn_rad2temp(wavenumber, radiances) - self.assertIsInstance(t__, da.Array) + assert isinstance(t__, da.Array) t__ = t__.compute() expected = np.array([290.3276916, 283.76115441, 302.4181330, 333.1414164]).reshape(2, 2) - self.assertAlmostEqual(t__[1, 1], expected[1, 1], 5) - self.assertAlmostEqual(t__[0, 0], expected[0, 0], 5) - self.assertAlmostEqual(t__[0, 1], expected[0, 1], 5) - self.assertAlmostEqual(t__[1, 0], expected[1, 0], 5) - - assertNumpyArraysEqual(t__, expected) + np.testing.assert_allclose(t__, expected) diff --git a/pyspectral/tests/test_rad_tb_conversions.py b/pyspectral/tests/test_rad_tb_conversions.py index 0ed6aa0..3c0688c 100644 --- a/pyspectral/tests/test_rad_tb_conversions.py +++ b/pyspectral/tests/test_rad_tb_conversions.py @@ -29,7 +29,7 @@ from pyspectral.radiance_tb_conversion import RadTbConverter, SeviriRadTbConverter, radiance2tb from pyspectral.utils import get_central_wave -TEST_TBS = np.array([200., 270., 300., 302., 350.], dtype='float32') +TEST_TBS = np.array([200., 270., 300., 302., 350.]) TRUE_RADS = np.array([856.937353205, 117420.385297, 479464.582505, 521412.9511, 2928735.18944], @@ -38,8 +38,7 @@ 2.559173e-06, 9.797091e-06, 1.061431e-05, - 5.531423e-05], - dtype='float64') + 5.531423e-05]) TEST_RSR = {'20': {}} TEST_RSR['20']['det-1'] = {} @@ -224,10 +223,10 @@ def __init__(self): self.rsr = TEST_RSR -class TestSeviriConversions(unittest.TestCase): +class TestSeviriConversions: """Testing the conversions between radiances and brightness temperatures.""" - def setUp(self): + def setup_method(self): """Set up.""" with patch('pyspectral.radiance_tb_conversion.RelativeSpectralResponse') as mymock: instance = mymock.return_value @@ -240,32 +239,38 @@ def setUp(self): self.sev2 = SeviriRadTbConverter('Meteosat-9', 'IR3.9') - def test_rad2tb(self): + @pytest.mark.parametrize("dtype", (np.float32, np.float64, float)) + def test_rad2tb(self, dtype): """Unit testing the radiance to brightness temperature conversion.""" - res = self.sev1.tb2radiance(TEST_TBS, lut=False) - self.assertTrue(np.allclose(TRUE_RADS_SEVIRI, res['radiance'])) + res = self.sev1.tb2radiance(TEST_TBS.astype(dtype), lut=False) + assert res['radiance'].dtype == dtype + np.testing.assert_allclose(TRUE_RADS_SEVIRI.astype(dtype), res['radiance'], atol=1e-8) - def test_conversion_simple(self): + @pytest.mark.parametrize("dtype", (np.float32, np.float64, float)) + def test_conversion_simple(self, dtype): """Test the conversion based on the non-linear approximation (SEVIRI). Test the tb2radiance function to convert radiances to Tb's using tabulated coefficients based on a non-linear approximation """ - retv = self.sev2.tb2radiance(TEST_TBS) + retv = self.sev2.tb2radiance(TEST_TBS.astype(dtype)) rads = retv['radiance'] # Units space = wavenumber (cm-1): - tbs = self.sev2.radiance2tb(rads) - self.assertTrue(np.allclose(TEST_TBS, tbs)) + tbs = self.sev2.radiance2tb(rads.astype(dtype)) + assert tbs.dtype == dtype + np.testing.assert_allclose(TEST_TBS.astype(dtype), tbs, rtol=1e-6) np.random.seed() tbs1 = 200.0 + np.random.random(50) * 150.0 - retv = self.sev2.tb2radiance(tbs1) + retv = self.sev2.tb2radiance(tbs1.astype(dtype)) rads = retv['radiance'] tbs = self.sev2.radiance2tb(rads) - self.assertTrue(np.allclose(tbs1, tbs)) + assert tbs.dtype == dtype + np.testing.assert_allclose(tbs1, tbs, rtol=1e-6) - def test_conversions_methods(self): + @pytest.mark.parametrize("dtype", (np.float32, np.float64, float)) + def test_conversions_methods(self, dtype): """Test the conversion methods. Using the two diferent conversion methods to verify that they give @@ -274,12 +279,14 @@ def test_conversions_methods(self): """ # Units space = wavenumber (cm-1): - retv2 = self.sev2.tb2radiance(TEST_TBS) - retv1 = self.sev1.tb2radiance(TEST_TBS) + retv2 = self.sev2.tb2radiance(TEST_TBS.astype(dtype)) + retv1 = self.sev1.tb2radiance(TEST_TBS.astype(dtype)) rads1 = retv1['radiance'] rads2 = retv2['radiance'] - self.assertTrue(np.allclose(rads1, rads2)) + assert rads1.dtype == dtype + assert rads2.dtype == dtype + np.testing.assert_allclose(rads1, rads2, atol=1e-8) class TestRadTbConversions(unittest.TestCase): diff --git a/pyspectral/tests/test_reflectance.py b/pyspectral/tests/test_reflectance.py index 06d33ea..3d0af56 100644 --- a/pyspectral/tests/test_reflectance.py +++ b/pyspectral/tests/test_reflectance.py @@ -203,13 +203,25 @@ def test_reflectance(self): refl = refl37.reflectance_from_tbs(sunz, tb3, tb4) self.assertTrue(hasattr(refl, 'mask')) + sunz = np.array([80.], dtype=np.float32) + tb3 = np.array([295.], dtype=np.float32) + tb4 = np.array([282.], dtype=np.float32) + refl = refl37.reflectance_from_tbs(sunz, tb3, tb4) + np.testing.assert_allclose(refl, np.array([0.45249779], np.float32)) + assert refl.dtype == np.float32 + try: import dask.array as da - sunz = da.from_array([50.], chunks=10) - tb3 = da.from_array([300.], chunks=10) - tb4 = da.from_array([285.], chunks=10) - refl = refl37.reflectance_from_tbs(sunz, tb3, tb4) - self.assertTrue(hasattr(refl, 'compute')) + import dask.config + + from pyspectral.tests.unittest_helpers import ComputeCountingScheduler + + with dask.config.set(scheduler=ComputeCountingScheduler(max_computes=0)): + sunz = da.from_array([50.], chunks=10) + tb3 = da.from_array([300.], chunks=10) + tb4 = da.from_array([285.], chunks=10) + refl = refl37.reflectance_from_tbs(sunz, tb3, tb4) + self.assertTrue(hasattr(refl, 'compute')) except ImportError: pass diff --git a/pyspectral/tests/unittest_helpers.py b/pyspectral/tests/unittest_helpers.py index 13d394c..9594a0c 100644 --- a/pyspectral/tests/unittest_helpers.py +++ b/pyspectral/tests/unittest_helpers.py @@ -28,3 +28,21 @@ def assertNumpyArraysEqual(self, other): raise AssertionError("Shapes don't match") if not np.allclose(self, other): raise AssertionError("Elements don't match!") + + +class ComputeCountingScheduler: + """Scheduler raising an exception if data are computed too many times.""" + + def __init__(self, max_computes=1): + """Set starting and maximum compute counts.""" + self.max_computes = max_computes + self.total_computes = 0 + + def __call__(self, dsk, keys, **kwargs): + """Compute dask task and keep track of number of times we do so.""" + import dask + self.total_computes += 1 + if self.total_computes > self.max_computes: + raise RuntimeError("Too many dask computations were scheduled: " + "{}".format(self.total_computes)) + return dask.get(dsk, keys, **kwargs)