Return 3.x um reflectance with the same dtype as the NIR data

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):

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

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):

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,108 +37,76 @@
 __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."""
         import dask
         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)

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):