[Git][debian-gis-team/pyspectral][upstream] New upstream version 0.13.0+ds
Antonio Valentino (@antonio.valentino)
gitlab at salsa.debian.org
Sat Dec 2 20:49:29 GMT 2023
Antonio Valentino pushed to branch upstream at Debian GIS Project / pyspectral
Commits:
14cabd37 by Antonio Valentino at 2023-12-02T15:18:46+00:00
New upstream version 0.13.0+ds
- - - - -
22 changed files:
- CHANGELOG.md
- CHANGELOG_RSR_DATA.rst
- continuous_integration/environment.yaml
- doc/37_reflectance.rst
- doc/rad_definitions.rst
- doc/usage.rst
- pyspectral/bandnames.py
- pyspectral/blackbody.py
- pyspectral/config.py
- pyspectral/etc/pyspectral.yaml
- pyspectral/near_infrared_reflectance.py
- pyspectral/radiance_tb_conversion.py
- pyspectral/solar.py
- pyspectral/tests/test_blackbody.py
- pyspectral/tests/test_rad_tb_conversions.py
- pyspectral/tests/test_rayleigh.py
- pyspectral/tests/test_reflectance.py
- pyspectral/tests/test_solarflux.py
- pyspectral/tests/unittest_helpers.py
- pyspectral/utils.py
- + rsr_convert_scripts/epic_reader.py
- setup.py
Changes:
=====================================
CHANGELOG.md
=====================================
@@ -1,3 +1,26 @@
+## Version <v0.13.0> (2023/11/27)
+
+### Issues Closed
+
+* [Issue 205](https://github.com/pytroll/pyspectral/issues/205) - Near infrared reflectance calculation does not preserve input dtype ([PR 206](https://github.com/pytroll/pyspectral/pull/206) by [@pnuu](https://github.com/pnuu))
+* [Issue 201](https://github.com/pytroll/pyspectral/issues/201) - blackbody functions not dask friendly, trigger early dask computation
+
+In this release 2 issues were closed.
+
+### Pull Requests Merged
+
+#### Features added
+
+* [PR 208](https://github.com/pytroll/pyspectral/pull/208) - Update the solar spectrum code to use the included reference spectra
+* [PR 207](https://github.com/pytroll/pyspectral/pull/207) - Remove pkg_resources usage in config handling
+* [PR 206](https://github.com/pytroll/pyspectral/pull/206) - Return 3.x um reflectance with the same dtype as the NIR data ([205](https://github.com/pytroll/pyspectral/issues/205))
+* [PR 204](https://github.com/pytroll/pyspectral/pull/204) - Bump up python versions
+* [PR 203](https://github.com/pytroll/pyspectral/pull/203) - Allow dataarrays and preserve dtype in rad2temp
+* [PR 202](https://github.com/pytroll/pyspectral/pull/202) - Update the solar spectrum code to use the included reference spectra
+
+In this release 6 pull requests were closed.
+
+
## Version <v0.12.5> (2023/09/21)
### Pull Requests Merged
=====================================
CHANGELOG_RSR_DATA.rst
=====================================
@@ -2,6 +2,20 @@ Changelog for the Relative Spectral Response data
=================================================
+
+Version <v1.2.4> (Sat Oct 21 12:25:00 PM CEST 2023)
+-------------------------------------------
+
+ * Normalized RSR responses for the EPIC sensor. Now values are sclaed to be
+ between 0 and 1.
+
+
+Version <v1.2.3> (Fri Oct 20 01:58:29 2023)
+-------------------------------------------
+
+ * Added RSR file for EPIC on DSCOVR (responses not normalized)
+
+
Version <v1.2.2> (Tue Nov 15 14:44:03 2022)
-------------------------------------------
=====================================
continuous_integration/environment.yaml
=====================================
@@ -17,8 +17,6 @@ dependencies:
- behave
- netcdf4
- h5py
- - h5netcdf
- - pyhdf
- mock
- pytest
- pytest-cov
=====================================
doc/37_reflectance.rst
=====================================
@@ -123,9 +123,9 @@ would be:
where :math:`S(\lambda)` is the spectral solar irradiance.
>>> from pyspectral.rsr_reader import RelativeSpectralResponse
- >>> from pyspectral.solar import (SolarIrradianceSpectrum, TOTAL_IRRADIANCE_SPECTRUM_2000ASTM)
+ >>> from pyspectral.solar import SolarIrradianceSpectrum
>>> viirs = RelativeSpectralResponse('Suomi-NPP', 'viirs')
- >>> solar_irr = SolarIrradianceSpectrum(TOTAL_IRRADIANCE_SPECTRUM_2000ASTM, dlambda=0.005)
+ >>> solar_irr = SolarIrradianceSpectrum(dlambda=0.005)
>>> sflux = solar_irr.inband_solarflux(viirs.rsr['M12'])
>>> print(np.round(sflux, 7))
2.2428119
=====================================
doc/rad_definitions.rst
=====================================
@@ -156,8 +156,8 @@ TOA Solar irradiance and solar constant
First, the TOA solar irradiance in wavelength space:
- >>> from pyspectral.solar import (SolarIrradianceSpectrum, TOTAL_IRRADIANCE_SPECTRUM_2000ASTM)
- >>> solar_irr = SolarIrradianceSpectrum(TOTAL_IRRADIANCE_SPECTRUM_2000ASTM, dlambda=0.0005)
+ >>> from pyspectral.solar import SolarIrradianceSpectrum
+ >>> solar_irr = SolarIrradianceSpectrum(dlambda=0.0005)
>>> print("Solar irradiance = {}".format(round(solar_irr.solar_constant(), 3)))
Solar irradiance = 1366.091
>>> solar_irr.plot('/tmp/solar_irradiance.png')
@@ -201,8 +201,8 @@ In python code it may look like this:
>>> from pyspectral.utils import convert2wavenumber, get_central_wave
>>> seviri = RelativeSpectralResponse('Meteosat-8', 'seviri')
>>> rsr, info = convert2wavenumber(seviri.rsr)
- >>> from pyspectral.solar import (SolarIrradianceSpectrum, TOTAL_IRRADIANCE_SPECTRUM_2000ASTM)
- >>> solar_irr = SolarIrradianceSpectrum(TOTAL_IRRADIANCE_SPECTRUM_2000ASTM, dlambda=0.0005, wavespace='wavenumber')
+ >>> from pyspectral.solar import SolarIrradianceSpectrum
+ >>> solar_irr = SolarIrradianceSpectrum(dlambda=0.0005, wavespace='wavenumber')
>>> print("Solar Irradiance (SEVIRI band VIS008) = {sflux:12.6f}".format(sflux=solar_irr.inband_solarflux(rsr['VIS0.8'])))
Solar Irradiance (SEVIRI band VIS008) = 63767.908405
=====================================
doc/usage.rst
=====================================
@@ -35,9 +35,9 @@ Now, you can work with the data as you wish, make some simple plot for instance:
A simple use case:
>>> from pyspectral.rsr_reader import RelativeSpectralResponse
- >>> from pyspectral.solar import (SolarIrradianceSpectrum, TOTAL_IRRADIANCE_SPECTRUM_2000ASTM)
+ >>> from pyspectral.solar import SolarIrradianceSpectrum
>>> modis = RelativeSpectralResponse('EOS-Aqua', 'modis')
- >>> solar_irr = SolarIrradianceSpectrum(TOTAL_IRRADIANCE_SPECTRUM_2000ASTM, dlambda=0.005)
+ >>> solar_irr = SolarIrradianceSpectrum(dlambda=0.005)
>>> sflux = solar_irr.inband_solarflux(modis.rsr['20'])
>>> print("Solar flux over Band: {sflux}".format(sflux=round(sflux, 6)))
Solar flux over Band: 2.002928
=====================================
pyspectral/bandnames.py
=====================================
@@ -1,7 +1,7 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
-# Copyright (c) 2021-2022 Adam.Dybbroe
+# Copyright (c) 2021-2023 Pytroll developers
#
#
# This program is free software: you can redistribute it and/or modify
@@ -70,7 +70,6 @@ for chan_num in range(1, 37):
# MODIS RSR files were made before 'chX' became standard in pyspectral
BANDNAMES['modis'] = {str(chan_num): str(chan_num) for chan_num in range(1, 37)}
-
BANDNAMES['seviri'] = {'VIS006': 'VIS0.6',
'VIS008': 'VIS0.8',
'IR_016': 'NIR1.6',
=====================================
pyspectral/blackbody.py
=====================================
@@ -28,96 +28,50 @@ try:
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
+ at 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)
-
- if np.isscalar(wavelength):
- wvl = np.array([wavelength, ], dtype='float64')
- else:
- wvl = np.array(wavelength, dtype='float64')
+ """Derive brightness temperatures from radiance using the inverse Planck function.
- 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))
+ 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.
- shape = rad.shape
- resshape = res.shape
+ Returns:
+ The derived temperature in Kelvin.
- 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))
+ at 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):
@@ -145,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
@@ -178,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"
@@ -197,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/config.py
=====================================
@@ -21,20 +21,13 @@
import logging
import os
+from collections.abc import Mapping
from os.path import expanduser
+from pathlib import Path
import yaml
from appdirs import AppDirs
-try:
- # python 3.3+
- from collections.abc import Mapping
-except ImportError:
- # deprecated (above can't be done in 2.7)
- from collections import Mapping
-
-import pkg_resources
-
try:
from yaml import UnsafeLoader
except ImportError:
@@ -43,16 +36,7 @@ except ImportError:
LOG = logging.getLogger(__name__)
-BUILTIN_CONFIG_FILE = pkg_resources.resource_filename(__name__,
- os.path.join('etc', 'pyspectral.yaml'))
-
-
-CONFIG_FILE = os.environ.get('PSP_CONFIG_FILE')
-
-if CONFIG_FILE is not None and (not os.path.exists(CONFIG_FILE) or
- not os.path.isfile(CONFIG_FILE)):
- raise IOError(str(CONFIG_FILE) + " pointed to by the environment " +
- "variable PSP_CONFIG_FILE is not a file or does not exist!")
+BUILTIN_CONFIG_FILE = Path(__file__).resolve().parent / "etc" / "pyspectral.yaml"
def recursive_dict_update(d, u):
@@ -72,15 +56,13 @@ def recursive_dict_update(d, u):
return d
-def get_config():
- """Get the configuration from file."""
- if CONFIG_FILE is not None:
- configfile = CONFIG_FILE
- else:
- configfile = BUILTIN_CONFIG_FILE
+def get_config(config_file: str | Path = None) -> dict:
+ """Get configuration options from YAML file."""
+ if config_file is None:
+ config_file = _get_env_or_builtin_config_path()
config = {}
- with open(configfile, 'r') as fp_:
+ with open(config_file, 'r') as fp_:
config = recursive_dict_update(config, yaml.load(fp_, Loader=UnsafeLoader))
app_dirs = AppDirs('pyspectral', 'pytroll')
@@ -91,3 +73,13 @@ def get_config():
os.makedirs(config['rayleigh_dir'], exist_ok=True)
return config
+
+
+def _get_env_or_builtin_config_path() -> Path:
+ config_file = os.environ.get('PSP_CONFIG_FILE')
+ if config_file is not None and not os.path.isfile(config_file):
+ raise IOError(f"{config_file} pointed to by the environment variable "
+ f"'PSP_CONFIG_FILE' is not a file or does not exist!")
+ if config_file is None:
+ return BUILTIN_CONFIG_FILE
+ return Path(config_file)
=====================================
pyspectral/etc/pyspectral.yaml
=====================================
@@ -50,6 +50,9 @@ download_from_internet: True
# path: /home/a000680/data/SpectralResponses/modis/aqua
# tb2rad_lut_filename: /path/to/radiance/tb/lut/data/tb2rad_lut_aqua_modis_ir3.7.npz
+#DSCOVR-epic:
+# path: C:/Users/Simon/Downloads/EPIC_Filters_Original_Data.xlsx
+
# NOAA-20-viirs:
# # J1_VIIRS_RSR_M1_Detector_Fused_V2.txt
=====================================
pyspectral/near_infrared_reflectance.py
=====================================
@@ -38,7 +38,7 @@ except ImportError:
from pyspectral.config import get_config
from pyspectral.radiance_tb_conversion import RadTbConverter
-from pyspectral.solar import TOTAL_IRRADIANCE_SPECTRUM_2000ASTM, SolarIrradianceSpectrum
+from pyspectral.solar import SolarIrradianceSpectrum
from pyspectral.utils import BANDNAMES, WAVE_LENGTH, get_bandname_from_wavelength
LOG = logging.getLogger(__name__)
@@ -161,8 +161,7 @@ class Calculator(RadTbConverter):
def _get_solarflux(self):
"""Derive the in-band solar flux from rsr over the Near IR band (3.7 or 3.9 microns)."""
solar_spectrum = \
- SolarIrradianceSpectrum(TOTAL_IRRADIANCE_SPECTRUM_2000ASTM,
- dlambda=0.0005,
+ SolarIrradianceSpectrum(dlambda=0.0005,
wavespace=self.wavespace)
self.solar_flux = solar_spectrum.inband_solarflux(self.rsr[self.bandname])
@@ -243,11 +242,11 @@ class Calculator(RadTbConverter):
# 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 +267,8 @@ class Calculator(RadTbConverter):
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 @@ class RadTbConverter(object):
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/solar.py
=====================================
@@ -2,7 +2,7 @@
#
# -*- coding: utf-8 -*-
#
-# Copyright (c) 2013-2022 Pytroll developers
+# Copyright (c) 2013-2023 Pytroll developers
#
#
# This program is free software: you can redistribute it and/or modify
@@ -26,17 +26,16 @@ various instrument bands given their relative spectral response functions
"""
import logging
+from pathlib import Path
import numpy as np
-from pkg_resources import resource_filename
LOG = logging.getLogger(__name__)
# STANDARD SPECTRA from Air Mass Zero: http://rredc.nrel.gov/solar/spectra/am0/
# * 2000 ASTM Standard Extraterrestrial Spectrum Reference E-490-00
# (119.5 - 1,000,000.0 nm)
-TOTAL_IRRADIANCE_SPECTRUM_2000ASTM = resource_filename(__name__,
- 'data/e490_00a.dat')
+TOTAL_IRRADIANCE_SPECTRUM_2000ASTM = Path(__file__).resolve().parent / "data" / "e490_00a.dat"
class SolarIrradianceSpectrum(object):
@@ -48,14 +47,26 @@ class SolarIrradianceSpectrum(object):
in units of W/m^2/micron
"""
- def __init__(self, filename, **options):
+ def __init__(self, filename=TOTAL_IRRADIANCE_SPECTRUM_2000ASTM, **options):
"""Initialize the top of atmosphere solar irradiance spectrum object from file.
+ By default, this will use the following spectra:
+ 2000 ASTM Standard Extraterrestrial Spectrum Reference E-490-00
+
+ To use a different spectra, specify the `filename` when initialising the class.
+
Input:
- filename: Filename of the solar irradiance spectrum
- dlambda:
- Delta wavelength: the step in wavelength defining the resolution on
- which to integrate/convolute.
+ filename: Filename of the solar irradiance spectrum (default: 2000 ASTM)
+ options:
+ dlambda:
+ Delta wavelength: the step in wavelength defining the resolution on
+ which to integrate/convolute.
+ Default is 0.005 if 'wavespace' is 'wavelength' and 2.0 if 'wavenumber'.
+ wavespace:
+ It is possible to specify if the solar irradiance spectrum should
+ be given in terms of wavelength (default) or in terms of
+ wavenumber. If the latter is desired 'wavespace' should be set to
+ 'wavenumber'.
"""
self.wavelength = None
=====================================
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,92 +37,65 @@ WN_RAD_11MICRON_301KELVIN = 0.00117547716523
__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, [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])
+ 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."""
@@ -129,15 +103,10 @@ class TestBlackbody(unittest.TestCase):
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
=====================================
@@ -19,24 +19,17 @@
"""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')
+TEST_TBS = np.array([200., 270., 300., 302., 350.])
TRUE_RADS = np.array([856.937353205, 117420.385297,
479464.582505, 521412.9511, 2928735.18944],
@@ -45,8 +38,7 @@ TRUE_RADS_SEVIRI = np.array([2.391091e-08,
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'] = {}
@@ -231,10 +223,10 @@ class RSRTestDataModis(object):
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
@@ -247,32 +239,38 @@ class TestSeviriConversions(unittest.TestCase):
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
@@ -281,12 +279,14 @@ class TestSeviriConversions(unittest.TestCase):
"""
# 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):
@@ -327,33 +327,106 @@ class TestRadTbConversions(unittest.TestCase):
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/tests/test_rayleigh.py
=====================================
@@ -219,7 +219,7 @@ class TestRayleighDask:
refl_corr = rayl.get_reflectance(sun_zenith, sat_zenith, azidiff, 1.2, redband_refl)
np.testing.assert_allclose(refl_corr, 0)
assert isinstance(refl_corr, da.Array)
- assert rsr_obj.not_called()
+ rsr_obj.assert_not_called()
class TestRayleigh:
@@ -385,7 +385,7 @@ class TestRayleigh:
refl_corr = rayl.get_reflectance(sun_zenith, sat_zenith, azidiff, 0.634, redband_refl)
np.testing.assert_allclose(refl_corr, TEST_RAYLEIGH_RESULT3)
assert isinstance(refl_corr, np.ndarray)
- assert rsr_obj.not_called()
+ rsr_obj.assert_not_called()
@patch('pyspectral.rayleigh.da', None)
def test_get_reflectance_wvl_outside_range(self, fake_lut_hdf5):
@@ -402,7 +402,7 @@ class TestRayleigh:
refl_corr = rayl.get_reflectance(sun_zenith, sat_zenith, azidiff, 1.2, redband_refl)
np.testing.assert_allclose(refl_corr, 0)
assert isinstance(refl_corr, np.ndarray)
- assert rsr_obj.not_called()
+ rsr_obj.assert_not_called()
def test_get_reflectance_no_lut(self, fake_lut_hdf5):
"""Test that missing a LUT causes an exception.."""
=====================================
pyspectral/tests/test_reflectance.py
=====================================
@@ -203,13 +203,25 @@ class TestReflectance(unittest.TestCase):
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 ComputeCountingScheduler
+
+ with dask.config.set(scheduler=ComputeCountingScheduler(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
=====================================
pyspectral/tests/test_solarflux.py
=====================================
@@ -24,7 +24,7 @@ import unittest
import numpy as np
-from pyspectral.solar import TOTAL_IRRADIANCE_SPECTRUM_2000ASTM, SolarIrradianceSpectrum
+from pyspectral.solar import SolarIrradianceSpectrum
TEST_RSR = {}
TEST_RSR['det-1'] = {}
@@ -59,8 +59,7 @@ class TestSolarflux(unittest.TestCase):
def setUp(self):
"""Set up."""
- self.solar_irr = SolarIrradianceSpectrum(TOTAL_IRRADIANCE_SPECTRUM_2000ASTM,
- dlambda=0.005)
+ self.solar_irr = SolarIrradianceSpectrum(dlambda=0.005)
self.rsr = TEST_RSR
def test_read(self):
=====================================
pyspectral/tests/unittest_helpers.py
=====================================
@@ -28,3 +28,21 @@ def assertNumpyArraysEqual(self, other):
raise AssertionError("Shapes don't match")
if not np.allclose(self, other):
raise AssertionError("Elements don't match!")
+
+
+class ComputeCountingScheduler:
+ """Scheduler raising an exception if data are 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)
=====================================
pyspectral/utils.py
=====================================
@@ -1,7 +1,7 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
-# Copyright (c) 2014-2022 Pytroll developers
+# Copyright (c) 2014-2023 Pytroll developers
#
#
# This program is free software: you can redistribute it and/or modify
@@ -23,6 +23,8 @@ import logging
import os
import tarfile
import warnings
+from functools import wraps
+from inspect import getfullargspec
import numpy as np
import requests
@@ -91,17 +93,18 @@ INSTRUMENTS = {'Envisat': 'aatsr',
'NOAA-21': 'viirs',
'Suomi-NPP': 'viirs',
'FY-3B': 'virr',
- 'FY-3C': 'virr'}
+ 'FY-3C': 'virr',
+ 'DSCOVR': 'epic'}
INSTRUMENT_TRANSLATION_DASH2SLASH = {'avhrr-1': 'avhrr/1',
'avhrr-2': 'avhrr/2',
'avhrr-3': 'avhrr/3'}
-HTTP_PYSPECTRAL_RSR = "https://zenodo.org/record/7311143/files/pyspectral_rsr_data.tgz"
+HTTP_PYSPECTRAL_RSR = "https://zenodo.org/records/10029746/files/pyspectral_rsr_data.tgz"
RSR_DATA_VERSION_FILENAME = "PYSPECTRAL_RSR_VERSION"
-RSR_DATA_VERSION = "v1.2.2"
+RSR_DATA_VERSION = "v1.2.4"
ATM_CORRECTION_LUT_VERSION = {}
ATM_CORRECTION_LUT_VERSION['antarctic_aerosol'] = {'version': 'v1.0.1',
@@ -579,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
=====================================
rsr_convert_scripts/epic_reader.py
=====================================
@@ -0,0 +1,120 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Copyright (c) 2023 Pytroll developers
+#
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+"""Read the DSCOVR-EPIC spectral response functions.
+
+Data from the NASA Goddard website:
+https://avdc.gsfc.nasa.gov/pub/DSCOVR/EPIC_Filter_Data/
+
+"""
+
+import logging
+import os
+
+import numpy as np
+import pandas as pd
+
+from pyspectral.raw_reader import InstrumentRSR
+from pyspectral.utils import convert2hdf5 as tohdf5
+
+LOG = logging.getLogger(__name__)
+
+
+EPIC_BAND_NAMES = {'B317': 'T(317.5) %',
+ 'B325': 'T(325) %',
+ 'B340': 'T(340) %',
+ 'B388': 'T(388) %',
+ 'B443': 'T(443) %',
+ 'B551': 'T(551) %',
+ 'B680': 'T(680) %',
+ 'B688': 'T(687.75) %',
+ 'B764': 'T(764) %',
+ 'B780': 'T(779.5) %'}
+
+#: Default time format
+_DEFAULT_TIME_FORMAT = '%Y-%m-%d %H:%M:%S'
+
+#: Default log format
+_DEFAULT_LOG_FORMAT = '[%(levelname)s: %(asctime)s : %(name)s] %(message)s'
+
+
+class EpicRSR(InstrumentRSR):
+ """Container for the DSCOVR EPIC relative spectral response data."""
+
+ def __init__(self, bandname, platform_name):
+ """Initialize the EPIC RSR class."""
+ super(EpicRSR, self).__init__(
+ bandname, platform_name, EPIC_BAND_NAMES.keys())
+
+ self.instrument = 'epic'
+ self._get_options_from_config()
+
+ LOG.debug("Filename: %s", str(self.path))
+ if os.path.exists(self.path):
+ self._load()
+
+ else:
+ LOG.warning("Couldn't find an existing file for this band: %s",
+ str(self.bandname))
+
+ # To be compatible with VIIRS....
+ self.filename = self.requested_band_filename
+
+ self.unit = 'micrometer'
+ self.wavespace = 'wavelength'
+
+ def _load(self, scale=10000.0):
+ """Load the EPIC relative spectral responses."""
+ df1 = pd.read_excel(self.path,
+ sheet_name='Data',
+ skiprows=4,
+ engine='openpyxl')
+ # Remove empty row from the data
+ df1.drop(df1.index[0], inplace=True)
+ # Column names don't match band names - so we use a dict to find correct
+ # columns. We also need to find the column with the wavelength data.
+ # This is the column before the actual RSR data for each band.
+ col_pos = df1.columns.get_loc(EPIC_BAND_NAMES[self.bandname])
+ wvl_data = df1.iloc[:, col_pos - 1]
+ srf_data = df1.iloc[:, col_pos]
+
+ # Not all bands have an identical number of RSR points, so we need to
+ # remove NaNs from the data.
+ wvl_data.dropna(inplace=True)
+ srf_data.dropna(inplace=True)
+
+ # Data is in nanometers, so we need to convert to micrometers.
+ self.rsr = {'wavelength': np.array(wvl_data) / 1000,
+ 'response': np.array(srf_data) / np.nanmax(srf_data)}
+
+
+if __name__ == "__main__":
+ import sys
+ LOG = logging.getLogger('epic_rsr')
+ handler = logging.StreamHandler(sys.stderr)
+
+ formatter = logging.Formatter(fmt=_DEFAULT_LOG_FORMAT,
+ datefmt=_DEFAULT_TIME_FORMAT)
+ handler.setFormatter(formatter)
+ handler.setLevel(logging.DEBUG)
+ LOG.setLevel(logging.DEBUG)
+ LOG.addHandler(handler)
+
+ for platform_name in ['DSCOVR', ]:
+ tohdf5(EpicRSR, platform_name, list(EPIC_BAND_NAMES.keys()))
=====================================
setup.py
=====================================
@@ -82,6 +82,6 @@ setup(name=NAME,
'pyspectral/data/MSG_SEVIRI_Spectral_Response_Characterisation.XLS'])],
test_suite='pyspectral.tests.suite',
tests_require=test_requires,
- python_requires='>=3.7',
+ python_requires='>=3.10',
zip_safe=False,
)
View it on GitLab: https://salsa.debian.org/debian-gis-team/pyspectral/-/commit/14cabd37b0f715e985b2502d04ad528226b6a72e
--
View it on GitLab: https://salsa.debian.org/debian-gis-team/pyspectral/-/commit/14cabd37b0f715e985b2502d04ad528226b6a72e
You're receiving this email because of your account on salsa.debian.org.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/pkg-grass-devel/attachments/20231202/0cda4635/attachment-0001.htm>
More information about the Pkg-grass-devel
mailing list