CHANGELOG.md

+## Version <RELEASE_VERSION> (2019/04/09)
+
+### Issues Closed
+
+* [Issue 66](https://github.com/pytroll/pyspectral/issues/66) - Throws a warning about non-existing directory - storing/reading cached radiance-tb look-up-tables ([PR 67](https://github.com/pytroll/pyspectral/pull/67))
+* [Issue 61](https://github.com/pytroll/pyspectral/issues/61) - can this program be used for user-defined sensor or rsr？ ([PR 62](https://github.com/pytroll/pyspectral/pull/62))
+* [Issue 58](https://github.com/pytroll/pyspectral/issues/58) - Use dask instead of numpy masked arrays ([PR 59](https://github.com/pytroll/pyspectral/pull/59))
+
+In this release 3 issues were closed.
+
+### Pull Requests Merged
+
+#### Bugs fixed
+
+* [PR 64](https://github.com/pytroll/pyspectral/pull/64) - Fix interp function in rayleigh correction to be serializable
+
+#### Features added
+
+* [PR 59](https://github.com/pytroll/pyspectral/pull/59) - Daskify NIR reflectance calculations ([58](https://github.com/pytroll/pyspectral/issues/58))
+
+In this release 2 pull requests were closed.
+
+
 ## Version <RELEASE_VERSION> (2018/12/04)
 
 ### Issues Closed

README.md

@@ -7,7 +7,6 @@ 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)
-[![Research software impact](http://depsy.org/api/package/pypi/pyspectral/badge.svg)](http://depsy.org/package/python/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)
 

RELEASING.md

@@ -17,7 +17,7 @@ Don't forget to commit!
 5. Create a tag with the new version number, starting with a 'v', eg:
 
 ```
-git tag v0.8.4 -m "Version 0.8.4
+git tag v0.8.4 -m "Version 0.8.4"
 ```
 
 See [semver.org](http://semver.org/) on how to write a version number.

bin/composite_rsr_plot.py

 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 
-# Copyright (c) 2016-2018 Adam.Dybbroe
+# Copyright (c) 2016-2019 Adam.Dybbroe
 
 # Author(s):
 
@@ -74,6 +74,9 @@ def get_arguments():
     parser.add_argument("-x", "--xlimits", nargs=2,
                         help=("x-axis boundaries for plot"),
                         default=None, type=float)
+    parser.add_argument("-y", "--ylimits", nargs=2,
+                        help=("y-axis boundaries for plot"),
+                        default=None, type=float)
     parser.add_argument("-t", "--minimum_response",
                         help=("Minimum response: Any response lower than " +
                               "this will be ignored when plotting"),
@@ -115,6 +118,7 @@ if __name__ == "__main__":
     sensors = args.sensor
     minimum_response = args.minimum_response
     xlimits = args.xlimits
+    ylimits = args.ylimits
     title = args.title
     if not title:
         title = 'Relative Spectral Responses'
@@ -211,6 +215,13 @@ if __name__ == "__main__":
 
     plt.xlim((wmin, wmax))
 
+    wmin, wmax = plt.ylim()
+    if ylimits:
+        wmin = ylimits[0]
+        wmax = ylimits[1]
+
+    plt.ylim((wmin, wmax))
+
     plt.title(title)
     plt.legend(loc='lower right')
     if filename:

debian/changelog

+pyspectral (0.8.7+ds-1) unstable; urgency=medium
+
+  * New upstream release.
+  * Install files in pyspectral/etc.
+  * Updated debian/copyright.
+  * debian/patches:
+    - refresh all patches
+
+ -- Antonio Valentino <antonio.valentino@tiscali.it>  Mon, 22 Apr 2019 09:28:29 +0000
+
 pyspectral (0.8.6+ds-1) unstable; urgency=medium
 
   * Initial version (Closes: #917351)

debian/copyright

@@ -16,7 +16,7 @@ Files-Excluded: pyspectral/data/modis
  doc/_static/mersi2_rsr_band_0040_0070_missingbands.png
 
 Files: *
-Copyright: 2013-2018 Adam Dybbroe
+Copyright: 2013-2019 Adam Dybbroe
                      PyTroll Community
                      Pytroll
 License: GPL-3+

debian/patches/0001-Fix-tests.patch

@@ -8,7 +8,7 @@ Disable tests that require the network or external data
  1 file changed, 4 insertions(+), 4 deletions(-)
 
 diff --git a/pyspectral/tests/__init__.py b/pyspectral/tests/__init__.py
-index a4e8886..826d17b 100644
+index b930fc2..686ce89 100644
 --- a/pyspectral/tests/__init__.py
 +++ b/pyspectral/tests/__init__.py
 @@ -53,10 +53,10 @@ def suite():

debian/python3-pyspectral.install

 usr/lib/python3*/dist-packages/*/*.py
+usr/lib/python3*/dist-packages/*/etc/*
 usr/lib/python3*/dist-packages/*/data/*
 usr/lib/python3*/dist-packages/*.egg-info/*

doc/37_reflectance.rst

@@ -45,7 +45,7 @@ expressed in :math:`W/m^2 sr^{-1} \mu m^{-1}`, or using SI units :math:`W/m^2 sr
   >>> tb37 = np.array([298.07385254, 297.15478516, 294.43276978, 281.67633057, 273.7923584])
   >>> viirs = RadTbConverter('Suomi-NPP', 'viirs', 'M12')
   >>> rad37 = viirs.tb2radiance(tb37)
-  >>> print([np.round(rad, 7) for rad in rad37['radiance'].data])
+  >>> print([np.round(rad, 7) for rad in rad37['radiance']])
   [369717.4765726, 355110.5207853, 314684.2788726, 173143.5424898, 116408.0007877]
   >>> rad37['unit']
   'W/m^2 sr^-1 m^-1'
@@ -58,7 +58,7 @@ In order to get the total radiance over the band one has to multiply with the eq
   >>> tb37 = np.array([298.07385254, 297.15478516, 294.43276978, 281.67633057, 273.7923584])
   >>> viirs = RadTbConverter('Suomi-NPP', 'viirs', 'M12')
   >>> rad37 = viirs.tb2radiance(tb37, normalized=False)
-  >>> print([np.round(rad, 8) for rad in rad37['radiance'].data])
+  >>> print([np.round(rad, 8) for rad in rad37['radiance']])
   [0.07037968, 0.06759909, 0.05990352, 0.03295972, 0.02215951]
   >>> rad37['unit']
   'W/m^2 sr^-1'
@@ -127,7 +127,7 @@ where :math:`S(\lambda)` is the spectral solar irradiance.
   >>> viirs = RelativeSpectralResponse('Suomi-NPP', 'viirs')
   >>> solar_irr = SolarIrradianceSpectrum(TOTAL_IRRADIANCE_SPECTRUM_2000ASTM, dlambda=0.005)
   >>> sflux = solar_irr.inband_solarflux(viirs.rsr['M12'])
-  >>> print(round(sflux, 7))
+  >>> print(np.round(sflux, 7))
   2.2428119
 
 Derive the reflective part of the observed 3.7 micron radiance
@@ -197,10 +197,10 @@ In Python this becomes:
   >>> tb37 = np.array([298.07385254, 297.15478516, 294.43276978, 281.67633057, 273.7923584])
   >>> tb11 = np.array([271.38806152, 271.38806152, 271.33453369, 271.98553467, 271.93609619])
   >>> m12r = refl_m12.reflectance_from_tbs(sunz, tb37, tb11)
-  >>> print(m12r.mask)
+  >>> print(np.any(np.isnan(m12r)))
   False
-  >>> print([round(refl, 8) for refl in m12r.data])
-  [0.21432927, 0.20285153, 0.17063976, 0.05408903, 0.00838111]
+  >>> print([np.round(refl, 7) for refl in m12r])
+  [0.2143293, 0.2028515, 0.1706398, 0.054089, 0.0083811]
   
 We can try decompose equation :eq:`refl37` above using the example of VIIRS M12 band:
 
@@ -215,19 +215,19 @@ We can try decompose equation :eq:`refl37` above using the example of VIIRS M12
   >>> rad11 = viirs.tb2radiance(tb11, normalized=False)
   >>> sflux = 2.242817881698326
   >>> nomin = rad37['radiance'] - rad11['radiance']
-  >>> print(nomin.mask)
+  >>> print(np.isnan(nomin))
   [False False False False False]
-  >>> print([np.round(val, 8) for val in nomin.data])
+  >>> print([np.round(val, 8) for val in nomin])
   [0.05083677, 0.04805618, 0.0404157, 0.01279279, 0.00204485]
   >>> denom = np.cos(np.deg2rad(sunz))/np.pi * sflux - rad11['radiance']
-  >>> print(denom.mask)
+  >>> print(np.isnan(denom))
   [False False False False False]
-  >>> print([np.round(val, 8) for val in denom.data])
+  >>> print([np.round(val, 8) for val in denom])
   [0.23646313, 0.23645682, 0.23650559, 0.23582015, 0.2358661]
   >>> res = nomin/denom
-  >>> print(res.mask)
+  >>> print(np.isnan(res))
   [False False False False False]
-  >>> print([round(val, 8) for val in res.data])
+  >>> print([np.round(val, 8) for val in res])
   [0.21498817, 0.2032345, 0.17088689, 0.05424807, 0.00866955]
 
 
@@ -252,8 +252,8 @@ Using the example of the VIIRS M12 band from above this gives the following spec
   >>> tb11 = np.array([271.38806152, 271.38806152, 271.33453369, 271.98553467, 271.93609619])
   >>> m12r = refl_m12.reflectance_from_tbs(sunz, tb37, tb11)
   >>> tb = refl_m12.emissive_part_3x()
-  >>> ['{tb:6.3f}'.format(tb=round(t, 4)) for t in tb.data]
+  >>> ['{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=round(r, 4)) for r in rad.data]
-  ['80285.149', '81458.022', '84749.639', '99761.400', '104582.030']
+  >>> ['{rad:6.3f}'.format(rad=np.round(r, 4)) for r in rad]
+  ['80285.149', '81458.022', '84749.639', '99761.401', '104582.030']

doc/platforms_supported.rst

@@ -76,7 +76,12 @@ have been included in PySpectral.
     * - FY-3D mersi-2
       - `rsr_mersi-2_FY-3D.h5`
       - CMA_ (Acquired via personal contact)
-
+    * - HY-1C cocts
+      - `rsr_cocts_HY-1C.h5`
+      - (Acquired via personal contact)
+    * - Metop-SG-A1 MetImage
+      - `rsr_metimage_Metop-SG-A1.h5`
+      - NWPSAF-MetImage_
 
 .. _Eumetsat: https://www.eumetsat.int/website/home/Data/Products/Calibration/MSGCalibration/index.html
 .. _GSICS: https://www.star.nesdis.noaa.gov/smcd/GCC/instrInfo-srf.php
@@ -89,4 +94,5 @@ have been included in PySpectral.
 .. _NASA-Landsat-OLI: https://landsat.gsfc.nasa.gov/wp-content/uploads/2013/06/Ball_BA_RSR.v1.1-1.xlsx
 .. _NESDIS: https://ncc.nesdis.noaa.gov/J1VIIRS/J1VIIRSSpectralResponseFunctions.php
 .. _CMA: http://www.cma.gov.cn/en2014/
+.. _NWPSAF-MetImage: https://nwpsaf.eu/downloads/rtcoef_rttov12/ir_srf/rtcoef_metopsg_1_metimage_srf.html
 

doc/usage.rst

@@ -22,9 +22,11 @@ Now, you can work with the data as you wish, make some simple plot for instance:
   >>> import matplotlib.pyplot as plt
   >>> dummy = plt.figure(figsize=(10, 5))
   >>> import numpy as np
-  >>> resp = np.ma.masked_less_equal(olci.rsr['Oa01']['det-1']['response'], 0.015)
-  >>> wvl = np.ma.masked_array(olci.rsr['Oa01']['det-1']['wavelength'], resp.mask)
-  >>> dummy = plt.plot(wvl.compressed(), resp.compressed())
+  >>> rsr = olci.rsr['Oa01']['det-1']['response']
+  >>> resp = np.where(rsr < 0.015, np.nan, rsr)
+  >>> wl_ = olci.rsr['Oa01']['det-1']['wavelength']
+  >>> wvl = np.where(np.isnan(resp), np.nan, wl_)
+  >>> dummy = plt.plot(wvl, resp)
   >>> plt.show() # doctest: +SKIP
 
   .. image:: _static/olci_ch1.png

pyspectral/blackbody.py

@@ -182,6 +182,7 @@ 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:

pyspectral/near_infrared_reflectance.py

@@ -28,6 +28,11 @@ window channel (usually around 11-12 microns).
 
 import os
 import numpy as np
+try:
+    from dask.array import where, logical_or
+except ImportError:
+    from numpy import where, logical_or
+
 from pyspectral.solar import (SolarIrradianceSpectrum,
                               TOTAL_IRRADIANCE_SPECTRUM_2000ASTM)
 from pyspectral.utils import BANDNAMES, get_bandname_from_wavelength
@@ -195,15 +200,25 @@ class Calculator(RadTbConverter):
                      absorption correction will be applied.
 
         """
+        # Check for dask arrays
+        if hasattr(tb_near_ir, 'compute') or hasattr(tb_thermal, 'compute'):
+            compute = False
+        else:
+            compute = True
+        if hasattr(tb_near_ir, 'mask') or hasattr(tb_thermal, 'mask'):
+            is_masked = True
+        else:
+            is_masked = False
+
         if np.isscalar(tb_near_ir):
             tb_nir = np.array([tb_near_ir, ])
         else:
-            tb_nir = np.array(tb_near_ir)
+            tb_nir = np.asanyarray(tb_near_ir)
 
         if np.isscalar(tb_thermal):
             tb_therm = np.array([tb_thermal, ])
         else:
-            tb_therm = np.array(tb_thermal)
+            tb_therm = np.asanyarray(tb_thermal)
 
         if tb_therm.shape != tb_nir.shape:
             errmsg = 'Dimensions do not match! {0} and {1}'.format(
@@ -221,7 +236,7 @@ class Calculator(RadTbConverter):
             if np.isscalar(tb_ir_co2):
                 tbco2 = np.array([tb_ir_co2, ])
             else:
-                tbco2 = np.array(tb_ir_co2)
+                tbco2 = np.asanyarray(tb_ir_co2)
 
         if not self.rsr:
             raise NotImplementedError("Reflectance calculations without "
@@ -239,9 +254,8 @@ class Calculator(RadTbConverter):
         if l_nir.ravel().shape[0] < 10:
             LOG.info('l_nir = %s', str(l_nir))
 
-        sunz = np.ma.masked_outside(sun_zenith, 0.0, 88.0)
-        sunzmask = sunz.mask
-        sunz = sunz.filled(88.)
+        sunzmask = (sun_zenith < 0.0) | (sun_zenith > 88.0)
+        sunz = where(sunzmask, 88.0, sun_zenith)
 
         mu0 = np.cos(np.deg2rad(sunz))
         # mu0 = np.where(np.less(mu0, 0.1), 0.1, mu0)
@@ -262,12 +276,18 @@ class Calculator(RadTbConverter):
         mask = (self._solar_radiance - thermal_emiss_one *
                 self._rad3x_correction) < EPSILON
 
-        np.logical_or(sunzmask, mask, out=mask)
-        np.logical_or(mask, np.isnan(tb_nir), out=mask)
+        logical_or(sunzmask, mask, out=mask)
+        logical_or(mask, np.isnan(tb_nir), out=mask)
 
-        self._r3x = np.ma.masked_where(mask, data)
+        self._r3x = where(mask, np.nan, data)
 
         # Reflectances should be between 0 and 1, but values above 1 is
         # perfectly possible and okay! (Multiply by 100 to get reflectances
         # in percent)
-        return self._r3x
+        if hasattr(self._r3x, 'compute') and compute:
+            res = self._r3x.compute()
+        else:
+            res = self._r3x
+        if is_masked:
+            res = np.ma.masked_array(res, mask=np.isnan(res))
+        return res

pyspectral/radiance_tb_conversion.py

@@ -228,7 +228,7 @@ class RadTbConverter(object):
         tb_ = np.arange(TB_MIN, TB_MAX, self.tb_resolution)
         retv = self.tb2radiance(tb_, normalized=normalized)
         rad = retv['radiance']
-        np.savez(filepath, tb=tb_, radiance=rad.compressed())
+        np.savez(filepath, tb=tb_, radiance=rad)
 
     @staticmethod
     def read_tb2rad_lut(filepath):

pyspectral/raw_reader.py

 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 
-# Copyright (c) 2016-2017 Adam.Dybbroe
+# Copyright (c) 2016-2019 Adam.Dybbroe
 
 # Author(s):
 

pyspectral/rayleigh.py

@@ -204,6 +204,12 @@ class Rayleigh(object):
         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 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:
@@ -277,19 +283,11 @@ class Rayleigh(object):
         minterp = MultilinearInterpolator(smin, smax, orders)
         minterp.set_values(f_3d_grid)
 
-        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)
-
         if HAVE_DASK:
-            ipn = map_blocks(_do_interp, minterp, sunzsec, azidiff,
-                             satzsec, dtype=raylwvl.dtype,
-                             chunks=azidiff.chunks)
+            ipn = map_blocks(self._do_interp, minterp, sunzsec, azidiff,
+                             satzsec, dtype=raylwvl.dtype, chunks=azidiff.chunks)
         else:
-            ipn = _do_interp(minterp, sunzsec, azidiff, satzsec)
+            ipn = self._do_interp(minterp, sunzsec, azidiff, satzsec)
 
         LOG.debug("Time - Interpolation: {0:f}".format(time.time() - tic))
 

pyspectral/tests/__init__.py

 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 
-# Copyright (c) 2013-2018 Adam.Dybbroe
+# Copyright (c) 2013-2019 Adam.Dybbroe
 
 # Author(s):
 

pyspectral/tests/test_reflectance.py

@@ -120,7 +120,7 @@ class TestReflectance(unittest.TestCase):
             refl37 = Calculator('EOS-Aqua', 'modis', '20')
 
         expected = 1.8563451e-07  # unit = 'm' (meter)
-        self.assertAlmostEqual(refl37.rsr_integral, expected)
+        np.testing.assert_allclose(refl37.rsr_integral, expected)
 
         with patch('pyspectral.radiance_tb_conversion.RelativeSpectralResponse') as mymock:
             instance = mymock.return_value
@@ -132,7 +132,7 @@ class TestReflectance(unittest.TestCase):
 
         expected = 13000.385  # SI units = 'm-1' (1/meter)
         res = refl37.rsr_integral
-        self.assertAlmostEqual(res / expected, 1.0, 6)
+        np.testing.assert_allclose(res / expected, 1.0, 6)
 
     def test_reflectance(self):
         """Test the derivation of the reflective part of a 3.7 micron band"""
@@ -162,28 +162,44 @@ class TestReflectance(unittest.TestCase):
         tb3 = np.array([290.])
         tb4 = np.array([282.])
         refl = refl37.reflectance_from_tbs(sunz, tb3, tb4)
-        self.assertAlmostEqual(refl.data[0], 0.251245010648, 6)
+        np.testing.assert_allclose(refl[0], 0.251245010648, 6)
 
         tb3x = refl37.emissive_part_3x()
-        self.assertAlmostEqual(tb3x, 276.213054, 6)
+        np.testing.assert_allclose(tb3x, 276.213054, 6)
 
         sunz = np.array([80.])
         tb3 = np.array([295.])
         tb4 = np.array([282.])
         refl = refl37.reflectance_from_tbs(sunz, tb3, tb4)
-        self.assertAlmostEqual(refl.data[0], 0.452497961, 6)
+        np.testing.assert_allclose(refl[0], 0.452497961, 6)
 
         tb3x = refl37.emissive_part_3x()
-        self.assertAlmostEqual(tb3x, 270.077268, 6)
+        np.testing.assert_allclose(tb3x, 270.077268, 6)
 
         sunz = np.array([50.])
         tb3 = np.array([300.])
         tb4 = np.array([285.])
         refl = refl37.reflectance_from_tbs(sunz, tb3, tb4)
-        self.assertAlmostEqual(refl.data[0], 0.1189217, 6)
+        np.testing.assert_allclose(refl[0], 0.1189217, 6)
 
         tb3x = refl37.emissive_part_3x()
-        self.assertAlmostEqual(tb3x, 282.455426, 6)
+        np.testing.assert_allclose(tb3x, 282.455426, 6)
+
+        sunz = np.array([50.])
+        tb3 = np.ma.masked_array([300.], mask=False)
+        tb4 = np.ma.masked_array([285.], mask=False)
+        refl = refl37.reflectance_from_tbs(sunz, tb3, tb4)
+        self.assertTrue(hasattr(refl, 'mask'))
+
+        try:
+            import dask.array as da
+            sunz = da.from_array([50.], chunks=10)
+            tb3 = da.from_array([300.], chunks=10)
+            tb4 = da.from_array([285.], chunks=10)
+            refl = refl37.reflectance_from_tbs(sunz, tb3, tb4)
+            self.assertTrue(hasattr(refl, 'compute'))
+        except ImportError:
+            pass
 
     def tearDown(self):
         """Clean up"""

pyspectral/utils.py

 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 
-# Copyright (c) 2014-2018 Pytroll
+# Copyright (c) 2014-2019 Pytroll
 
 # Author(s):
 
@@ -123,9 +123,9 @@ INSTRUMENTS = {'NOAA-19': 'avhrr/3',
                'Feng-Yun 3D': 'mersi-2'
                }
 
-HTTP_PYSPECTRAL_RSR = "https://zenodo.org/record/1491277/files/pyspectral_rsr_data.tgz"
+HTTP_PYSPECTRAL_RSR = "https://zenodo.org/record/2617441/files/pyspectral_rsr_data.tgz"
 RSR_DATA_VERSION_FILENAME = "PYSPECTRAL_RSR_VERSION"
-RSR_DATA_VERSION = "v1.0.3"
+RSR_DATA_VERSION = "v1.0.5"
 
 ATM_CORRECTION_LUT_VERSION = {}
 ATM_CORRECTION_LUT_VERSION['antarctic_aerosol'] = {'version': 'v1.0.1',