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

pyspectral/blackbody.py

@@ -99,24 +99,17 @@ 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')
-
+    temperature = _scalar_tuple_list_to_numpy(temperature)
     shape = temperature.shape
-    if np.isscalar(wave):
-        wln = np.array([wave, ], dtype='float64')
-    else:
-        wln = np.array(wave, dtype='float64')
+    wave = _scalar_tuple_list_to_numpy(wave)
 
     if wavelength:
         const = 2 * H_PLANCK * C_SPEED ** 2
-        nom = const / wln ** 5
-        arg1 = H_PLANCK * C_SPEED / (K_BOLTZMANN * wln)
+        nom = const / wave ** 5
+        arg1 = H_PLANCK * C_SPEED / (K_BOLTZMANN * wave)
     else:
-        nom = 2 * H_PLANCK * (C_SPEED ** 2) * (wln ** 3)
-        arg1 = H_PLANCK * C_SPEED * wln / K_BOLTZMANN
+        nom = 2 * H_PLANCK * (C_SPEED ** 2) * (wave ** 3)
+        arg1 = H_PLANCK * C_SPEED * wave / K_BOLTZMANN
 
     with np.errstate(divide='ignore', invalid='ignore'):
         # use dask functions when needed
@@ -132,7 +125,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"
@@ -153,7 +146,7 @@ def planck(wave, temperature, wavelength=True):
         denom = np.exp(exp_arg) - 1
         rad = nom / denom
         radshape = rad.shape
-        if wln.shape[0] == 1:
+        if wave.shape[0] == 1:
             if temperature.shape[0] == 1:
                 return rad[0, 0]
             else:
@@ -168,6 +161,14 @@ def planck(wave, temperature, wavelength=True):
                     return rad.reshape((shape[0], shape[1], radshape[1]))
 
 
+def _scalar_tuple_list_to_numpy(val):
+    if np.isscalar(val):
+        return np.array([val, ], dtype='float64')
+    if isinstance(val, (list, tuple)):
+        val = np.array(val, dtype='float64')
+    return val
+
+
 def blackbody_wn(wavenumber, temp):
     """Derive the Planck radiation as a function of wavenumber.
 

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 @@
             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,12 @@
 
 """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
 
 RAD_11MICRON_300KELVIN = 9573176.935507433
 RAD_11MICRON_301KELVIN = 9714686.576498277
@@ -53,75 +53,78 @@ def __call__(self, dsk, keys, **kwargs):
         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)
+        assert black.shape[0] == 2
+        np.testing.assert_almost_equal(black[0], RAD_11MICRON_300KELVIN)
+        np.testing.assert_almost_equal(black[1], RAD_11MICRON_301KELVIN)
 
         temp1 = blackbody_rad2temp(wavel, black[0])
-        self.assertAlmostEqual(temp1, 300.0, 4)
+        np.testing.assert_almost_equal(temp1, 300.0, 4)
         temp2 = blackbody_rad2temp(wavel, black[1])
-        self.assertAlmostEqual(temp2, 301.0, 4)
+        np.testing.assert_almost_equal(temp2, 301.0, 4)
 
         black = blackbody(13. * 1E-6, 200.)
-        self.assertTrue(np.isscalar(black))
+        assert np.isscalar(black)
 
         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)
+        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)
+        assert isinstance(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)
+        tb_therm = da.array([[300., 301], [299, 298], [279, 286]])
         with dask.config.set(scheduler=CustomScheduler(0)):
             black = blackbody((10. * 1E-6, 11.e-6), tb_therm)
-        self.assertIsInstance(black, da.Array)
+        assert isinstance(black, da.Array)
+        assert black.dtype == np.float64
 
-        tb_therm = da.from_array([[300., 301], [0., 298], [279, 286]], chunks=2)
+    @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=CustomScheduler(0)):
-            black = blackbody((10. * 1E-6, 11.e-6), tb_therm)
-        self.assertIsInstance(black, da.Array)
+            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)
+        assert black.shape[0] == 2
+        np.testing.assert_almost_equal(black[0], WN_RAD_11MICRON_300KELVIN)
+        np.testing.assert_almost_equal(black[1], WN_RAD_11MICRON_301KELVIN)
 
         temp1 = blackbody_wn_rad2temp(wavenumber, black[0])
-        self.assertAlmostEqual(temp1, 300.0, 4)
+        np.testing.assert_almost_equal(temp1, 300.0, 4)
         temp2 = blackbody_wn_rad2temp(wavenumber, black[1])
-        self.assertAlmostEqual(temp2, 301.0, 4)
+        np.testing.assert_almost_equal(temp2, 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_almost_equal(t__[0], expected[0])
+        np.testing.assert_almost_equal(t__[1], expected[1])
 
         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)
+        np.testing.assert_almost_equal(t__[1, 1], expected[1, 1], 5)
+        np.testing.assert_almost_equal(t__[0, 0], expected[0, 0], 5)
+        np.testing.assert_almost_equal(t__[0, 1], expected[0, 1], 5)
+        np.testing.assert_almost_equal(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."""
@@ -131,13 +134,13 @@ def test_blackbody_wn_dask(self):
         radiances = da.from_array([0.001, 0.0009, 0.0012, 0.0018], chunks=2).reshape(2, 2)
         with dask.config.set(scheduler=CustomScheduler(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)
+        np.testing.assert_almost_equal(t__[1, 1], expected[1, 1], 5)
+        np.testing.assert_almost_equal(t__[0, 0], expected[0, 0], 5)
+        np.testing.assert_almost_equal(t__[0, 1], expected[0, 1], 5)
+        np.testing.assert_almost_equal(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):

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 CustomScheduler
+
+            with dask.config.set(scheduler=CustomScheduler(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