[Git][debian-gis-team/pyspectral][upstream] New upstream version 0.11.0+ds
Antonio Valentino (@antonio.valentino)
gitlab at salsa.debian.org
Sun Mar 6 08:47:02 GMT 2022
Antonio Valentino pushed to branch upstream at Debian GIS Project / pyspectral
Commits:
98af4161 by Antonio Valentino at 2022-03-06T08:26:03+00:00
New upstream version 0.11.0+ds
- - - - -
11 changed files:
- + AUTHORS.md
- CHANGELOG.md
- pyspectral/config.py
- pyspectral/etc/pyspectral.yaml
- pyspectral/near_infrared_reflectance.py
- pyspectral/rayleigh.py
- pyspectral/tests/data.py
- pyspectral/tests/test_rayleigh.py
- pyspectral/utils.py
- rsr_convert_scripts/README.rst
- + rsr_convert_scripts/msu_gs_reader.py
Changes:
=====================================
AUTHORS.md
=====================================
@@ -0,0 +1,19 @@
+# Project Contributors
+
+The following people have made contributions to this project:
+
+<!--- Use your GitHub account or any other personal reference URL --->
+<!--- If you wish to not use your real name, please use your github username --->
+<!--- The list should be alphabetical by last name if possible, with github usernames at the bottom --->
+<!--- See https://gist.github.com/djhoese/52220272ec73b12eb8f4a29709be110d for auto-generating parts of this list --->
+
+- [K.-Michael Aye](https://github.com/michaelaye)
+- [Adam Dybbroe (adybbroe)](https://github.com/adybbroe)
+- [David Hoese (djhoese)](https://github.com/djhoese)
+- [Christian Kliche](https://github.com/chk)
+- [Panu Lahtinen (pnuu)](https://github.com/pnuu)
+- [Simon R. Proud (simonrp84)](https://github.com/simonrp84)
+- [Martin Raspaud (mraspaud)](https://github.com/mraspaud)
+- [Ronald Scheirer](https://github.com/)
+- [Xin Zhang (zxdawn)](https://github.com/zxdawn)
+- [Rolf Helge ...](https://github.com/HelgeCPH)
\ No newline at end of file
=====================================
CHANGELOG.md
=====================================
@@ -1,3 +1,16 @@
+## Version <v0.11.0> (2022/03/01)
+
+
+### Pull Requests Merged
+
+#### Features added
+
+* [PR 144](https://github.com/pytroll/pyspectral/pull/144) - Add converter for Electro-L N2 MSU/GS spectral response functions
+* [PR 133](https://github.com/pytroll/pyspectral/pull/133) - Refactor rayleigh correction to be more dask friendly
+
+In this release 2 pull requests were closed.
+
+
## Version <v0.10.6> (2021/12/22)
### Issues Closed
=====================================
pyspectral/config.py
=====================================
@@ -87,5 +87,7 @@ def get_config():
user_datadir = app_dirs.user_data_dir
config['rsr_dir'] = expanduser(config.get('rsr_dir', user_datadir))
config['rayleigh_dir'] = expanduser(config.get('rayleigh_dir', user_datadir))
+ os.makedirs(config['rsr_dir'], exist_ok=True)
+ os.makedirs(config['rayleigh_dir'], exist_ok=True)
return config
=====================================
pyspectral/etc/pyspectral.yaml
=====================================
@@ -226,6 +226,19 @@ download_from_internet: True
# ch13: FY4A_AGRI_SRF_CH13.txt
# ch14: FY4A_AGRI_SRF_CH14.txt
+#Electro-L-N2-msu-gs:
+# path: /path/to/rsrs/rtcoef_electro-l_2_msugs_srf/
+# ch1: rtcoef_electro-l_2_msugs_srf_ch01.txt
+# ch2: rtcoef_electro-l_2_msugs_srf_ch02.txt
+# ch3: rtcoef_electro-l_2_msugs_srf_ch03.txt
+# ch4: rtcoef_electro-l_2_msugs_srf_ch04.txt
+# ch5: rtcoef_electro-l_2_msugs_srf_ch05.txt
+# ch6: rtcoef_electro-l_2_msugs_srf_ch06.txt
+# ch7: rtcoef_electro-l_2_msugs_srf_ch07.txt
+# ch8: rtcoef_electro-l_2_msugs_srf_ch08.txt
+# ch9: rtcoef_electro-l_2_msugs_srf_ch09.txt
+# ch10: rtcoef_electro-l_2_msugs_srf_ch10.txt
+
#FY-3B-virr:
# path: /Users/davidh/repos/git/pyspectral/virr_srf/FY3B-VIRR
# ch1: ch1.prn
=====================================
pyspectral/near_infrared_reflectance.py
=====================================
@@ -31,6 +31,7 @@ window channel (usually around 11-12 microns).
import os
import logging
+import tempfile
import numpy as np
try:
from dask.array import where, logical_or, asanyarray, isnan
@@ -40,7 +41,6 @@ except ImportError:
from pyspectral.solar import (SolarIrradianceSpectrum,
TOTAL_IRRADIANCE_SPECTRUM_2000ASTM)
from pyspectral.utils import BANDNAMES, get_bandname_from_wavelength
-from pyspectral.utils import TB2RAD_DIR
from pyspectral.utils import WAVE_LENGTH
from pyspectral.radiance_tb_conversion import RadTbConverter
from pyspectral.config import get_config
@@ -108,6 +108,7 @@ class Calculator(RadTbConverter):
self.lutfile = None
platform_sensor = platform_name + '-' + instrument
+ tb2rad_dir = options.get('tb2rad_dir', tempfile.gettempdir())
if platform_sensor in options and 'tb2rad_lut_filename' in options[platform_sensor]:
if isinstance(options[platform_sensor]['tb2rad_lut_filename'], dict):
for item in options[platform_sensor]['tb2rad_lut_filename']:
@@ -126,14 +127,14 @@ class Calculator(RadTbConverter):
if self.lutfile and not os.path.exists(os.path.dirname(self.lutfile)):
LOG.warning(
"Directory %s does not exist! Check config file", os.path.dirname(self.lutfile))
- self.lutfile = os.path.join(TB2RAD_DIR, os.path.basename(self.lutfile))
+ self.lutfile = os.path.join(tb2rad_dir, os.path.basename(self.lutfile))
if self.lutfile is None:
LOG.info("No lut filename available in config file. "
"Will generate filename automatically")
lutname = 'tb2rad_lut_{0}_{1}_{band}'.format(
self.platform_name.lower(), self.instrument.lower(), band=self.bandname.lower())
- self.lutfile = os.path.join(TB2RAD_DIR, lutname + '.npz')
+ self.lutfile = os.path.join(tb2rad_dir, lutname + '.npz')
LOG.info("lut filename: " + str(self.lutfile))
if not os.path.exists(self.lutfile):
=====================================
pyspectral/rayleigh.py
=====================================
@@ -24,48 +24,52 @@
"""Atmospheric correction of shortwave imager bands in the wavelength range 400 to 800 nm."""
import os
-import time
import logging
+import numbers
import h5py
import numpy as np
try:
- from dask.array import (where, zeros, clip, rad2deg, deg2rad, cos, arccos,
- atleast_2d, Array, map_blocks, from_array, nan_to_num)
import dask.array as da
- HAVE_DASK = True
- # try:
- # # use serializable h5py wrapper to make sure files are closed properly
- # import h5pickle as h5py
- # except ImportError:
- # pass
except ImportError:
- from numpy import where, zeros, clip, rad2deg, deg2rad, cos, arccos, atleast_2d, nan_to_num
da = None
- map_blocks = None
- from_array = None
- Array = None
- HAVE_DASK = False
from geotiepoints.multilinear import MultilinearInterpolator
from pyspectral.rsr_reader import RelativeSpectralResponse
-from pyspectral.utils import (INSTRUMENTS, RAYLEIGH_LUT_DIRS,
+from pyspectral.utils import (INSTRUMENTS, get_rayleigh_lut_dir,
AEROSOL_TYPES, ATMOSPHERES,
BANDNAMES,
ATM_CORRECTION_LUT_VERSION,
- download_luts, get_central_wave,
- get_bandname_from_wavelength)
-
+ download_luts, get_central_wave)
from pyspectral.config import get_config
LOG = logging.getLogger(__name__)
-class BandFrequencyOutOfRange(ValueError):
- """Exception when the band frequency is out of the visible range."""
+def _map_blocks_or_direct_call(func, *args, **kwargs):
+ """Call dask's map_blocks or call func directly if dask is not available."""
+ if da is None:
+ kwargs.pop("meta", None)
+ kwargs.pop("dtype", None)
+ kwargs.pop("chunks", None)
+ return func(*args, **kwargs)
+ return da.map_blocks(func, *args, **kwargs)
+
+
+def _clip_angles_inside_coordinate_range(zenith_angle, zenith_secant_max):
+ """Clipping solar- or satellite zenith angles to be inside the allowed coordinate range.
- pass
+ The "allowed" coordinate range is 0 to the angle corresponiding to a
+ maximum secant. The maximum secant is here the maximum available in the
+ Rayleigh LUT file.
+
+ Also the nan's in the angle array are filled with zeros (0)!
+
+ """
+ clip_angle = np.nan_to_num(np.rad2deg(np.arccos(1. / zenith_secant_max)))
+ zang = np.nan_to_num(zenith_angle)
+ return np.clip(zang, 0, clip_angle)
class RayleighConfigBaseClass(object):
@@ -74,12 +78,9 @@ class RayleighConfigBaseClass(object):
def __init__(self, aerosol_type, atm_type='us-standard'):
"""Initialize class and determine LUT to use."""
options = get_config()
- self.do_download = False
+ self.do_download = 'download_from_internet' in options and options['download_from_internet']
self._lutfiles_version_uptodate = False
- if 'download_from_internet' in options and options['download_from_internet']:
- self.do_download = True
-
if atm_type not in ATMOSPHERES:
raise AttributeError('Atmosphere type not supported! ' +
'Need to be one of {}'.format(str(ATMOSPHERES)))
@@ -99,14 +100,13 @@ class RayleighConfigBaseClass(object):
self.lutfiles_version_uptodate = True
def _get_lutfiles_version(self):
- """
- Get LUT file version.
+ """Get LUT file version.
Check the version of the atm correction luts from the version file in the
specific aerosol correction directory.
"""
- basedir = RAYLEIGH_LUT_DIRS[self._aerosol_type]
+ basedir = get_rayleigh_lut_dir(self._aerosol_type)
lutfiles_version_path = os.path.join(basedir,
ATM_CORRECTION_LUT_VERSION[self._aerosol_type]['filename'])
@@ -157,8 +157,7 @@ class Rayleigh(RayleighConfigBaseClass):
self.sensor = sensor.replace('/', '')
- rayleigh_dir = RAYLEIGH_LUT_DIRS[aerosol_type]
-
+ rayleigh_dir = get_rayleigh_lut_dir(aerosol_type)
ext = atm_type.replace(' ', '_')
lutname = "rayleigh_lut_{0}.h5".format(ext)
self.reflectance_lut_filename = os.path.join(rayleigh_dir, lutname)
@@ -174,216 +173,180 @@ class Rayleigh(RayleighConfigBaseClass):
str(self.reflectance_lut_filename))
LOG.debug('LUT filename: %s', str(self.reflectance_lut_filename))
- self._rayl = None
- self._wvl_coord = None
- self._azid_coord = None
- self._satz_sec_coord = None
- self._sunz_sec_coord = None
-
- def get_effective_wavelength(self, bandname):
- """Get the effective wavelength with Rayleigh scattering in mind."""
- try:
- rsr = RelativeSpectralResponse(self.platform_name, self.sensor)
- except OSError:
- LOG.exception(
- "No spectral responses for this platform and sensor: %s %s", self.platform_name, self.sensor)
- if isinstance(bandname, (float, int)):
- LOG.warning(
- "Effective wavelength is set to the requested band wavelength = %f", bandname)
- return bandname
-
- msg = ("Can't get effective wavelength for band %s on platform %s and sensor %s" %
- (str(bandname), self.platform_name, self.sensor))
- raise KeyError(msg)
-
- if isinstance(bandname, str):
- bandname = BANDNAMES.get(self.sensor, BANDNAMES['generic']).get(bandname, bandname)
- elif isinstance(bandname, (float, int)):
- if not(0.4 < bandname < 0.8):
- raise BandFrequencyOutOfRange(
- 'Requested band frequency should be between 0.4 and 0.8 microns!')
- bandname = get_bandname_from_wavelength(self.sensor, bandname, rsr.rsr)
-
- wvl, resp = (rsr.rsr[bandname]['det-1']['wavelength'],
- rsr.rsr[bandname]['det-1']['response'])
-
- cwvl = get_central_wave(wvl, resp, weight=1. / wvl**4)
- LOG.debug("Band name: %s Effective wavelength: %f", bandname, cwvl)
-
- return cwvl
-
- def get_reflectance_lut(self):
- """Get reflectance LUT.
-
- If not already cached (read previously) read the file with Look-Up
- Tables of reflectances as a function of wavelength, satellite zenith
- secant, azimuth difference angle, and sun zenith secant.
-
- """
- if self._rayl is None:
- lut_vars = get_reflectance_lut_from_file(self.reflectance_lut_filename)
- self._rayl = lut_vars[0]
- self._wvl_coord = lut_vars[1]
- self._azid_coord = lut_vars[2]
- self._satz_sec_coord = lut_vars[3]
- self._sunz_sec_coord = lut_vars[4]
- return self._rayl, self._wvl_coord, self._azid_coord,\
- self._satz_sec_coord, self._sunz_sec_coord
-
- @staticmethod
- def _do_interp(minterp, sunzsec, azidiff, satzsec):
- interp_points2 = np.vstack((sunzsec.ravel(), 180 - azidiff.ravel(), satzsec.ravel()))
- res = minterp(interp_points2)
- return res.reshape(sunzsec.shape)
- def _clip_angles_inside_coordinate_range(self, zenith_angle, zenith_secant_max):
- """Clipping solar- or satellite zenith angles to be inside the allowed coordinate range.
+ def _get_effective_wavelength_and_band_name(self, band_name_or_wavelength):
+ """Get the effective wavelength in nanometers and name of the band/channel.
- The "allowed" coordinate range is 0 to the angle corresponiding to a
- maximum secant. The maximum secant is here the maximum available in the
- Rayleigh LUT file.
+ The provided argument can be either the name of a band/channel that
+ pyspectral knows about in its RSR files or it can be a wavelength
+ specified in micrometers. If a wavelength is provided it is converted
+ to nanometers and returned. If a band name is provided then it is used
+ to lookup the spectral response and determine an effective wavelength
+ useful for rayleigh correction.
- Also the nan's in the angle array are filled with zeros (0)!
+ Returns:
+ A two element tuple of (wavelength in nanometers, name of band). If
+ a wavelength is provided then the band name is a string
+ representation of the wavelength.
"""
- clip_angle = nan_to_num(rad2deg(arccos(1. / zenith_secant_max)))
- zang = nan_to_num(zenith_angle)
- return clip(zang, 0, clip_angle)
-
- def get_reflectance(self, sun_zenith, sat_zenith, azidiff, bandname, redband=None):
- """Get the reflectance from the three sun-sat angles."""
# Get wavelength in nm for band:
- if isinstance(bandname, float):
+ if isinstance(band_name_or_wavelength, numbers.Real):
LOG.warning('A wavelength is provided instead of band name - ' +
'disregard the relative spectral responses and assume ' +
- 'it is the effective wavelength: %f (micro meter)', bandname)
- wvl = bandname * 1000.0
+ 'it is the effective wavelength: %f (micro meter)', band_name_or_wavelength)
+ wvl = band_name_or_wavelength
+ band_name = f'{wvl:f}um'
else:
- wvl = self.get_effective_wavelength(bandname)
- wvl = wvl * 1000.0
-
- rayl, wvl_coord, azid_coord, satz_sec_coord, sunz_sec_coord = self.get_reflectance_lut()
-
- # force dask arrays
- compute = False
- if HAVE_DASK and not isinstance(sun_zenith, Array):
- compute = True
- sun_zenith = from_array(sun_zenith, chunks=sun_zenith.shape)
- sat_zenith = from_array(sat_zenith, chunks=sat_zenith.shape)
- azidiff = from_array(azidiff, chunks=azidiff.shape)
- if redband is not None:
- redband = from_array(redband, chunks=redband.shape)
-
- sun_zenith = self._clip_angles_inside_coordinate_range(sun_zenith, sunz_sec_coord.max())
- sunzsec = 1. / cos(deg2rad(sun_zenith))
+ band_name = band_name_or_wavelength
+ wvl = _get_rsr_wavelength_from_band_name(
+ self.platform_name,
+ self.sensor,
+ band_name,
+ )
+ band_name_map = BANDNAMES.get(self.sensor, BANDNAMES['generic'])
+ band_name_or_wavelength = band_name_map.get(band_name, band_name)
+ LOG.debug("Band name: %s Effective wavelength: %fum", band_name_or_wavelength, wvl)
+ return wvl * 1000.0, band_name
- sat_zenith = self._clip_angles_inside_coordinate_range(sat_zenith, satz_sec_coord.max())
- satzsec = 1. / cos(deg2rad(sat_zenith))
+ @staticmethod
+ def _interp_rayleigh_refl_by_angles(sun_zenith, sat_zenith, azidiff,
+ rayleigh_refl, reflectance_lut_filename):
+ azid_coord, satz_sec_coord, sunz_sec_coord = get_reflectance_lut_from_file(
+ reflectance_lut_filename)
- shape = sun_zenith.shape
+ sun_zenith = _clip_angles_inside_coordinate_range(sun_zenith, sunz_sec_coord.max())
+ sunzsec = 1. / np.cos(np.deg2rad(sun_zenith))
- if not wvl_coord.min() < wvl < wvl_coord.max():
- LOG.warning(
- "Effective wavelength for band %s outside 400-800 nm range!",
- str(bandname))
- LOG.info(
- "Set the rayleigh/aerosol reflectance contribution to zero!")
- if HAVE_DASK:
- chunks = sun_zenith.chunks if redband is None else redband.chunks
- res = zeros(shape, chunks=chunks)
- return res.compute() if compute else res
-
- return zeros(shape)
-
- idx = np.searchsorted(wvl_coord, wvl)
- wvl1 = wvl_coord[idx - 1]
- wvl2 = wvl_coord[idx]
-
- fac = (wvl2 - wvl) / (wvl2 - wvl1)
- raylwvl = fac * rayl[idx - 1, :, :, :] + (1 - fac) * rayl[idx, :, :, :]
- tic = time.time()
+ sat_zenith = _clip_angles_inside_coordinate_range(sat_zenith, satz_sec_coord.max())
+ satzsec = 1. / np.cos(np.deg2rad(sat_zenith))
smin = [sunz_sec_coord[0], azid_coord[0], satz_sec_coord[0]]
smax = [sunz_sec_coord[-1], azid_coord[-1], satz_sec_coord[-1]]
orders = [
len(sunz_sec_coord), len(azid_coord), len(satz_sec_coord)]
- f_3d_grid = atleast_2d(raylwvl.ravel())
+ f_3d_grid = np.atleast_2d(rayleigh_refl.ravel())
- if HAVE_DASK and isinstance(smin[0], Array):
- # compute all of these at the same time before passing to the interpolator
- # otherwise they are computed separately
- smin, smax, orders, f_3d_grid = da.compute(smin, smax, orders, f_3d_grid)
minterp = MultilinearInterpolator(smin, smax, orders)
minterp.set_values(f_3d_grid)
+ interp_points2 = np.vstack((sunzsec.ravel(), 180 - azidiff.ravel(), satzsec.ravel()))
+ res = minterp(interp_points2)
+ res *= 100
+ return res.reshape(sunzsec.shape)
- if HAVE_DASK:
- ipn = map_blocks(self._do_interp, minterp, sunzsec, azidiff,
- satzsec, dtype=raylwvl.dtype, chunks=azidiff.chunks)
- else:
- ipn = self._do_interp(minterp, sunzsec, azidiff, satzsec)
+ def get_reflectance(self, sun_zenith, sat_zenith, azidiff,
+ band_name_or_wavelength, redband=None):
+ """Get the reflectance from the three sun-sat angles."""
+ # if the user gave us non-dask arrays we'll use the dask
+ # version of the algorithm but return numpy arrays back
+ compute = da is not None and not isinstance(sun_zenith, da.Array)
- LOG.debug("Time - Interpolation: {0:f}".format(time.time() - tic))
+ wvl, band_name = self._get_effective_wavelength_and_band_name(band_name_or_wavelength)
+ try:
+ rayleigh_refl = _get_wavelength_adjusted_lut_rayleigh_reflectance(
+ self.reflectance_lut_filename, wvl)
+ except ValueError:
+ LOG.warning("Effective wavelength for band %s outside "
+ "nominal 400-800 nm range!", str(band_name))
+ LOG.info("Setting the rayleigh/aerosol reflectance contribution to zero!")
+ repr_arr = sun_zenith if redband is None else redband
+ zeros_like = np.zeros_like if isinstance(repr_arr, np.ndarray) else da.zeros_like
+ res = zeros_like(repr_arr)
+ else:
+ res = _map_blocks_or_direct_call(self._interp_rayleigh_refl_by_angles,
+ sun_zenith, sat_zenith, azidiff, rayleigh_refl,
+ self.reflectance_lut_filename,
+ meta=np.array((), dtype=rayleigh_refl.dtype),
+ dtype=rayleigh_refl.dtype,
+ chunks=getattr(azidiff, "chunks", None))
- ipn *= 100
- res = ipn
if redband is not None:
- res = where(redband < 20., res,
- (1 - (redband - 20) / 80) * res)
+ res = _map_blocks_or_direct_call(self._relax_rayleigh_refl_correction_where_cloudy,
+ redband, res,
+ meta=np.array((), dtype=res.dtype),
+ dtype=res.dtype,
+ chunks=getattr(res, "chunks", None))
- res = clip(nan_to_num(res), 0, 100)
+ res = np.clip(np.nan_to_num(res), 0, 100)
if compute:
res = res.compute()
-
return res
+ @staticmethod
+ def _relax_rayleigh_refl_correction_where_cloudy(redband, rayleigh_refl):
+ return np.where(redband < 20., rayleigh_refl,
+ (1 - (redband - 20) / 80) * rayleigh_refl)
+
@staticmethod
def reduce_rayleigh_highzenith(zenith, rayref, thresh_zen, maxzen, strength):
"""Reduce the Rayleigh correction amount at high zenith angles.
+
This linearly scales the Rayleigh reflectance, `rayref`, for solar or satellite zenith angles, `zenith`,
above a threshold angle, `thresh_zen`. Between `thresh_zen` and `maxzen` the Rayleigh reflectance will
be linearly scaled, from one at `thresh_zen` to zero at `maxzen`.
"""
LOG.info("Reducing Rayleigh effect at high zenith angles.")
- factor = 1. - strength * where(zenith < thresh_zen, 0, (zenith - thresh_zen) / (maxzen - thresh_zen))
+ factor = 1. - strength * np.where(zenith < thresh_zen, 0, (zenith - thresh_zen) / (maxzen - thresh_zen))
# For low zenith factor can be greater than one, so we need to clip it into a sensible range.
- factor = clip(factor, 0, 1)
+ factor = np.clip(factor, 0, 1)
return rayref * factor
-def get_reflectance_lut_from_file(filename):
+def _get_rsr_wavelength_from_band_name(platform_name, sensor, band_name):
+ try:
+ rsr = RelativeSpectralResponse(platform_name, sensor)
+ except OSError:
+ LOG.exception(
+ "No spectral responses for this platform and sensor: %s %s", platform_name, sensor)
+ msg = ("Can't get effective wavelength for band %s on platform %s and sensor %s" %
+ (str(band_name), platform_name, sensor))
+ raise KeyError(msg)
+
+ wvl = rsr.rsr[band_name]['det-1']['wavelength']
+ resp = rsr.rsr[band_name]['det-1']['response']
+ cwvl = get_central_wave(wvl, resp, weight=1. / wvl**4)
+ return cwvl
+
+
+def _get_wavelength_adjusted_lut_rayleigh_reflectance(lut_filename, wvl):
+ with h5py.File(lut_filename, 'r') as h5f:
+ rayleigh_refl = h5f["reflectance"]
+ wvl_coord = h5f["wavelengths"][:]
+ if not wvl_coord.min() < wvl < wvl_coord.max():
+ raise ValueError("Wavelength out of range for available LUT wavelengths")
+ wavelength_index, wavelength_factor = _get_wavelength_index_and_factor(wvl_coord, wvl)
+ raylwvl = (wavelength_factor * rayleigh_refl[wavelength_index - 1, :, :, :] +
+ (1 - wavelength_factor) * rayleigh_refl[wavelength_index, :, :, :])
+ return raylwvl
+
+
+def _get_wavelength_index_and_factor(wvl_coord, wvl):
+ wavelength_index = np.searchsorted(wvl_coord, wvl)
+ wvl1 = wvl_coord[wavelength_index - 1]
+ wvl2 = wvl_coord[wavelength_index]
+ wavelength_factor = (wvl2 - wvl) / (wvl2 - wvl1)
+ return wavelength_index, wavelength_factor
+
+
+def get_reflectance_lut_from_file(lut_filename):
"""Get reflectance LUT.
Read the Look-Up Tables from file with reflectances as a function of
wavelength, satellite zenith secant, azimuth difference angle, and sun
- zenith secant
+ zenith secant.
"""
- h5f = h5py.File(filename, 'r')
-
- tab = h5f['reflectance']
- wvl = h5f['wavelengths']
- azidiff = h5f['azimuth_difference']
- satellite_zenith_secant = h5f['satellite_zenith_secant']
- sun_zenith_secant = h5f['sun_zenith_secant']
-
- if HAVE_DASK:
- tab = from_array(tab, chunks=(10, 10, 10, 10))
- wvl = wvl[:] # no benefit to dask-ifying this
- azidiff = from_array(azidiff, chunks=(1000,))
- satellite_zenith_secant = from_array(satellite_zenith_secant,
- chunks=(1000,))
- sun_zenith_secant = from_array(sun_zenith_secant,
- chunks=(1000,))
- else:
+ with h5py.File(lut_filename, 'r') as h5f:
+ azidiff = h5f['azimuth_difference']
+ satellite_zenith_secant = h5f['satellite_zenith_secant']
+ sun_zenith_secant = h5f['sun_zenith_secant']
+
# load all of the data we are going to use in to memory
- tab = tab[:]
- wvl = wvl[:]
azidiff = azidiff[:]
satellite_zenith_secant = satellite_zenith_secant[:]
sun_zenith_secant = sun_zenith_secant[:]
- h5f.close()
- return tab, wvl, azidiff, satellite_zenith_secant, sun_zenith_secant
+ return azidiff, satellite_zenith_secant, sun_zenith_secant
def check_and_download(**kwargs):
=====================================
pyspectral/tests/data.py
=====================================
@@ -268,4 +268,4 @@ TEST_RAYLEIGH_LUT = np.array([[[[0.08674087, 0.09136106, 0.09747749, 0.104837
TEST_RAYLEIGH_AZID_COORD = np.array([100., 110., 120., 130., 140., 150., 160.], dtype='float32')
TEST_RAYLEIGH_SATZ_COORD = np.array([1., 1.25, 1.5, 1.75, 2., 2.25], dtype='float32')
TEST_RAYLEIGH_SUNZ_COORD = np.array([1., 1.25, 1.5, 1.75, 2., 2.25, 2.5, 2.75], dtype='float32')
-TEST_RAYLEIGH_WVL_COORD = np.array([440., 445.], dtype='float32')
+TEST_RAYLEIGH_WVL_COORD = np.array([631., 636.], dtype='float32')
=====================================
pyspectral/tests/test_rayleigh.py
=====================================
@@ -21,37 +21,32 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Unittest for the rayleigh correction utilities."""
-
+import contextlib
import os
-import sys
import numpy as np
import dask.array as da
+import h5py
+import pytest
from pyspectral import rayleigh
-from pyspectral.rayleigh import BandFrequencyOutOfRange
from pyspectral.tests.data import (
TEST_RAYLEIGH_LUT,
TEST_RAYLEIGH_AZID_COORD,
TEST_RAYLEIGH_SUNZ_COORD,
TEST_RAYLEIGH_SATZ_COORD,
TEST_RAYLEIGH_WVL_COORD)
-from pyspectral.utils import RSR_DATA_VERSION
+from pyspectral.utils import ATM_CORRECTION_LUT_VERSION
-import unittest
from unittest.mock import patch
-TEST_RAYLEIGH_RESULT1 = np.array([10.40727436, 8.69775471], dtype='float32')
-TEST_RAYLEIGH_RESULT2 = np.array([9.71696059, 8.51415689], dtype='float32')
-TEST_RAYLEIGH_RESULT3 = np.array([5.61532271, 8.69267476], dtype='float32')
-TEST_RAYLEIGH_RESULT4 = np.array([0.0, 8.69775471], dtype='float32')
+TEST_RAYLEIGH_RESULT1 = np.array([10.339923, 8.64748], dtype='float32')
+TEST_RAYLEIGH_RESULT2 = np.array([9.653559, 8.464997], dtype='float32')
+TEST_RAYLEIGH_RESULT3 = np.array([5.488735, 8.533125], dtype='float32')
+TEST_RAYLEIGH_RESULT4 = np.array([0.0, 8.64748], dtype='float32')
TEST_RAYLEIGH_RESULT_R1 = np.array([16.66666667, 20.83333333, 25.], dtype='float32')
TEST_RAYLEIGH_RESULT_R2 = np.array([0., 6.25, 12.5], dtype='float32')
TEST_ZENITH_ANGLES_RESULTS = np.array([68.67631374, 68.67631374, 32., 0.])
-# Mock some modules, so we don't need them for tests.
-
-# sys.modules['pyresample'] = MagicMock()
-
class RelativeSpectralResponseTestData(object):
"""Create the class instance to hold the RSR test data."""
@@ -87,288 +82,214 @@ class RelativeSpectralResponseTestData(object):
self.rsr[chname]['det-1']['response'] = ch3_resp
-class TestRayleighDask(unittest.TestCase):
+ at pytest.fixture
+def fake_lut_hdf5(tmp_path):
+ """Create a fake LUT HDF5 file and necessary mocks to use it."""
+ aerosol_type = "marine_clean_aerosol"
+ base_dir = tmp_path / aerosol_type
+ base_dir.mkdir()
+ _create_fake_lut_version_file(base_dir, aerosol_type)
+ for atm_type in ("midlatitude_summer", "subarctic_winter", "tropical", "us-standard"):
+ _create_fake_lut_hdf5_file(base_dir, atm_type)
+ fake_config = {
+ "rsr_dir": str(tmp_path),
+ "rayleigh_dir": str(tmp_path),
+ "download_from_internet": False,
+ }
+ with patch('pyspectral.rayleigh.get_config', lambda: fake_config), \
+ patch('pyspectral.utils.get_config', lambda: fake_config):
+ yield None
+
+
+def _create_fake_lut_version_file(base_dir, aerosol_type):
+ lut_version = ATM_CORRECTION_LUT_VERSION[aerosol_type]["version"]
+ lutfiles_path = str(base_dir / ATM_CORRECTION_LUT_VERSION[aerosol_type]["filename"])
+ with open(lutfiles_path, "w") as version_file:
+ version_file.write(lut_version)
+
+
+def _create_fake_lut_hdf5_file(base_dir, atm_type) -> str:
+ filename = str(base_dir / f"rayleigh_lut_{atm_type}.h5")
+ with h5py.File(filename, "w") as h5f:
+ h5f["reflectance"] = TEST_RAYLEIGH_LUT
+ h5f["wavelengths"] = TEST_RAYLEIGH_WVL_COORD
+ h5f["azimuth_difference"] = TEST_RAYLEIGH_AZID_COORD
+ h5f["sun_zenith_secant"] = TEST_RAYLEIGH_SUNZ_COORD
+ h5f["satellite_zenith_secant"] = TEST_RAYLEIGH_SATZ_COORD
+ return filename
+
+
+def _create_rayleigh():
+ rayl = rayleigh.Rayleigh('NOAA-20', 'viirs', atmosphere='midlatitude summer')
+ return rayl
+
+
+ at contextlib.contextmanager
+def mocked_rsr():
+ """Mock the RSR class used by the Rayleigh class with fake data."""
+ with patch('pyspectral.rayleigh.RelativeSpectralResponse') as mymock:
+ instance = mymock.return_value
+ instance.rsr = RelativeSpectralResponseTestData().rsr
+ instance.unit = '1e-6 m'
+ instance.si_scale = 1e-6
+ yield mymock
+
+
+class TestRayleighDask:
"""Class for testing pyspectral.rayleigh - with dask-arrays as input."""
- def setUp(self):
- """Set up the test."""
- self.cwvl = 0.4440124
- self.rsr = RelativeSpectralResponseTestData()
- self._res1 = da.from_array(TEST_RAYLEIGH_RESULT1, chunks=2)
- self._res2 = da.from_array(TEST_RAYLEIGH_RESULT2, chunks=2)
- self._res3 = da.from_array(TEST_RAYLEIGH_RESULT3, chunks=2)
-
- self.rayl = da.from_array(TEST_RAYLEIGH_LUT, chunks=(10, 10, 10, 10))
- self.wvl_coord = TEST_RAYLEIGH_WVL_COORD
- self.azid_coord = da.from_array(TEST_RAYLEIGH_AZID_COORD, chunks=(1000,))
- self.sunz_sec_coord = da.from_array(TEST_RAYLEIGH_SUNZ_COORD,
- chunks=(1000,))
- self.satz_sec_coord = da.from_array(TEST_RAYLEIGH_SATZ_COORD,
- chunks=(1000,))
-
- # mymock:
- with patch('pyspectral.rayleigh.RelativeSpectralResponse') as mymock:
- instance = mymock.return_value
- instance.rsr = RelativeSpectralResponseTestData().rsr
- instance.unit = '1e-6 m'
- instance.si_scale = 1e-6
-
- self.viirs_rayleigh = rayleigh.Rayleigh('NOAA-20', 'viirs', atmosphere='midlatitude summer')
-
- @patch('os.path.exists')
- @patch('pyspectral.utils.download_luts')
- @patch('pyspectral.rayleigh.get_reflectance_lut_from_file')
- @patch('pyspectral.rsr_reader.RelativeSpectralResponse.'
- '_get_rsr_data_version')
- @patch('pyspectral.rayleigh.Rayleigh.get_effective_wavelength')
- def test_dask_cliping_angles_with_nans(self, get_effective_wvl,
- get_rsr_version, get_reflectance_lut_from_file,
- download_luts, exists):
- """Test the cliping of angles outside coordinate range - with nan's in input."""
- get_reflectance_lut_from_file.return_value = (self.rayl, self.wvl_coord, self.azid_coord,
- self.satz_sec_coord, self.sunz_sec_coord)
- download_luts.return_code = None
- exists.return_code = True
- get_rsr_version.return_code = RSR_DATA_VERSION
- get_effective_wvl.return_value = self.cwvl
-
+ def test_dask_clipping_angles_with_nans(self):
+ """Test the clipping of angles outside coordinate range - with nan's in input."""
+ from pyspectral.rayleigh import _clip_angles_inside_coordinate_range
zenith_angle = da.array([79., 69., 32., np.nan])
- result = self.viirs_rayleigh._clip_angles_inside_coordinate_range(zenith_angle, 2.75)
-
+ result = _clip_angles_inside_coordinate_range(zenith_angle, 2.75)
np.testing.assert_allclose(result, TEST_ZENITH_ANGLES_RESULTS)
- @patch('os.path.exists')
- @patch('pyspectral.utils.download_luts')
- @patch('pyspectral.rayleigh.get_reflectance_lut_from_file')
- @patch('pyspectral.rsr_reader.RelativeSpectralResponse.'
- '_get_rsr_data_version')
- @patch('pyspectral.rayleigh.Rayleigh.get_effective_wavelength')
- def test_get_reflectance_dask_redband_outside_clip(self, get_effective_wvl,
- get_rsr_version, get_reflectance_lut_from_file,
- download_luts, exists):
+ def test_get_reflectance_dask_redband_outside_clip(self, fake_lut_hdf5):
"""Test getting the reflectance correction with dask inputs - using red band reflections outside 20-100."""
- get_reflectance_lut_from_file.return_value = (self.rayl, self.wvl_coord, self.azid_coord,
- self.satz_sec_coord, self.sunz_sec_coord)
- download_luts.return_code = None
- exists.return_code = True
- get_rsr_version.return_code = RSR_DATA_VERSION
- get_effective_wvl.return_value = self.cwvl
-
sun_zenith = da.array([67., 32.])
sat_zenith = da.array([45., 18.])
azidiff = da.array([150., 110.])
redband_refl = da.array([108., -0.5])
- refl_corr = self.viirs_rayleigh.get_reflectance(sun_zenith, sat_zenith, azidiff, 'M2', redband_refl)
+ viirs_rayl = _create_rayleigh()
+ with mocked_rsr():
+ refl_corr = viirs_rayl.get_reflectance(sun_zenith, sat_zenith, azidiff, 'ch3', redband_refl)
np.testing.assert_allclose(refl_corr, TEST_RAYLEIGH_RESULT4)
- self.assertIsInstance(refl_corr, da.Array)
-
- @patch('os.path.exists')
- @patch('pyspectral.utils.download_luts')
- @patch('pyspectral.rayleigh.get_reflectance_lut_from_file')
- @patch('pyspectral.rsr_reader.RelativeSpectralResponse.'
- '_get_rsr_data_version')
- @patch('pyspectral.rayleigh.Rayleigh.get_effective_wavelength')
- def test_get_reflectance_dask(self, get_effective_wvl,
- get_rsr_version, get_reflectance_lut_from_file,
- download_luts, exists):
- """Test getting the reflectance correction with dask inputs."""
- get_reflectance_lut_from_file.return_value = (self.rayl, self.wvl_coord, self.azid_coord,
- self.satz_sec_coord, self.sunz_sec_coord)
- download_luts.return_code = None
- exists.return_code = True
- get_rsr_version.return_code = RSR_DATA_VERSION
- get_effective_wvl.return_value = self.cwvl
+ assert isinstance(refl_corr, da.Array)
+ def test_get_reflectance_dask(self, fake_lut_hdf5):
+ """Test getting the reflectance correction with dask inputs."""
sun_zenith = da.array([67., 32.])
sat_zenith = da.array([45., 18.])
azidiff = da.array([150., 110.])
redband_refl = da.array([14., 5.])
- refl_corr = self.viirs_rayleigh.get_reflectance(sun_zenith, sat_zenith, azidiff, 'M2', redband_refl)
+ rayl = _create_rayleigh()
+ with mocked_rsr():
+ refl_corr = rayl.get_reflectance(sun_zenith, sat_zenith, azidiff, 'ch3', redband_refl)
np.testing.assert_allclose(refl_corr, TEST_RAYLEIGH_RESULT1)
- self.assertIsInstance(refl_corr, da.Array)
+ assert isinstance(refl_corr, da.Array)
sun_zenith = da.array([60., 20.])
sat_zenith = da.array([49., 26.])
azidiff = da.array([140., 130.])
redband_refl = da.array([12., 8.])
- refl_corr = self.viirs_rayleigh.get_reflectance(sun_zenith, sat_zenith, azidiff, 'M2', redband_refl)
+ with mocked_rsr():
+ refl_corr = rayl.get_reflectance(sun_zenith, sat_zenith, azidiff, 'ch3', redband_refl)
np.testing.assert_allclose(refl_corr, TEST_RAYLEIGH_RESULT2)
- self.assertIsInstance(refl_corr, da.Array)
-
- @patch('os.path.exists')
- @patch('pyspectral.utils.download_luts')
- @patch('pyspectral.rayleigh.get_reflectance_lut_from_file')
- @patch('pyspectral.rsr_reader.RelativeSpectralResponse.'
- '_get_rsr_data_version')
- @patch('pyspectral.rayleigh.Rayleigh.get_effective_wavelength')
- def test_get_reflectance_numpy_dask(self, get_effective_wvl,
- get_rsr_version, get_reflectance_lut_from_file,
- download_luts, exists):
- """Test getting the reflectance correction with dask inputs."""
- get_reflectance_lut_from_file.return_value = (self.rayl, self.wvl_coord, self.azid_coord,
- self.satz_sec_coord, self.sunz_sec_coord)
- download_luts.return_code = None
- exists.return_code = True
- get_rsr_version.return_code = RSR_DATA_VERSION
- get_effective_wvl.return_value = self.cwvl
+ assert isinstance(refl_corr, da.Array)
+ def test_get_reflectance_numpy_dask(self, fake_lut_hdf5):
+ """Test getting the reflectance correction with dask inputs."""
sun_zenith = np.array([67., 32.])
sat_zenith = np.array([45., 18.])
azidiff = np.array([150., 110.])
redband_refl = np.array([14., 5.])
- refl_corr = self.viirs_rayleigh.get_reflectance(sun_zenith, sat_zenith, azidiff, 'M2', redband_refl)
+ rayl = _create_rayleigh()
+ with mocked_rsr():
+ refl_corr = rayl.get_reflectance(sun_zenith, sat_zenith, azidiff, 'ch3', redband_refl)
np.testing.assert_allclose(refl_corr, TEST_RAYLEIGH_RESULT1)
- self.assertIsInstance(refl_corr, np.ndarray)
+ assert isinstance(refl_corr, np.ndarray)
sun_zenith = np.array([60., 20.])
sat_zenith = np.array([49., 26.])
azidiff = np.array([140., 130.])
redband_refl = np.array([12., 8.])
- refl_corr = self.viirs_rayleigh.get_reflectance(sun_zenith, sat_zenith, azidiff, 'M2', redband_refl)
+ with mocked_rsr():
+ refl_corr = rayl.get_reflectance(sun_zenith, sat_zenith, azidiff, 'ch3', redband_refl)
np.testing.assert_allclose(refl_corr, TEST_RAYLEIGH_RESULT2)
- self.assertIsInstance(refl_corr, np.ndarray)
-
-
-class TestRayleigh(unittest.TestCase):
+ assert isinstance(refl_corr, np.ndarray)
+
+ def test_get_reflectance_wvl_outside_range(self, fake_lut_hdf5):
+ """Test getting the reflectance correction with wavelength outside correction range."""
+ with mocked_rsr() as rsr_obj:
+ rsr_obj.side_effect = IOError("No rsr data in pyspectral for this platform and sensor")
+ sun_zenith = da.from_array([50., 10.])
+ sat_zenith = da.from_array([39., 16.])
+ azidiff = da.from_array([170., 110.])
+ redband_refl = da.from_array([29., 12.])
+ rayl = rayleigh.Rayleigh('UFO', 'unknown', atmosphere='midlatitude summer')
+
+ # we gave a direct wavelength so the RSR object shouldn't be needed
+ 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()
+
+
+class TestRayleigh:
"""Class for testing pyspectral.rayleigh."""
- def setUp(self):
- """Set up the test."""
- self.cwvl = 0.4440124
- self.rsr = RelativeSpectralResponseTestData()
- self._res1 = da.from_array(TEST_RAYLEIGH_RESULT1, chunks=2)
- self._res2 = da.from_array(TEST_RAYLEIGH_RESULT2, chunks=2)
- self._res3 = da.from_array(TEST_RAYLEIGH_RESULT3, chunks=2)
-
- self.rayl = TEST_RAYLEIGH_LUT
- self.wvl_coord = TEST_RAYLEIGH_WVL_COORD
- self.azid_coord = TEST_RAYLEIGH_AZID_COORD
- self.sunz_sec_coord = TEST_RAYLEIGH_SUNZ_COORD
- self.satz_sec_coord = TEST_RAYLEIGH_SATZ_COORD
-
- # mymock:
- with patch('pyspectral.rayleigh.RelativeSpectralResponse') as mymock:
- instance = mymock.return_value
- instance.rsr = RelativeSpectralResponseTestData().rsr
- instance.unit = '1e-6 m'
- instance.si_scale = 1e-6
-
- self.viirs_rayleigh = rayleigh.Rayleigh('NOAA-20', 'viirs', atmosphere='midlatitude summer')
-
- def test_get_effective_wavelength(self):
+ def test_get_effective_wavelength_and_band_name_with_floats(self, fake_lut_hdf5):
"""Test getting the effective wavelength."""
- # mymock:
- with patch('pyspectral.rayleigh.RelativeSpectralResponse') as mymock:
- instance = mymock.return_value
- instance.rsr = RelativeSpectralResponseTestData().rsr
-
- this = rayleigh.Rayleigh('Himawari-8', 'ahi')
- with self.assertRaises(BandFrequencyOutOfRange):
- this.get_effective_wavelength(0.9)
-
- # Only ch3 (~0.63) testdata implemented yet...
- ewl = this.get_effective_wavelength(0.65)
- self.assertAlmostEqual(ewl, 0.6356167)
-
- # mymock:
- with patch('pyspectral.rayleigh.RelativeSpectralResponse') as mymock:
- instance = mymock.side_effect = IOError(
- 'Fake that there is no spectral response file...')
-
- this = rayleigh.Rayleigh('Himawari-8', 'ahi')
- ewl = this.get_effective_wavelength(0.7)
- self.assertEqual(ewl, 0.7)
- ewl = this.get_effective_wavelength(0.9)
- self.assertEqual(ewl, 0.9)
- ewl = this.get_effective_wavelength(0.455)
- self.assertEqual(ewl, 0.455)
-
- @patch('os.path.exists')
- @patch('pyspectral.utils.download_luts')
- def test_rayleigh_init(self, download_luts, exists):
+ this = rayleigh.Rayleigh('Himawari-8', 'ahi')
+ # Only ch3 (~0.63) testdata implemented yet...
+ ewl, band_name = this._get_effective_wavelength_and_band_name(0.65)
+ np.testing.assert_allclose(ewl, 650) # 635.6167)
+ assert isinstance(band_name, str)
+
+ this = rayleigh.Rayleigh('Himawari-8', 'ahi')
+ ewl, band_name = this._get_effective_wavelength_and_band_name(0.7)
+ assert ewl == 700.0
+ assert isinstance(band_name, str)
+ ewl, band_name = this._get_effective_wavelength_and_band_name(0.9)
+ assert ewl == 900.0
+ assert isinstance(band_name, str)
+ ewl, band_name = this._get_effective_wavelength_and_band_name(0.455)
+ assert ewl == 455.0
+ assert isinstance(band_name, str)
+
+ def test_rayleigh_init(self, fake_lut_hdf5):
"""Test creating the Rayleigh object."""
- download_luts.return_code = None
- exists.return_code = True
-
- # mymock:
with patch('pyspectral.rayleigh.RelativeSpectralResponse') as mymock:
instance = mymock.return_value
instance.rsr = RelativeSpectralResponseTestData().rsr
- with self.assertRaises(AttributeError):
- this = rayleigh.Rayleigh('Himawari-8', 'ahi', atmosphere='unknown')
+ with pytest.raises(AttributeError):
+ rayleigh.Rayleigh('Himawari-8', 'ahi', atmosphere='unknown')
- with self.assertRaises(AttributeError):
- this = rayleigh.Rayleigh('Himawari-8', 'ahi', aerosol_type='unknown')
+ with pytest.raises(AttributeError):
+ rayleigh.Rayleigh('Himawari-8', 'ahi', aerosol_type='unknown')
this = rayleigh.Rayleigh('Himawari-8', 'ahi', atmosphere='subarctic winter')
- self.assertTrue(os.path.basename(this.reflectance_lut_filename).endswith('subarctic_winter.h5'))
- self.assertTrue(this.sensor == 'ahi')
+ assert os.path.basename(this.reflectance_lut_filename).endswith('subarctic_winter.h5')
+ assert this.sensor == 'ahi'
this = rayleigh.Rayleigh('NOAA-19', 'avhrr/3', atmosphere='tropical')
- self.assertTrue(this.sensor == 'avhrr3')
-
- @patch('os.path.exists')
- @patch('pyspectral.utils.download_luts')
- @patch('pyspectral.rayleigh.get_reflectance_lut_from_file')
- @patch('pyspectral.rsr_reader.RelativeSpectralResponse.'
- '_get_rsr_data_version')
- @patch('pyspectral.rayleigh.Rayleigh.get_effective_wavelength')
- def test_cliping_angles_with_nans(self, get_effective_wvl,
- get_rsr_version, get_reflectance_lut_from_file,
- download_luts, exists):
- """Test the cliping of angles outside coordinate range - with nan's in input."""
- get_reflectance_lut_from_file.return_value = (self.rayl, self.wvl_coord, self.azid_coord,
- self.satz_sec_coord, self.sunz_sec_coord)
- download_luts.return_code = None
- exists.return_code = True
- get_rsr_version.return_code = RSR_DATA_VERSION
- get_effective_wvl.return_value = self.cwvl
+ assert this.sensor == 'avhrr3'
+ def test_clipping_angles_with_nans(self):
+ """Test the clipping of angles outside coordinate range - with nan's in input."""
+ from pyspectral.rayleigh import _clip_angles_inside_coordinate_range
zenith_angle = np.array([79., 69., 32., np.nan])
- result = self.viirs_rayleigh._clip_angles_inside_coordinate_range(zenith_angle, 2.75)
-
+ result = _clip_angles_inside_coordinate_range(zenith_angle, 2.75)
np.testing.assert_allclose(result, TEST_ZENITH_ANGLES_RESULTS)
- def test_rayleigh_reduction(self):
+ def test_rayleigh_reduction(self, fake_lut_hdf5):
"""Test the code that reduces Rayleigh correction for high zenith angles."""
-
# Test the Rayleigh reduction code
sun_zenith = np.array([70., 65., 60.])
in_rayleigh = [50, 50, 50]
# Test case where no reduction is done.
- retv = self.viirs_rayleigh.reduce_rayleigh_highzenith(sun_zenith, in_rayleigh, 70., 90., 1)
- self.assertTrue(np.allclose(retv, in_rayleigh))
+ retv = rayleigh.Rayleigh.reduce_rayleigh_highzenith(sun_zenith, in_rayleigh, 70., 90., 1)
+ np.testing.assert_allclose(retv, in_rayleigh)
# Test case where moderate reduction is performed.
- retv = self.viirs_rayleigh.reduce_rayleigh_highzenith(sun_zenith, in_rayleigh, 30., 90., 1)
- self.assertTrue(np.allclose(retv, TEST_RAYLEIGH_RESULT_R1))
+ retv = rayleigh.Rayleigh.reduce_rayleigh_highzenith(sun_zenith, in_rayleigh, 30., 90., 1)
+ np.testing.assert_allclose(retv, TEST_RAYLEIGH_RESULT_R1)
# Test case where extreme reduction is performed.
- retv = self.viirs_rayleigh.reduce_rayleigh_highzenith(sun_zenith, in_rayleigh, 30., 90., 1.5)
- self.assertTrue(np.allclose(retv, TEST_RAYLEIGH_RESULT_R2))
-
- @patch('pyspectral.rayleigh.HAVE_DASK', False)
- @patch('os.path.exists')
- @patch('pyspectral.utils.download_luts')
- @patch('pyspectral.rayleigh.get_reflectance_lut_from_file')
- @patch('pyspectral.rsr_reader.RelativeSpectralResponse._get_rsr_data_version')
- @patch('pyspectral.rayleigh.Rayleigh.get_effective_wavelength')
- def test_get_reflectance_redband_outside_clip(self, get_effective_wvl,
- get_rsr_version, get_reflectance_lut_from_file,
- download_luts, exists):
- """Test getting the reflectance correction - using red band reflections outside 20 to 100."""
- get_reflectance_lut_from_file.return_value = (self.rayl, self.wvl_coord, self.azid_coord,
- self.satz_sec_coord, self.sunz_sec_coord)
- download_luts.return_code = None
- exists.return_code = True
- get_rsr_version.return_code = RSR_DATA_VERSION
- get_effective_wvl.return_value = self.cwvl
+ retv = rayleigh.Rayleigh.reduce_rayleigh_highzenith(sun_zenith, in_rayleigh, 30., 90., 1.5)
+ np.testing.assert_allclose(retv, TEST_RAYLEIGH_RESULT_R2)
+ @patch('pyspectral.rayleigh.da', None)
+ def test_get_reflectance_redband_outside_clip(self, fake_lut_hdf5):
+ """Test getting the reflectance correction - using red band reflections outside 20 to 100."""
sun_zenith = np.array([67., 32.])
sat_zenith = np.array([45., 18.])
azidiff = np.array([150., 110.])
redband_refl = np.array([100., 20.])
- refl_corr1 = self.viirs_rayleigh.get_reflectance(
- sun_zenith, sat_zenith, azidiff, 'M2', redband_refl)
+ rayl = _create_rayleigh()
+ with mocked_rsr():
+ refl_corr1 = rayl.get_reflectance(
+ sun_zenith, sat_zenith, azidiff, 'ch3', redband_refl)
np.testing.assert_allclose(refl_corr1, TEST_RAYLEIGH_RESULT4)
@@ -377,71 +298,94 @@ class TestRayleigh(unittest.TestCase):
rints_high = rng.integers(low=100, high=200, size=2).astype('float')
redband_refl = np.array([rints_high[0], rints_low[0]])
- refl_corr2 = self.viirs_rayleigh.get_reflectance(
- sun_zenith, sat_zenith, azidiff, 'M2', redband_refl)
+ with mocked_rsr():
+ refl_corr2 = rayl.get_reflectance(
+ sun_zenith, sat_zenith, azidiff, 'ch3', redband_refl)
redband_refl = np.array([rints_high[1], rints_low[1]])
- refl_corr3 = self.viirs_rayleigh.get_reflectance(
- sun_zenith, sat_zenith, azidiff, 'M2', redband_refl)
+ with mocked_rsr():
+ refl_corr3 = rayl.get_reflectance(
+ sun_zenith, sat_zenith, azidiff, 'ch3', redband_refl)
+ assert isinstance(refl_corr1, np.ndarray)
+ assert isinstance(refl_corr2, np.ndarray)
np.testing.assert_allclose(refl_corr1, refl_corr2)
np.testing.assert_allclose(refl_corr2, refl_corr3)
- @patch('pyspectral.rayleigh.HAVE_DASK', False)
- @patch('os.path.exists')
- @patch('pyspectral.utils.download_luts')
- @patch('pyspectral.rayleigh.get_reflectance_lut_from_file')
- @patch('pyspectral.rsr_reader.RelativeSpectralResponse._get_rsr_data_version')
- @patch('pyspectral.rayleigh.Rayleigh.get_effective_wavelength')
- def test_get_reflectance(self, get_effective_wvl,
- get_rsr_version, get_reflectance_lut_from_file, download_luts, exists):
+ @patch('pyspectral.rayleigh.da', None)
+ def test_get_reflectance(self, fake_lut_hdf5):
"""Test getting the reflectance correction."""
- get_reflectance_lut_from_file.return_value = (self.rayl, self.wvl_coord, self.azid_coord,
- self.satz_sec_coord, self.sunz_sec_coord)
- download_luts.return_code = None
- exists.return_code = True
- get_rsr_version.return_code = RSR_DATA_VERSION
- get_effective_wvl.return_value = self.cwvl
-
sun_zenith = np.array([67., 32.])
sat_zenith = np.array([45., 18.])
azidiff = np.array([150., 110.])
redband_refl = np.array([14., 5.])
- refl_corr = self.viirs_rayleigh.get_reflectance(
- sun_zenith, sat_zenith, azidiff, 'M2', redband_refl)
+ rayl = _create_rayleigh()
+ with mocked_rsr():
+ refl_corr = rayl.get_reflectance(
+ sun_zenith, sat_zenith, azidiff, 'ch3', redband_refl)
np.testing.assert_allclose(refl_corr, TEST_RAYLEIGH_RESULT1)
sun_zenith = np.array([60., 20.])
sat_zenith = np.array([49., 26.])
azidiff = np.array([140., 130.])
redband_refl = np.array([12., 8.])
- refl_corr = self.viirs_rayleigh.get_reflectance(
- sun_zenith, sat_zenith, azidiff, 'M2', redband_refl)
+ with mocked_rsr():
+ refl_corr = rayl.get_reflectance(
+ sun_zenith, sat_zenith, azidiff, 'ch3', redband_refl)
+ assert isinstance(refl_corr, np.ndarray)
np.testing.assert_allclose(refl_corr, TEST_RAYLEIGH_RESULT2)
- @patch('pyspectral.rayleigh.HAVE_DASK', False)
- @patch('os.path.exists')
- @patch('pyspectral.utils.download_luts')
- @patch('pyspectral.rayleigh.get_reflectance_lut_from_file')
- def test_get_reflectance_no_rsr(self, get_reflectance_lut_from_file, download_luts, exists):
+ @patch('pyspectral.rayleigh.da', None)
+ def test_get_reflectance_no_rsr(self, fake_lut_hdf5):
"""Test getting the reflectance correction, simulating that we have no RSR data."""
- get_reflectance_lut_from_file.return_value = (self.rayl, self.wvl_coord, self.azid_coord,
- self.satz_sec_coord, self.sunz_sec_coord)
- download_luts.return_code = None
- exists.return_code = True
+ with mocked_rsr() as rsr_obj:
+ rsr_obj.side_effect = IOError("No rsr data in pyspectral for this platform and sensor")
+ sun_zenith = np.array([50., 10.])
+ sat_zenith = np.array([39., 16.])
+ azidiff = np.array([170., 110.])
+ redband_refl = np.array([29., 12.])
+ ufo = rayleigh.Rayleigh('UFO', 'unknown', atmosphere='midlatitude summer')
+ with pytest.raises(KeyError):
+ ufo.get_reflectance(sun_zenith, sat_zenith, azidiff, 'ch3', redband_refl)
- with patch('pyspectral.rayleigh.RelativeSpectralResponse') as mymock:
- instance = mymock.return_value
- mymock.side_effect = IOError("No rsr data in pyspectral for this platform and sensor")
- instance.rsr = None
- instance.unit = '1e-6 m'
- instance.si_scale = 1e-6
+ @patch('pyspectral.rayleigh.da', None)
+ def test_get_reflectance_float_wavelength(self, fake_lut_hdf5):
+ """Test getting the reflectance correction."""
+ with mocked_rsr() as rsr_obj:
+ rsr_obj.side_effect = IOError("No rsr data in pyspectral for this platform and sensor")
sun_zenith = np.array([50., 10.])
sat_zenith = np.array([39., 16.])
azidiff = np.array([170., 110.])
redband_refl = np.array([29., 12.])
- ufo = rayleigh.Rayleigh('UFO', 'unknown')
+ # rayl = _create_rayleigh()
+ rayl = rayleigh.Rayleigh('UFO', 'unknown', atmosphere='midlatitude summer')
- refl_corr = ufo.get_reflectance(sun_zenith, sat_zenith, azidiff, 0.441, redband_refl)
+ # we gave a direct wavelength so the RSR object shouldn't be needed
+ 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()
+
+ @patch('pyspectral.rayleigh.da', None)
+ def test_get_reflectance_wvl_outside_range(self, fake_lut_hdf5):
+ """Test getting the reflectance correction with wavelength outside correction range."""
+ with mocked_rsr() as rsr_obj:
+ rsr_obj.side_effect = IOError("No rsr data in pyspectral for this platform and sensor")
+ sun_zenith = np.array([50., 10.])
+ sat_zenith = np.array([39., 16.])
+ azidiff = np.array([170., 110.])
+ redband_refl = np.array([29., 12.])
+ rayl = rayleigh.Rayleigh('UFO', 'unknown', atmosphere='midlatitude summer')
+
+ # we gave a direct wavelength so the RSR object shouldn't be needed
+ 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()
+
+ def test_get_reflectance_no_lut(self, fake_lut_hdf5):
+ """Test that missing a LUT causes an exception.."""
+ # test LUT doesn't have a subartic_summer file
+ with pytest.raises(IOError):
+ rayleigh.Rayleigh('UFO', 'unknown', atmosphere='subarctic summer')
=====================================
pyspectral/utils.py
=====================================
@@ -1,12 +1,8 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
-# Copyright (c) 2014-2021 Pytroll developers
+# Copyright (c) 2014-2022 Pytroll developers
#
-# Author(s):
-#
-# Adam.Dybbroe <adam.dybbroe at smhi.se>
-# Panu Lahtinen <panu.lahtinen at fmi.fi>
#
# 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
@@ -25,7 +21,6 @@
import os
import logging
-import tempfile
import numpy as np
from pyspectral.config import get_config
from pyspectral.bandnames import BANDNAMES
@@ -70,10 +65,10 @@ INSTRUMENTS = {'NOAA-19': 'avhrr/3',
'MTG-I1': 'fci'
}
-HTTP_PYSPECTRAL_RSR = "https://zenodo.org/record/5797957/files/pyspectral_rsr_data.tgz"
+HTTP_PYSPECTRAL_RSR = "https://zenodo.org/record/6026563/files/pyspectral_rsr_data.tgz"
RSR_DATA_VERSION_FILENAME = "PYSPECTRAL_RSR_VERSION"
-RSR_DATA_VERSION = "v1.0.17"
+RSR_DATA_VERSION = "v1.0.18"
ATM_CORRECTION_LUT_VERSION = {}
ATM_CORRECTION_LUT_VERSION['antarctic_aerosol'] = {'version': 'v1.0.1',
@@ -116,22 +111,12 @@ for atype in AEROSOL_TYPES:
url = "{prefix}_{name}.tgz".format(prefix=URL_PREFIX, name=name)
HTTPS_RAYLEIGH_LUTS[atype] = url
-CONF = get_config()
-LOCAL_RSR_DIR = CONF.get('rsr_dir')
-LOCAL_RAYLEIGH_DIR = CONF.get('rayleigh_dir')
-
-try:
- os.makedirs(LOCAL_RSR_DIR)
-except OSError:
- if not os.path.isdir(LOCAL_RSR_DIR):
- raise
-RAYLEIGH_LUT_DIRS = {}
-for sub_dir_name in HTTPS_RAYLEIGH_LUTS:
- dirname = os.path.join(LOCAL_RAYLEIGH_DIR, sub_dir_name)
- RAYLEIGH_LUT_DIRS[sub_dir_name] = dirname
-
-TB2RAD_DIR = CONF.get('tb2rad_dir', tempfile.gettempdir())
+def get_rayleigh_lut_dir(aerosol_type):
+ """Get the rayleight LUT directory for the specified aerosol type."""
+ conf = get_config()
+ local_rayleigh_dir = conf.get('rayleigh_dir')
+ return os.path.join(local_rayleigh_dir, aerosol_type)
def convert2wavenumber(rsr):
@@ -296,7 +281,9 @@ def download_rsr(**kwargs):
except ImportError:
TQDM_LOADED = False
- dest_dir = kwargs.get('dest_dir', LOCAL_RSR_DIR)
+ config = get_config()
+ local_rsr_dir = config.get('rsr_dir')
+ dest_dir = kwargs.get('dest_dir', local_rsr_dir)
dry_run = kwargs.get('dry_run', False)
LOG.info("Download RSR files and store in directory %s", dest_dir)
@@ -344,14 +331,13 @@ def download_luts(**kwargs):
aerosol_types = HTTPS_RAYLEIGH_LUTS.keys()
chunk_size = 1024 * 1024 # 1 MB
-
for subname in aerosol_types:
LOG.debug('Aerosol type: %s', subname)
http = HTTPS_RAYLEIGH_LUTS[subname]
LOG.debug('URL = %s', http)
- subdir_path = RAYLEIGH_LUT_DIRS[subname]
+ subdir_path = get_rayleigh_lut_dir(subname)
try:
LOG.debug('Create directory: %s', subdir_path)
if not dry_run:
=====================================
rsr_convert_scripts/README.rst
=====================================
@@ -172,9 +172,13 @@ the pyspectral.yaml file:
filename: J1_VIIRS_Detector_RSR_V2/J1_VIIRS_RSR_DNBLGS_Detector_Fused_V2S.txt
bands: [DNB]
+.. code::
+
+ %> python msu_gs_reader.py
-Adam Dybbroe
-Sat Dec 1 17:39:48 2018
+Converts RSRs for the MSU-GS sensor aboard Electro-L N2. RSRs were retrieved from the NWP-SAF.
+Filenames look like:
+``rtcoef_electro-l_2_msugs_srf_ch01.txt``
.. code::
@@ -185,3 +189,6 @@ for FY-3B come as ``.prn`` text files for each channel (ex. ``ch1.prn``). For
FY-3C they come as ``.txt`` text files for channels 1, 2, 6, 7, 8, 9, and 10
only with names like ``FY3C_VIRR_CH01.txt``.
+PyTroll developers
+
+
=====================================
rsr_convert_scripts/msu_gs_reader.py
=====================================
@@ -0,0 +1,105 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Copyright (c) 2022 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 Electro-L N2 MSU-GS spectral response functions. Data from the NWPSAF:
+https://nwp-saf.eumetsat.int/downloads/rtcoef_rttov13/ir_srf/rtcoef_electro-l_2_msugs_srf.html
+
+"""
+
+import os
+import numpy as np
+from pyspectral.utils import convert2hdf5 as tohdf5
+from pyspectral.raw_reader import InstrumentRSR
+import logging
+
+LOG = logging.getLogger(__name__)
+
+MSUGS_BAND_NAMES = ['ch1', 'ch2', 'ch3', 'ch4', 'ch5',
+ 'ch6', 'ch7', 'ch8', 'ch9', 'ch10']
+
+#: 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 MsugsRSR(InstrumentRSR):
+ """Container for the Electro-L N2 MSU-GS relative spectral response data"""
+
+ def __init__(self, bandname, platform_name):
+
+ super(MsugsRSR, self).__init__(
+ bandname, platform_name, MSUGS_BAND_NAMES)
+
+ self.instrument = 'msu-gs'
+ self._get_options_from_config()
+ self._get_bandfilenames()
+
+ LOG.debug("Filenames: %s", str(self.filenames))
+ if self.filenames[bandname] and os.path.exists(self.filenames[bandname]):
+ self.requested_band_filename = self.filenames[bandname]
+ 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 MSU-GS RSR data for the band requested"""
+ data = np.genfromtxt(self.requested_band_filename,
+ unpack=True,
+ names=['wavenumber',
+ 'response'],
+ skip_header=4)
+
+ # Data are wavenumbers in cm-1:
+ wavelength = 1. / data[0][::-1] * scale
+ response = data[1][::-1]
+
+ # The real MSU-GS likely has multiple detectors. However, for now
+ # we store the single rsr available in the NWPSAF coefficient files:
+ self.rsr = {'wavelength': wavelength, 'response': response}
+
+
+def main():
+ """Main"""
+ for platform_name in ['Electro-L-N2', ]:
+ tohdf5(MsugsRSR, platform_name, MSUGS_BAND_NAMES)
+
+
+if __name__ == "__main__":
+ import sys
+
+ LOG = logging.getLogger('msu_gs_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)
+
+ main()
View it on GitLab: https://salsa.debian.org/debian-gis-team/pyspectral/-/commit/98af4161de0345eb30c8fe3a3a3ccc264a940b31
--
View it on GitLab: https://salsa.debian.org/debian-gis-team/pyspectral/-/commit/98af4161de0345eb30c8fe3a3a3ccc264a940b31
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/20220306/42dd5a53/attachment-0001.htm>
More information about the Pkg-grass-devel
mailing list