Merge branch 'main' into fix-nir-reflectance-dtype

pyspectral/blackbody.py

@@ -28,96 +28,50 @@
 except ImportError:
     da = np
 
+from pyspectral.utils import use_map_blocks_on
+
 LOG = logging.getLogger(__name__)
 
 H_PLANCK = 6.62606957 * 1e-34  # SI-unit = [J*s]
 K_BOLTZMANN = 1.3806488 * 1e-23  # SI-unit = [J/K]
 C_SPEED = 2.99792458 * 1e8  # SI-unit = [m/s]
 
+PLANCK_C1 = H_PLANCK * C_SPEED / K_BOLTZMANN
+PLANCK_C2 = 2 * H_PLANCK * C_SPEED**2
+
 EPSILON = 0.000001
 
 
+@use_map_blocks_on("radiance")
 def blackbody_rad2temp(wavelength, radiance):
-    """Derive brightness temperatures from radiance using the Planck function.
-
-    Wavelength space. Assumes SI units as input and returns
-    temperature in Kelvin
-
-    """
-    mask = False
-    if np.isscalar(radiance):
-        rad = np.array([radiance, ], dtype='float64')
-    else:
-        rad = np.array(radiance, dtype='float64')
-        if np.ma.is_masked(radiance):
-            mask = radiance.mask
-
-    rad = np.ma.masked_array(rad, mask=mask)
-    rad = np.ma.masked_less_equal(rad, 0)
+    """Derive brightness temperatures from radiance using the inverse Planck function.
 
-    if np.isscalar(wavelength):
-        wvl = np.array([wavelength, ], dtype='float64')
-    else:
-        wvl = np.array(wavelength, dtype='float64')
+    Args:
+        wavelength: The wavelength to use, in SI units (metre).
+        radiance: The radiance to derive temperatures from, in SI units (W/m^2 sr^-1). Scalar or arrays are accepted.
 
-    const1 = H_PLANCK * C_SPEED / K_BOLTZMANN
-    const2 = 2 * H_PLANCK * C_SPEED**2
-    res = const1 / (wvl * np.log(np.divide(const2, rad * wvl**5) + 1.0))
+    Returns:
+        The derived temperature in Kelvin.
 
-    shape = rad.shape
-    resshape = res.shape
-
-    if wvl.shape[0] == 1:
-        if rad.shape[0] == 1:
-            return res[0]
-        else:
-            return res[::].reshape(shape)
-    else:
-        if rad.shape[0] == 1:
-            return res[0, :]
-        else:
-            if len(shape) == 1:
-                return np.reshape(res, (shape[0], resshape[1]))
-            else:
-                return np.reshape(res, (shape[0], shape[1], resshape[1]))
+    """
+    with np.errstate(invalid='ignore'):
+        return PLANCK_C1 / (wavelength * np.log(PLANCK_C2 / (radiance * wavelength**5) + 1.0))
 
 
+@use_map_blocks_on("radiance")
 def blackbody_wn_rad2temp(wavenumber, radiance):
-    """Derive brightness temperatures from radiance using the Planck function.
+    """Derive brightness temperatures from radiance using the inverse Planck function.
+
+    Args:
+        wavenumber: The wavenumber to use, in SI units (1/metre).
+        radiance: The radiance to derive temperatures from, in SI units (W/m^2 sr^-1). Scalar or arrays are accepted.
 
