From fdfeea04c69832fe96fc5ceea04e31cc8d638411 Mon Sep 17 00:00:00 2001 From: Martin Raspaud Date: Fri, 17 Nov 2023 08:47:07 +0100 Subject: [PATCH 1/4] Allow dataarrays and preserve dtype in rad2temp --- pyspectral/blackbody.py | 82 ++++----------------- pyspectral/tests/test_blackbody.py | 2 +- pyspectral/tests/test_rad_tb_conversions.py | 56 +++++++++----- 3 files changed, 54 insertions(+), 86 deletions(-) diff --git a/pyspectral/blackbody.py b/pyspectral/blackbody.py index b2af8725..bc02e4eb 100644 --- a/pyspectral/blackbody.py +++ b/pyspectral/blackbody.py @@ -38,86 +38,36 @@ def blackbody_rad2temp(wavelength, radiance): - """Derive brightness temperatures from radiance using the Planck function. + """Derive brightness temperatures from radiance using the inverse Planck function. - Wavelength space. Assumes SI units as input and returns - temperature in Kelvin + 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. - """ - 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) - - if np.isscalar(wavelength): - wvl = np.array([wavelength, ], dtype='float64') - else: - wvl = np.array(wavelength, dtype='float64') + Returns: + The derived temperature in Kelvin. + """ 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)) + return const1 / (wavelength * np.log(const2 / (radiance * wavelength**5) + 1.0)) - 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])) 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. - Wavenumber space + Args: + wavenumber: The wavenumber to use, in SI units. + radiance: The radiance to derive temperatures from, in SI units. Scalar or arrays are accepted. - """ - 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') + Returns: + The derived temperature in Kelvin. + """ 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])) + return const1 * wavenumber / np.log((const2 * wavenumber**3) / radiance + 1.0) def planck(wave, temperature, wavelength=True): diff --git a/pyspectral/tests/test_blackbody.py b/pyspectral/tests/test_blackbody.py index e9341ced..a40e2d5a 100644 --- a/pyspectral/tests/test_blackbody.py +++ b/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]) diff --git a/pyspectral/tests/test_rad_tb_conversions.py b/pyspectral/tests/test_rad_tb_conversions.py index 2a8c41d6..8628f2ec 100644 --- a/pyspectral/tests/test_rad_tb_conversions.py +++ b/pyspectral/tests/test_rad_tb_conversions.py @@ -22,18 +22,14 @@ import sys 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 +import unittest + +from unittest.mock import patch TEST_TBS = np.array([200., 270., 300., 302., 350.], dtype='float32') @@ -327,33 +323,55 @@ 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_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)) + tbs = radiance2tb(rad, central_wavelength) + assert tbs.shape == rad.shape + + xr = pytest.importorskip("xarray") + rad = xr.DataArray(TRUE_RADS, dims=["x"]) + tbs = radiance2tb(rad, central_wavelength) + assert isinstance(tbs, xr.DataArray) + + rad = xr.DataArray(TRUE_RADS.astype(np.float32), dims=["x"]) + tbs = radiance2tb(rad, central_wavelength) + assert tbs.dtype == np.float32 From 14945c833caacb241a37e5d3e66bad1bb6d81c60 Mon Sep 17 00:00:00 2001 From: Martin Raspaud Date: Fri, 17 Nov 2023 08:50:57 +0100 Subject: [PATCH 2/4] Fix linting errors --- pyspectral/blackbody.py | 1 - pyspectral/tests/test_rad_tb_conversions.py | 8 ++------ 2 files changed, 2 insertions(+), 7 deletions(-) diff --git a/pyspectral/blackbody.py b/pyspectral/blackbody.py index bc02e4eb..4c607121 100644 --- a/pyspectral/blackbody.py +++ b/pyspectral/blackbody.py @@ -53,7 +53,6 @@ def blackbody_rad2temp(wavelength, radiance): return const1 / (wavelength * np.log(const2 / (radiance * wavelength**5) + 1.0)) - def blackbody_wn_rad2temp(wavenumber, radiance): """Derive brightness temperatures from radiance using the inverse Planck function. diff --git a/pyspectral/tests/test_rad_tb_conversions.py b/pyspectral/tests/test_rad_tb_conversions.py index 8628f2ec..101b52fa 100644 --- a/pyspectral/tests/test_rad_tb_conversions.py +++ b/pyspectral/tests/test_rad_tb_conversions.py @@ -19,7 +19,8 @@ """Testing the radiance to brightness temperature conversion.""" -import sys +import unittest +from unittest.mock import patch import numpy as np import pytest @@ -27,11 +28,6 @@ from pyspectral.radiance_tb_conversion import RadTbConverter, SeviriRadTbConverter, radiance2tb from pyspectral.utils import get_central_wave -import unittest - -from unittest.mock import patch - - TEST_TBS = np.array([200., 270., 300., 302., 350.], dtype='float32') TRUE_RADS = np.array([856.937353205, 117420.385297, From 900fc7d941c794d3bcca6d728e01e7688b7a8bd8 Mon Sep 17 00:00:00 2001 From: Martin Raspaud Date: Fri, 17 Nov 2023 17:00:53 +0100 Subject: [PATCH 3/4] Suppress log warnings for inverse planck functions --- pyspectral/blackbody.py | 23 +++--- pyspectral/tests/test_rad_tb_conversions.py | 82 +++++++++++++++++---- pyspectral/utils.py | 31 ++++++++ 3 files changed, 112 insertions(+), 24 deletions(-) diff --git a/pyspectral/blackbody.py b/pyspectral/blackbody.py index 4c607121..ed9faf1a 100644 --- a/pyspectral/blackbody.py +++ b/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): diff --git a/pyspectral/tests/test_rad_tb_conversions.py b/pyspectral/tests/test_rad_tb_conversions.py index 101b52fa..0ed6aa0d 100644 --- a/pyspectral/tests/test_rad_tb_conversions.py +++ b/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) diff --git a/pyspectral/utils.py b/pyspectral/utils.py index 6b6cccb6..edc825bb 100644 --- a/pyspectral/utils.py +++ b/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 From ec59df798105210a0c183136ec160cb5850a0155 Mon Sep 17 00:00:00 2001 From: Martin Raspaud Date: Mon, 20 Nov 2023 09:07:45 +0100 Subject: [PATCH 4/4] Improve docstring --- pyspectral/utils.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/pyspectral/utils.py b/pyspectral/utils.py index edc825bb..02f92bf5 100644 --- a/pyspectral/utils.py +++ b/pyspectral/utils.py @@ -585,7 +585,10 @@ def are_instruments_identical(name1, name2): def use_map_blocks_on(argument_to_run_map_blocks_on): - """Use map blocks on a given argument.""" + """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)