Suppress log warnings for inverse planck functions

pyspectral/blackbody.py

@@ -28,45 +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 inverse Planck function.
 
     Args:
         wavelength: The wavelength to use, in SI units (metre).
-        radiance: The radiance to derive temperatures from, in SI units. Scalar or arrays are accepted.
+        radiance: The radiance to derive temperatures from, in SI units (W/m^2 sr^-1). Scalar or arrays are accepted.
 
     Returns:
         The derived temperature in Kelvin.
 
     """
-    const1 = H_PLANCK * C_SPEED / K_BOLTZMANN
-    const2 = 2 * H_PLANCK * C_SPEED**2
-    return const1 / (wavelength * np.log(const2 / (radiance * wavelength**5) + 1.0))
+    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 inverse Planck function.
 
     Args:
-        wavenumber: The wavenumber to use, in SI units.
-        radiance: The radiance to derive temperatures from, in SI units. Scalar or arrays are accepted.
+        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.
 
     Returns:
         The derived temperature in Kelvin.
 
     """
-    const1 = H_PLANCK * C_SPEED / K_BOLTZMANN
-    const2 = 2 * H_PLANCK * C_SPEED**2
-    return const1 * wavenumber / np.log((const2 * wavenumber**3) / radiance + 1.0)
+    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_rad_tb_conversions.py

@@ -20,6 +20,7 @@
 """Testing the radiance to brightness temperature conversion."""
 
 import unittest
+import warnings
 from unittest.mock import patch
 
 import numpy as np
@@ -350,24 +351,75 @@ def test_rad2tb(self):
         res = self.modis.tb2radiance(200.1, lut=False)
         assert res['radiance'] == pytest.approx(865.09759706)
 
-    def test_rad2tb_arrays_and_types(self):
-        """Unit testing the radiance to brightness temperature conversion."""
-        # 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 = np.vstack((TRUE_RADS, TRUE_RADS))
+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 tbs.shape == rad.shape
+        assert np.isnan(tbs[0])
+
+    da = pytest.importorskip("dask.array")
 
-        xr = pytest.importorskip("xarray")
-        rad = xr.DataArray(TRUE_RADS, dims=["x"])
+    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 isinstance(tbs, xr.DataArray)
+        assert np.isnan(tbs[0])
 
-        rad = xr.DataArray(TRUE_RADS.astype(np.float32), dims=["x"])
+    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 tbs.dtype == np.float32
+        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,32 @@ 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."""
+    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