CHANGELOG.md

+## Version <RELEASE_VERSION> (2019/08/30)
+
+### Issues Closed
+
+* [Issue 73](https://github.com/pytroll/pyspectral/issues/73) - Fix blackbody code to work with dask arrays ([PR 74](https://github.com/pytroll/pyspectral/pull/74))
+
+In this release 1 issue was closed.
+
+### Pull Requests Merged
+
+#### Bugs fixed
+
+* [PR 80](https://github.com/pytroll/pyspectral/pull/80) - Fix doc tests for python 2&3
+* [PR 79](https://github.com/pytroll/pyspectral/pull/79) - Fix rsr zenodo version
+* [PR 74](https://github.com/pytroll/pyspectral/pull/74) - Fix dask compatibility in blackbody functions ([73](https://github.com/pytroll/pyspectral/issues/73))
+
+#### Features added
+
+* [PR 78](https://github.com/pytroll/pyspectral/pull/78) - Add FY-3B VIRR and FY-3C VIRR RSRs
+* [PR 77](https://github.com/pytroll/pyspectral/pull/77) - Add FY-4A AGRI support
+
+In this release 5 pull requests were closed.
+
+
 ## Version <RELEASE_VERSION> (2019/06/07)
 
 ### Issues Closed

README.md

 PySpectral
 ==========
 
-[![Codacy Badge](https://api.codacy.com/project/badge/Grade/9f039d7d640846ca89be8a78fa11e1f6)](https://www.codacy.com/app/adybbroe/pyspectral?utm_source=github.com&utm_medium=referral&utm_content=pytroll/pyspectral&utm_campaign=badger)
 [![Build Status](https://travis-ci.org/pytroll/pyspectral.png?branch=master)](https://travis-ci.org/pytroll/pyspectral)
 [![Build status](https://ci.appveyor.com/api/projects/status/5lm42n0l65l5o9xn?svg=true)](https://ci.appveyor.com/project/pytroll/pyspectral)
 [![Coverage Status](https://coveralls.io/repos/github/pytroll/pyspectral/badge.svg?branch=master)](https://coveralls.io/github/pytroll/pyspectral?branch=master)
-[![Code Health](https://landscape.io/github/pytroll/pyspectral/master/landscape.png)](https://landscape.io/github/pytroll/pyspectral/master)
 [![PyPI version](https://badge.fury.io/py/pyspectral.svg)](https://badge.fury.io/py/pyspectral)
 [![Code Climate](https://codeclimate.com/github/pytroll/pyspectral/badges/gpa.svg)](https://codeclimate.com/github/pytroll/pyspectral)
-[![Scrutinizer Code Quality](https://scrutinizer-ci.com/g/pytroll/pyspectral/badges/quality-score.png?b=master)](https://scrutinizer-ci.com/g/pytroll/pyspectral/?branch=master)
 
 Given a passive sensor on a meteorological satellite PySpectral provides the
 relative spectral response (rsr) function(s) and offer you some basic

debian/changelog

+pyspectral (0.9.0+ds-1) unstable; urgency=medium
+
+  * New upstream release.
+
+ -- Antonio Valentino <antonio.valentino@tiscali.it>  Mon, 02 Sep 2019 06:04:23 +0000
+
 pyspectral (0.8.9+ds-2) unstable; urgency=medium
 
   * Bump Standards-Version to 4.4.0, no changes.

doc/37_reflectance.rst

@@ -46,7 +46,7 @@ expressed in :math:`W/m^2 sr^{-1} \mu m^{-1}`, or using SI units :math:`W/m^2 sr
   >>> viirs = RadTbConverter('Suomi-NPP', 'viirs', 'M12')
   >>> rad37 = viirs.tb2radiance(tb37)
   >>> print([np.round(rad, 7) for rad in rad37['radiance']])
-  [369717.4765726, 355110.5207853, 314684.2788726, 173143.5424898, 116408.0007877]
+  [369717.4972296, 355110.6414922, 314684.3507084, 173143.4836477, 116408.0022674]
   >>> rad37['unit']
   'W/m^2 sr^-1 m^-1'
 
@@ -59,7 +59,7 @@ In order to get the total radiance over the band one has to multiply with the eq
   >>> viirs = RadTbConverter('Suomi-NPP', 'viirs', 'M12')
   >>> rad37 = viirs.tb2radiance(tb37, normalized=False)
   >>> print([np.round(rad, 8) for rad in rad37['radiance']])
-  [0.07037968, 0.06759909, 0.05990352, 0.03295972, 0.02215951]
+  [0.07037968, 0.06759911, 0.05990353, 0.03295971, 0.02215951]
   >>> rad37['unit']
   'W/m^2 sr^-1'
 
@@ -218,17 +218,17 @@ We can try decompose equation :eq:`refl37` above using the example of VIIRS M12
   >>> print(np.isnan(nomin))
   [False False False False False]
   >>> print([np.round(val, 8) for val in nomin])
-  [0.05083677, 0.04805618, 0.0404157, 0.01279279, 0.00204485]
+  [0.05083677, 0.0480562, 0.04041571, 0.01279277, 0.00204485]
   >>> denom = np.cos(np.deg2rad(sunz))/np.pi * sflux - rad11['radiance']
   >>> print(np.isnan(denom))
   [False False False False False]
   >>> print([np.round(val, 8) for val in denom])
-  [0.23646313, 0.23645682, 0.23650559, 0.23582015, 0.2358661]
+  [0.23646312, 0.23645681, 0.23650559, 0.23582014, 0.23586609]
   >>> res = nomin/denom
   >>> print(np.isnan(res))
   [False False False False False]
   >>> print([np.round(val, 8) for val in res])
-  [0.21498817, 0.2032345, 0.17088689, 0.05424807, 0.00866955]
+  [0.21498817, 0.20323458, 0.17088693, 0.05424801, 0.00866952]
 
 
 Derive the emissive part of the 3.7 micron band
@@ -255,5 +255,5 @@ Using the example of the VIIRS M12 band from above this gives the following spec
   >>> ['{tb:6.3f}'.format(tb=np.round(t, 4)) for t in tb]
   ['266.996', '267.262', '267.991', '271.033', '271.927']
   >>> rad = refl_m12.emissive_part_3x(tb=False)
-  >>> ['{rad:6.3f}'.format(rad=np.round(r, 3)) for r in rad]
-  ['80285.149', '81458.022', '84749.639', '99761.400', '104582.030']
+  >>> ['{rad:6.3f}'.format(rad=np.round(r, 3)) for r in rad.compute()]
+  ['80285.150', '81458.022', '84749.638', '99761.401', '104582.031']

doc/rad_definitions.rst

@@ -229,7 +229,7 @@ And using wavelength representation:
    >>> wvl = 1./wavenumber
    >>> rad = blackbody(wvl, [300., 301])
    >>> print("{0:10.3f} {1:10.3f}".format(rad[0], rad[1]))
-   9573178.886 9714689.259
+   9573177.494 9714687.157
 
 Which are the spectral radiances in SI units around :math:`11 \mu m` at
 temperatures 300 and 301 Kelvin. In units of :math:`mW/m^2\ m^{-1} sr^{-1}` this becomes:

pyspectral/blackbody.py

 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 
-# Copyright (c) 2013-2018 Adam.Dybbroe
+# Copyright (c) 2013-2019 Adam.Dybbroe
 
 # Author(s):
 
@@ -23,8 +23,13 @@
 
 """Planck radiation equation"""
 import numpy as np
-
 import logging
+
+try:
+    import dask.array as da
+except ImportError:
+    da = np
+
 LOG = logging.getLogger(__name__)
 
 H_PLANCK = 6.62606957 * 1e-34  # SI-unit = [J*s]
@@ -84,37 +89,38 @@ def blackbody_wn_rad2temp(wavenumber, radiance):
     function. Wavenumber space"""
 
     if np.isscalar(radiance):
-        rad = np.array([radiance, ], dtype='float64')
-    else:
-        rad = np.array(radiance, dtype='float64')
+        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')
-    else:
+        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, rad) + 1.0)
+    res = const1 * wavnum / np.log(
+        np.divide(const2 * wavnum**3, radiance) + 1.0)
 
-    shape = rad.shape
+    shape = radiance.shape
     resshape = res.shape
 
     if wavnum.shape[0] == 1:
-        if rad.shape[0] == 1:
+        if radiance.shape[0] == 1:
             return res[0]
         else:
             return res[::].reshape(shape)
     else:
-        if rad.shape[0] == 1:
+        if radiance.shape[0] == 1:
             return res[0, :]
         else:
             if len(shape) == 1:
-                return np.reshape(res, (shape[0], resshape[1]))
+                return res.reshape((shape[0], resshape[1]))
             else:
-                return np.reshape(res, (shape[0], shape[1], resshape[1]))
+                return res.reshape((shape[0], shape[1], resshape[1]))
 
 
-def planck(wave, temp, wavelength=True):
+def planck(wave, temperature, wavelength=True):
     """The Planck radiation or Blackbody radiation as a function of wavelength
     or wavenumber. SI units.
     _planck(wave, temperature, wavelength=True)
@@ -136,12 +142,12 @@ def planck(wave, temp, wavelength=True):
     units = ['wavelengths', 'wavenumbers']
     if wavelength:
         LOG.debug("Using {0} when calculating the Blackbody radiance".format(
-            units[(wavelength == True) - 1]))
+            units[(wavelength is True) - 1]))
 
-    if np.isscalar(temp):
-        temperature = np.array([temp, ], dtype='float64')
-    else:
-        temperature = np.array(temp, dtype='float64')
+    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):
@@ -157,13 +163,19 @@ def planck(wave, temp, wavelength=True):
         nom = 2 * H_PLANCK * (C_SPEED ** 2) * (wln ** 3)
         arg1 = H_PLANCK * C_SPEED * wln / K_BOLTZMANN
 
-    arg2 = np.where(np.greater(np.abs(temperature), EPSILON),
-                    np.array(1. / temperature), -9).reshape(-1, 1)
-    arg2 = np.ma.masked_array(arg2, mask=arg2 == -9)
-    LOG.debug("Max and min - arg1: %s  %s", str(arg1.max()), str(arg1.min()))
-    LOG.debug("Max and min - arg2: %s  %s", str(arg2.max()), str(arg2.min()))
+    # use dask functions when needed
+    np_ = np if isinstance(temperature, np.ndarray) else da
+    arg2 = np_.where(np.greater(np.abs(temperature), EPSILON),
+                     (1. / temperature), np.nan).reshape(-1, 1)
+    if isinstance(arg2, np.ndarray):
+        # don't compute min/max if we have dask arrays
+        LOG.debug("Max and min - arg1: %s  %s",
+                  str(np.nanmax(arg1)), str(np.nanmin(arg1)))
+        LOG.debug("Max and min - arg2: %s  %s",
+                  str(np.nanmax(arg2)), str(np.nanmin(arg2)))
+
     try:
-        exp_arg = np.multiply(arg1.astype('float32'), arg2.astype('float32'))
+        exp_arg = np.multiply(arg1.astype('float64'), arg2.astype('float64'))
     except MemoryError:
         LOG.warning(("Dimensions used in numpy.multiply probably reached "
                      "limit!\n"
@@ -171,9 +183,9 @@ def planck(wave, temp, wavelength=True):
                      "and try running again"))
         raise
 
-    LOG.debug("Max and min before exp: %s  %s", str(exp_arg.max()),
-              str(exp_arg.min()))
-    if exp_arg.min() < 0:
+    if isinstance(exp_arg, np.ndarray) and exp_arg.min() < 0:
+        LOG.debug("Max and min before exp: %s  %s",
+                  str(exp_arg.max()), str(exp_arg.min()))
         LOG.warning("Something is fishy: \n" +
                     "\tDenominator might be zero or negative in radiance derivation:")
         dubious = np.where(exp_arg < 0)[0]
@@ -182,7 +194,6 @@ def planck(wave, temp, wavelength=True):
 
     denom = np.exp(exp_arg) - 1
     rad = nom / denom
-    rad = np.where(rad.mask, np.nan, rad.data)
     radshape = rad.shape
     if wln.shape[0] == 1:
         if temperature.shape[0] == 1:
@@ -194,9 +205,9 @@ def planck(wave, temp, wavelength=True):
             return rad[0, :]
         else:
             if len(shape) == 1:
-                return np.reshape(rad, (shape[0], radshape[1]))
+                return rad.reshape((shape[0], radshape[1]))
             else:
-                return np.reshape(rad, (shape[0], shape[1], radshape[1]))
+                return rad.reshape((shape[0], shape[1], radshape[1]))
 
 
 def blackbody_wn(wavenumber, temp):

pyspectral/config.py

@@ -28,6 +28,11 @@ import os
 from os.path import expanduser
 from appdirs import AppDirs
 import yaml
+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
 

pyspectral/etc/pyspectral.yaml

@@ -209,6 +209,46 @@ download_from_internet: True
 #   ch15: GOES-R_ABI_FM2_SRF_CWG_ch15.txt
 #   ch16: GOES-R_ABI_FM2_SRF_CWG_ch16.txt
 
+# FY-4A-agri:
+#     path: /path/to/original/fy4a/agri/data
+#     ch1: FY4A_AGRI_SRF_CH01.txt
+#     ch2: FY4A_AGRI_SRF_CH02.txt
+#     ch3: FY4A_AGRI_SRF_CH03.txt
+#     ch4: FY4A_AGRI_SRF_CH04.txt
+#     ch5: FY4A_AGRI_SRF_CH05.txt
+#     ch6: FY4A_AGRI_SRF_CH06.txt
+#     ch7: FY4A_AGRI_SRF_CH07.txt
+#     ch8: FY4A_AGRI_SRF_CH08.txt
+#     ch9: FY4A_AGRI_SRF_CH09.txt
+#     ch10: FY4A_AGRI_SRF_CH10.txt
+#     ch11: FY4A_AGRI_SRF_CH11.txt
+#     ch12: FY4A_AGRI_SRF_CH12.txt
+#     ch13: FY4A_AGRI_SRF_CH13.txt
+#     ch14: FY4A_AGRI_SRF_CH14.txt
+
+#FY-3B-virr:
+#  path: /Users/davidh/repos/git/pyspectral/virr_srf/FY3B-VIRR
+#  ch1: ch1.prn
+#  ch2: ch2.prn
+#  ch3: ch3.prn
+#  ch4: ch4.prn
+#  ch5: ch5.prn
+#  ch6: ch6.prn
+#  ch7: ch7.prn
+#  ch8: ch8.prn
+#  ch9: ch9.prn
+#  ch10: ch10.prn
+
+#FY-3C-virr:
+#  path: /Users/davidh/repos/git/pyspectral/virr_srf/FY3C_VIRR_SRF
+#  ch1: FY3C_VIRR_CH01.txt
+#  ch2: FY3C_VIRR_CH02.txt
+#  ch6: FY3C_VIRR_CH06.txt
+#  ch7: FY3C_VIRR_CH07.txt
+#  ch8: FY3C_VIRR_CH08.txt
+#  ch9: FY3C_VIRR_CH09.txt
+#  ch10: FY3C_VIRR_CH10.txt
+
 # FY-3D-mersi-2:
 #   path: /path/to/original/fy3d/mersi2/data
 #   ch1:  FY3D_MERSI_SRF_CH01_Pub.txt

pyspectral/rsr_reader.py

@@ -27,16 +27,15 @@ import os
 import numpy as np
 from glob import glob
 from os.path import expanduser
-
-import logging
-LOG = logging.getLogger(__name__)
-
 from pyspectral.config import get_config
 from pyspectral.utils import WAVE_NUMBER
 from pyspectral.utils import WAVE_LENGTH
 from pyspectral.utils import (INSTRUMENTS, download_rsr)
 from pyspectral.utils import (RSR_DATA_VERSION_FILENAME, RSR_DATA_VERSION)
 
+import logging
+LOG = logging.getLogger(__name__)
+
 
 class RSRDataBaseClass(object):
 
@@ -168,11 +167,18 @@ class RelativeSpectralResponse(RSRDataBaseClass):
 
         no_detectors_message = False
         with h5py.File(self.filename, 'r') as h5f:
-            self.band_names = [b.decode('utf-8') for b in h5f.attrs['band_names'].tolist()]
-            self.description = h5f.attrs['description'].decode('utf-8')
+            self.band_names = h5f.attrs['band_names'].tolist()
+            self.description = h5f.attrs['description']
+            if not isinstance(self.band_names[0], str):
+                # byte array in python 3
+                self.band_names = [x.decode('utf-8') for x in self.band_names]
+                self.description = self.description.decode('utf-8')
+
             if not self.platform_name:
                 try:
-                    self.platform_name = h5f.attrs['platform_name'].decode('utf-8')
+                    self.platform_name = h5f.attrs['platform_name']
+                    if not isinstance(self.platform_name, str):
+                        self.platform_name = self.platform_name.decode('utf-8')
                 except KeyError:
                     LOG.warning("No platform_name in HDF5 file")
                     try:
@@ -186,6 +192,8 @@ class RelativeSpectralResponse(RSRDataBaseClass):
             if not self.instrument:
                 try:
                     self.instrument = h5f.attrs['sensor'].decode('utf-8')
+                    if not isinstance(self.instrument, str):
+                        self.instrument = self.instrument.decode('utf-8')
                 except KeyError:
                     LOG.warning("No sensor name specified in HDF5 file")
                     self.instrument = INSTRUMENTS.get(self.platform_name)

pyspectral/tests/test_atm_correction_ir.py

 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 
-# Copyright (c) 2017 Adam.Dybbroe
+# Copyright (c) 2017 - 2019 Pytroll
 
 # Author(s):
 
-#   Adam.Dybbroe <a000680@c20671.ad.smhi.se>
+#   Adam.Dybbroe <adam.dybbroe@smhi.se>
 
 # 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
@@ -19,27 +19,17 @@
 
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
-
-"""Unit tests of the atmospherical correction in the ir spectral range
-"""
+"""Unit tests of the atmospherical correction in the ir spectral range."""
 
 
+import numpy as np
+from pyspectral.atm_correction_ir import AtmosphericalCorrection
 import sys
 if sys.version_info < (2, 7):
     import unittest2 as unittest
 else:
     import unittest
 
-import numpy as np
-from pyspectral.atm_correction_ir import AtmosphericalCorrection
-#from mock import patch
-
-from pyspectral.tests.unittest_helpers import assertNumpyArraysEqual
-
-# Mock some modules, so we don't need them for tests.
-
-#sys.modules['pyresample'] = MagicMock()
-
 SATZ = np.ma.array([[48.03,  48.03002,  48.03004,  48.03006,  48.03008,  48.0301,
                      48.03012,  48.03014,  48.03016,  48.03018],
                     [48.09,  48.09002,  48.09004,  48.09006,  48.09008,  48.0901,
@@ -125,28 +115,18 @@ RES = np.ma.array([[286.03159412,  286.03162417,  286.03165421,  286.03168426,
 
 
 class TestAtmCorrection(unittest.TestCase):
-
-    """Class for testing pyspectral.atm_correction_ir"""
-
-    def setUp(self):
-        """Setup the test"""
-        pass
+    """Class for testing pyspectral.atm_correction_ir."""
 
     def test_get_correction(self):
         """Test getting the atm correction"""
 
         this = AtmosphericalCorrection('EOS-Terra', 'modis')
         atm_corr = this.get_correction(SATZ, None, TBS)
-        assertNumpyArraysEqual(TBS, atm_corr)
-
-    def tearDown(self):
-        """Clean up"""
-        pass
+        np.testing.assert_almost_equal(TBS, atm_corr)
 
 
 def suite():
-    """The test suite for test_atm_correction_ir.
-    """
+    """Create the test suite for test_atm_correction_ir."""
     loader = unittest.TestLoader()
     mysuite = unittest.TestSuite()
     mysuite.addTest(loader.loadTestsFromTestCase(TestAtmCorrection))

pyspectral/tests/test_blackbody.py

 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 
-# Copyright (c) 2013, 2014, 2015, 2016, 2017 Adam.Dybbroe
+# Copyright (c) 2013 - 2019 Pytroll
 
 # Author(s):
 
-#   Adam.Dybbroe <a000680@c14526.ad.smhi.se>
+#   Adam.Dybbroe <adam.dybbroe@smhi.se>
 
 # 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
@@ -31,10 +31,8 @@ from pyspectral.tests.unittest_helpers import assertNumpyArraysEqual
 import unittest
 import numpy as np
 
-#RAD_11MICRON_300KELVIN = 9572498.1141643394
-RAD_11MICRON_300KELVIN = 9573177.8811719529
-#RAD_11MICRON_301KELVIN = 9713997.9623772576
-RAD_11MICRON_301KELVIN = 9714688.2959563732
+RAD_11MICRON_300KELVIN = 9573176.935507433
+RAD_11MICRON_301KELVIN = 9714686.576498277
 
 # Radiances in wavenumber space (SI-units)
 WN_RAD_11MICRON_300KELVIN = 0.00115835441353
@@ -43,18 +41,28 @@ WN_RAD_11MICRON_301KELVIN = 0.00117547716523
 __unittest = True
 
 
-class TestBlackbody(unittest.TestCase):
+class CustomScheduler(object):
+    """Custom dask scheduler that raises an exception if dask is computed too many times."""
 
-    """Unit testing the blackbody function"""
+    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)
 
-    def setUp(self):
-        """Set up"""
-        return
+
+class TestBlackbody(unittest.TestCase):
+    """Unit testing the blackbody function"""
 
     def test_blackbody(self):
-        """Calculate the blackbody radiation from wavelengths and
-        temperatures
-        """
+        """Calculate the blackbody radiation from wavelengths and temperatures."""
         wavel = 11. * 1E-6
         black = blackbody((wavel, ), [300., 301])
         self.assertEqual(black.shape[0], 2)
@@ -71,14 +79,28 @@ class TestBlackbody(unittest.TestCase):
 
         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)
 
         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):
+        """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)
 
     def test_blackbody_wn(self):
-        """Calculate the blackbody radiation from wavenumbers and
-        temperatures
-        """
+        """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)
@@ -106,9 +128,24 @@ class TestBlackbody(unittest.TestCase):
 
         assertNumpyArraysEqual(t__, expected)
 
-    def tearDown(self):
-        """Clean up"""
-        return
+    def test_blackbody_wn_dask(self):
+        """Test that blackbody rad2temp preserves dask arrays."""
+        import dask
+        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)):
+            t__ = blackbody_wn_rad2temp(wavenumber, radiances)
+        self.assertIsInstance(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)
 
 
 def suite():

pyspectral/tests/test_rad_tb_conversions.py

 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 
-# Copyright (c) 2014-2018 Adam.Dybbroe
+# Copyright (c) 2014-2019 Adam.Dybbroe
 
 # Author(s):
 
@@ -347,10 +347,10 @@ class TestRadTbConversions(unittest.TestCase):
         self.assertTrue(np.allclose(TRUE_RADS * integral, res['radiance']))
 
         res = self.modis.tb2radiance(237., lut=False)
-        self.assertAlmostEqual(16570.592171157, res['radiance'])
+        self.assertAlmostEqual(16570.579551068, res['radiance'])
 
         res = self.modis.tb2radiance(277., lut=False)
-        self.assertAlmostEqual(167544.3823631, res['radiance'])
+        self.assertAlmostEqual(167544.39368663222, res['radiance'])
 
         res = self.modis.tb2radiance(1.1, lut=False)
         self.assertAlmostEqual(0.0, res['radiance'])
@@ -362,7 +362,7 @@ class TestRadTbConversions(unittest.TestCase):
         self.assertAlmostEqual(5.3940515573e-06, res['radiance'])
 
         res = self.modis.tb2radiance(200.1, lut=False)
-        self.assertAlmostEqual(865.09776189, res['radiance'])
+        self.assertAlmostEqual(865.09759706, res['radiance'])
 
     def tearDown(self):
         """Clean up"""

pyspectral/tests/test_reflectance.py

@@ -139,7 +139,10 @@ class TestReflectance(unittest.TestCase):
 
         with patch('pyspectral.radiance_tb_conversion.RelativeSpectralResponse') as mymock:
             instance = mymock.return_value
-            instance.rsr = TEST_RSR
+            # VIIRS doesn't have a channel '20' like MODIS so the generic
+            # mapping this test will end up using will find 'ch20' for VIIRS
+            viirs_rsr = {'ch20': TEST_RSR['20'], '99': TEST_RSR['99']}
+            instance.rsr = viirs_rsr
             instance.unit = '1e-6 m'
             instance.si_scale = 1e-6
 
@@ -148,7 +151,7 @@ class TestReflectance(unittest.TestCase):
 
             refl37 = Calculator('Suomi-NPP', 'viirs', 3.7)
             self.assertEqual(refl37.bandwavelength, 3.7)
-            self.assertEqual(refl37.bandname, '20')
+            self.assertEqual(refl37.bandname, 'ch20')
 
         with patch('pyspectral.radiance_tb_conversion.RelativeSpectralResponse') as mymock:
             instance = mymock.return_value

pyspectral/tests/test_utils.py

@@ -115,8 +115,7 @@ class TestUtils(unittest.TestCase):
         self.rsr = RsrTestData()
 
     def test_convert2wavenumber(self):
-        """Testing the conversion of rsr from wavelength to wavenumber
-        """
+        """Testing the conversion of rsr from wavelength to wavenumber."""
         newrsr, info = utils.convert2wavenumber(TEST_RSR)
         unit = info['unit']
         self.assertEqual(unit, 'cm-1')
@@ -127,11 +126,7 @@ class TestUtils(unittest.TestCase):
         self.assertTrue(np.allclose(wvn_res, wvn))
 
     def test_get_bandname_from_wavelength(self):
-        """Test that it is possible to get the right bandname provided the wavelength
-        in micro meters
-
-        """
-
+        """Test the right bandname is found provided the wavelength in micro meters."""
         x = utils.get_bandname_from_wavelength('abi', 0.4, self.rsr.rsr)
         self.assertEqual(x, 'ch1')
         with self.assertRaises(AttributeError):
@@ -146,19 +141,16 @@ class TestUtils(unittest.TestCase):
         x = utils.get_bandname_from_wavelength('abi', 1.0, self.rsr.rsr)
         self.assertEqual(x, None)
 
+        # uses generic channel mapping where '20' -> 'ch20'
         bandname = utils.get_bandname_from_wavelength('abi', 3.7, TEST_RSR)
-        self.assertEqual(bandname, '20')
+        self.assertEqual(bandname, 'ch20')
 
         bandname = utils.get_bandname_from_wavelength('abi', 3.0, TEST_RSR)
         self.assertIsNone(bandname)
 
-    def tearDown(self):
-        """Clean up"""
-        pass
-
 
 def suite():
-    """The suite for test_utils."""
+    """Create the suite for test_utils."""
     loader = unittest.TestLoader()
     mysuite = unittest.TestSuite()
     mysuite.addTest(loader.loadTestsFromTestCase(TestUtils))

pyspectral/utils.py

@@ -78,9 +78,19 @@ BANDNAMES['generic'] = {'VIS006': 'VIS0.6',
                         'C15': 'ch15',
                         'C16': 'ch16',
                         }
+# handle arbitrary channel numbers
+for chan_num in range(1, 37):
+    BANDNAMES['generic'][str(chan_num)] = 'ch{:d}'.format(chan_num)
 
-BANDNAMES['avhrr-3'] = {'3b': 'ch3b',
-                        '3a': 'ch3a'}
+# 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['avhrr-3'] = {'1': 'ch1',
+                        '2': 'ch2',
+                        '3b': 'ch3b',
+                        '3a': 'ch3a',
+                        '4': 'ch4',
+                        '5': 'ch5'}
 
 BANDNAMES['ahi'] = {'B01': 'ch1',
                     'B02': 'ch2',
@@ -120,12 +130,16 @@ INSTRUMENTS = {'NOAA-19': 'avhrr/3',
                'Suomi-NPP': 'viirs',
                'NOAA-20': 'viirs',
                'FY-3D': 'mersi-2',
-               'Feng-Yun 3D': 'mersi-2'
+               'FY-3C': 'virr',
+               'FY-3B': 'virr',
+               'Feng-Yun 3D': 'mersi-2',
+               'FY-4A': 'agri'
                }
 
-HTTP_PYSPECTRAL_RSR = "https://zenodo.org/record/2653487/files/pyspectral_rsr_data.tgz"
+
+HTTP_PYSPECTRAL_RSR = "https://zenodo.org/record/3381130/files/pyspectral_rsr_data.tgz"
 RSR_DATA_VERSION_FILENAME = "PYSPECTRAL_RSR_VERSION"
-RSR_DATA_VERSION = "v1.0.6"
+RSR_DATA_VERSION = "v1.0.9"
 
 ATM_CORRECTION_LUT_VERSION = {}
 ATM_CORRECTION_LUT_VERSION['antarctic_aerosol'] = {'version': 'v1.0.1',

pyspectral/version.py

@@ -46,9 +46,9 @@ def get_keywords():
     # setup.py/versioneer.py will grep for the variable names, so they must
     # each be defined on a line of their own. _version.py will just call
     # get_keywords().
-    git_refnames = " (tag: v0.8.9)"
-    git_full = "b3a1f58f49e2ab7c81c331e379065f24eb21f1ba"
-    git_date = "2019-06-07 09:27:42 +0200"
+    git_refnames = " (HEAD -> master, tag: v0.9.0)"
+    git_full = "2372396547ecace7fbcefc50d5dcbaa32e4a5177"
+    git_date = "2019-08-30 18:04:19 +0200"
     keywords = {"refnames": git_refnames, "full": git_full, "date": git_date}
     return keywords
 

rsr_convert_scripts/README.rst

@@ -176,3 +176,12 @@ the pyspectral.yaml file:
 Adam Dybbroe
 Sat Dec  1 17:39:48 2018
 
+.. code::
+
+    %> python virr_rsr.py
+
+Converting the FY-3B or FY-3C VIRR spectral responses to HDF5. Original files
+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``.
+

rsr_convert_scripts/agri_rsr.py

+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+# Copyright (c) 2018, 2019 Pytroll
+
+# Author(s):
+
+#   Xin.Zhang <xinzhang1215@gmail.com>
+#   Adam.Dybbroe <adam.dybbroe@smhi.se>
+
+# 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 FY-4A AGRI relative spectral responses. Data from
+http://fy4.nsmc.org.cn/portal/cn/fycv/srf.html
+"""
+import os
+import numpy as np
+from pyspectral.utils import INSTRUMENTS
+from pyspectral.utils import convert2hdf5 as tohdf5
+from pyspectral.raw_reader import InstrumentRSR
+from pyspectral.utils import logging_on, get_logger
+
+FY4A_BAND_NAMES = ['ch1', 'ch2', 'ch3', 'ch4', 'ch5', 'ch6', 'ch7', 'ch8',
+                   'ch9', 'ch10', 'ch11', 'ch12', 'ch13', 'ch14']
+BANDNAME_SCALE2MICROMETERS = {'ch1': 0.001,
+                              'ch2': 0.001,
+                              'ch3': 0.001,
+                              'ch4': 1.0,
+                              'ch5': 1.0,
+                              'ch6': 1.0,
+                              'ch7': 1.0,
+                              'ch8': 1.0,
+                              'ch9': 1.0,
+                              'ch10': 1.0,
+                              'ch11': 1.0,
+                              'ch12': 1.0,
+                              'ch13': 1.0,
+                              'ch14': 1.0}
+
+
+class AGRIRSR(InstrumentRSR):
+    """Container for the FY-4 AGRI RSR data"""
+
+    def __init__(self, bandname, platform_name):
+        """Initialise the FY-4 AGRI relative spectral response data"""
+        super(AGRIRSR, self).__init__(bandname, platform_name, FY4A_BAND_NAMES)
+
+        self.instrument = INSTRUMENTS.get(platform_name, 'agri')
+
+        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]
+            scale = BANDNAME_SCALE2MICROMETERS.get(bandname)
+            if scale:
+                self._load(scale=scale)
+            else:
+                LOG.error(
+                    "Failed determine the scale used to convert to wavelength in micrometers - channel = %s", bandname)
+                raise AttributeError('no scale for bandname %s', bandname)
+
+        else:
+            LOG.warning("Couldn't find an existing file for this band: %s",
+                        str(self.bandname))
+
+        self.filename = self.requested_band_filename
+
+    def _load(self, scale=0.001):
+        """Load the AGRI RSR data for the band requested
+
+           Wavelength is given in nanometers.
+        """
+        data = np.genfromtxt(self.requested_band_filename,
+                             unpack=True,
+                             names=['wavelength',
+                                    'response'],
+                             skip_header=0)
+
+        wavelength = data['wavelength'] * scale
+        response = data['response']
+
+        self.rsr = {'wavelength': wavelength, 'response': response}
+
+
+def main():
+    """Main"""
+    for platform_name in ["FY-4A", ]:
+        tohdf5(AGRIRSR, platform_name, FY4A_BAND_NAMES)
+
+
+if __name__ == "__main__":
+    LOG = get_logger(__name__)
+    logging_on()
+
+    main()

rsr_convert_scripts/virr_rsr.py

+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Copyright (c) 2019 Adam.Dybbroe
+#
+# Author(s):
+#
+#   David Hoese <david.hoese@ssec.wisc.edu>
+#
+# 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 VIRR relative spectral responses.
+
+Data from http://gsics.nsmc.org.cn/portal/en/fycv/srf.html
+
+"""
+
+import os
+import numpy as np
+from pyspectral.utils import INSTRUMENTS
+from pyspectral.utils import convert2hdf5 as tohdf5
+from pyspectral.raw_reader import InstrumentRSR
+
+import logging
+LOG = logging.getLogger(__name__)
+
+VIRR_BAND_NAMES = {
+    'FY-3B': ['ch{:d}'.format(x) for x in range(1, 11)],
+    'FY-3C': ['ch1', 'ch2'] + ['ch{:d}'.format(x) for x in range(6, 11)],
+}
+
+
+class VirrRSR(InstrumentRSR):
+    """Container for the FY-3B/FY-3C VIRR RSR data."""
+
+    def __init__(self, bandname, platform_name):
+        """Verify that file exists and can be read."""
+        super(VirrRSR, self).__init__(bandname, platform_name, VIRR_BAND_NAMES[platform_name])
+
+        self.instrument = INSTRUMENTS.get(platform_name, 'virr')
+        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
+
+    def _load(self, scale=0.001):
+        """Load the VIRR RSR data for the band requested.
+
+        Wavelength is given in nanometers.
+        """
+        data = np.genfromtxt(self.requested_band_filename,
+                             unpack=True,
+                             names=['wavelength',
+                                    'response'],
+                             skip_header=0)
+
+        wavelength = data['wavelength'] * scale
+        response = data['response']
+
+        self.rsr = {'wavelength': wavelength, 'response': response}
+
+
+def main():
+    """Main"""
+    for platform_name, band_names in VIRR_BAND_NAMES.items():
+        tohdf5(VirrRSR, platform_name, band_names)
+
+
+if __name__ == "__main__":
+    main()