-    Wavenumber space
+    Returns:
+        The derived temperature in Kelvin.
 
     """
-    if np.isscalar(radiance):
-        radiance = np.array([radiance], dtype='float64')
-    elif isinstance(radiance, (list, tuple)):
-        radiance = np.array(radiance, dtype='float64')
-    if np.isscalar(wavenumber):
-        wavnum = np.array([wavenumber], dtype='float64')
-    elif isinstance(wavenumber, (list, tuple)):
-        wavnum = np.array(wavenumber, dtype='float64')
-
-    const1 = H_PLANCK * C_SPEED / K_BOLTZMANN
-    const2 = 2 * H_PLANCK * C_SPEED**2
-    res = const1 * wavnum / np.log(
-        np.divide(const2 * wavnum**3, radiance) + 1.0)
-
-    shape = radiance.shape
-    resshape = res.shape
-
-    if wavnum.shape[0] == 1:
-        if radiance.shape[0] == 1:
-            return res[0]
-        else:
-            return res[::].reshape(shape)
-    else:
-        if radiance.shape[0] == 1:
-            return res[0, :]
-        else:
-            if len(shape) == 1:
-                return res.reshape((shape[0], resshape[1]))
-            else:
-                return res.reshape((shape[0], shape[1], resshape[1]))
+    with np.errstate(invalid='ignore'):
+        return PLANCK_C1 * wavenumber / np.log((PLANCK_C2 * wavenumber**3) / radiance + 1.0)
 
 
 def planck(wave, temperature, wavelength=True):

pyspectral/tests/test_blackbody.py

@@ -107,7 +107,7 @@ def test_blackbody_wn(self):
         temp2 = blackbody_wn_rad2temp(wavenumber, black[1])
         self.assertAlmostEqual(temp2, 301.0, 4)
 
-        t__ = blackbody_wn_rad2temp(wavenumber, [0.001, 0.0009])
+        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])

pyspectral/tests/test_rad_tb_conversions.py

@@ -19,23 +19,16 @@
 
 """Testing the radiance to brightness temperature conversion."""
 
-import sys
+import unittest
+import warnings
+from unittest.mock import patch
 
 import numpy as np
+import pytest
 
-from pyspectral.radiance_tb_conversion import RadTbConverter, SeviriRadTbConverter
+from pyspectral.radiance_tb_conversion import RadTbConverter, SeviriRadTbConverter, radiance2tb
 from pyspectral.utils import get_central_wave
 
-if sys.version_info < (2, 7):
-    import unittest2 as unittest
-else:
-    import unittest
-if sys.version_info < (3,):
-    from mock import patch
-else:
-    from unittest.mock import patch
-
-
 TEST_TBS = np.array([200., 270., 300., 302., 350.], dtype='float32')
 
 TRUE_RADS = np.array([856.937353205, 117420.385297,
@@ -327,33 +320,106 @@ def test_get_bandname(self, download_rsr, load, isfile, exists):
     def test_rad2tb(self):
         """Unit testing the radiance to brightness temperature conversion."""
         res = self.modis.tb2radiance(TEST_TBS, lut=False)
-        self.assertTrue(np.allclose(TRUE_RADS, res['radiance']))
+        np.testing.assert_allclose(TRUE_RADS, res['radiance'], rtol=1e-5, atol=1e-8)
 
         res = self.modis2.tb2radiance(TEST_TBS, lut=False)
-        self.assertTrue(np.allclose(TRUE_RADS, res['radiance']))
+        np.testing.assert_allclose(TRUE_RADS, res['radiance'], rtol=1e-5, atol=1e-8)
 
         rad = res['radiance']
         tbs = self.modis.radiance2tb(rad)
-        self.assertTrue(np.allclose(TEST_TBS, tbs, atol=0.25))
+        np.testing.assert_allclose(TEST_TBS, tbs, atol=0.25)
 
         res = self.modis.tb2radiance(TEST_TBS, lut=False, normalized=False)
         integral = self.modis.rsr_integral
-        self.assertTrue(np.allclose(TRUE_RADS * integral, res['radiance']))
+        np.testing.assert_allclose(TRUE_RADS * integral, res['radiance'], rtol=1e-5, atol=1e-8)
 
         res = self.modis.tb2radiance(237., lut=False)
-        self.assertAlmostEqual(16570.579551068, res['radiance'])
+        assert res['radiance'] == pytest.approx(16570.579551068)
 
         res = self.modis.tb2radiance(277., lut=False)
-        self.assertAlmostEqual(167544.39368663222, res['radiance'])
+        assert res['radiance'] == pytest.approx(167544.39368663222)
 
         res = self.modis.tb2radiance(1.1, lut=False)
-        self.assertAlmostEqual(0.0, res['radiance'])
+        assert res['radiance'] == pytest.approx(0.0)
 
         res = self.modis.tb2radiance(11.1, lut=False)
-        self.assertAlmostEqual(0.0, res['radiance'])
+        assert res['radiance'] == pytest.approx(0.0)
 
         res = self.modis.tb2radiance(100.1, lut=False)
-        self.assertAlmostEqual(5.3940515573e-06, res['radiance'])
+        assert res['radiance'] == pytest.approx(5.3940515573e-06, abs=1e-7)
 
         res = self.modis.tb2radiance(200.1, lut=False)
-        self.assertAlmostEqual(865.09759706, res['radiance'])
+        assert res['radiance'] == pytest.approx(865.09759706)
+
+
+def test_rad2tb_types():
+    """Test radiance to brightness temperature conversion preserves shape and type."""
+    # 1d rads
+    rad = TRUE_RADS
+    central_wavelength = TEST_RSR['20']['det-1']['central_wavelength'] * 1e-6
+    tbs = radiance2tb(rad, central_wavelength)
+    np.testing.assert_allclose(tbs, TEST_TBS, atol=0.25)
+
+    # 2d rads
+    rad_2d = np.vstack((TRUE_RADS, TRUE_RADS))
+    tbs = radiance2tb(rad_2d, central_wavelength)
+    assert tbs.shape == rad_2d.shape
+
+    tbs = radiance2tb(rad_2d.astype(np.float32), central_wavelength)
+    assert tbs.dtype == np.float32
+
+
+def test_rad2tb_xarray_types():
+    """Test radiance to brightness temperature conversion works with xarray."""
+    xr = pytest.importorskip("xarray")
+
+    central_wavelength = 3.78e-6
+
+    rad_2d = np.vstack((TRUE_RADS, TRUE_RADS))
+    rad = xr.DataArray(rad_2d, dims=["y", "x"])
+    tbs = radiance2tb(rad, central_wavelength)
+    assert isinstance(tbs, xr.DataArray)
+
+    rad = xr.DataArray(rad_2d.astype(np.float32), dims=["y", "x"])
+    tbs = radiance2tb(rad, central_wavelength)
+    assert tbs.dtype == np.float32
+
+
+def test_rad2tb_does_not_warn_about_invalid_log_values():
+    """Test that radiance to temp does not warn about invalid log values."""
+    xr = pytest.importorskip("xarray")
+
+    central_wavelength = 3.78e-6
+
+    rad = xr.DataArray(np.array([-1, 2, 3]).astype(np.float32), dims=["x"])
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        tbs = radiance2tb(rad, central_wavelength)
+        assert np.isnan(tbs[0])
+
+    da = pytest.importorskip("dask.array")
+
+    rad = da.from_array(np.array([-1, 2, 3]), chunks=2).astype(np.float32)
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        tbs = da.map_blocks(radiance2tb, rad, wavelength=central_wavelength)
+        assert np.isnan(tbs[0])
+
+    rad = xr.DataArray(da.from_array(np.array([-1, 2, 3]), chunks=2).astype(np.float32), dims=["x"])
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        tbs = xr.map_blocks(radiance2tb, rad, kwargs=dict(wavelength=central_wavelength))
+        assert np.isnan(tbs[0])
+
+    rad = da.from_array(np.array([-1, 2, 3]), chunks=2).astype(np.float32)
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        tbs = radiance2tb(rad, central_wavelength)
+        assert np.isnan(tbs[0])
+
+    rad = xr.DataArray(da.from_array(np.array([-1, 2, 3]), chunks=2).astype(np.float32), dims=["x"])
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        tbs = radiance2tb(rad, central_wavelength)
+        assert np.isnan(tbs[0])
+        assert isinstance(tbs, xr.DataArray)

pyspectral/utils.py

@@ -23,6 +23,8 @@
 import os
 import tarfile
 import warnings
+from functools import wraps
+from inspect import getfullargspec
 
 import numpy as np
 import requests
@@ -580,3 +582,35 @@ def are_instruments_identical(name1, name2):
     if name1 == INSTRUMENT_TRANSLATION_DASH2SLASH.get(name2):
         return True
     return False
+
+
+def use_map_blocks_on(argument_to_run_map_blocks_on):
+    """Use map blocks on a given argument.
+
+    This decorator assumes only one of the arguments of the decorated function is chunked.
+    """
+    def decorator(f):
+        argspec = getfullargspec(f)
+        argument_index = argspec.args.index(argument_to_run_map_blocks_on)
+
+        @wraps(f)
+        def wrapper(*args, **kwargs):
+            array = args[argument_index]
+            chunks = getattr(array, "chunks", None)
+            if chunks is None:
+                return f(*args, **kwargs)
+            import dask.array as da
+            import xarray as xr
+            if isinstance(array, da.Array):
+                return da.map_blocks(f, *args, **kwargs)
+            elif isinstance(array, xr.DataArray):
+                new_array = array.copy()
+                new_args = list(args)
+                new_args[argument_index] = array.data
+                new_data = da.map_blocks(f, *new_args, **kwargs)
+                new_array.data = new_data
+                return new_array
+            else:
+                raise NotImplementedError(f"Don't know how to map_blocks on {type(array)}")
+        return wrapper
+    return decorator