Merge pull request #59 from pnuu/feature-daskify

Daskify NIR reflectance calculations

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,9 +197,9 @@ 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])
+  >>> print([np.round(refl, 8) for refl in m12r])
   [0.21432927, 0.20285153, 0.17063976, 0.05408903, 0.00838111]
 
 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]
+  >>> ['{rad:6.3f}'.format(rad=np.round(r, 4)) for r in rad]
   ['80285.149', '81458.022', '84749.639', '99761.400', '104582.030']

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 @@
 
 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 @@ def reflectance_from_tbs(self, sun_zenith, tb_near_ir, tb_thermal, **kwargs):
                      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 @@ def reflectance_from_tbs(self, sun_zenith, tb_near_ir, tb_thermal, **kwargs):
             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 @@ def reflectance_from_tbs(self, sun_zenith, tb_near_ir, tb_thermal, **kwargs):
         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 @@ def reflectance_from_tbs(self, sun_zenith, tb_near_ir, tb_thermal, **kwargs):
         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 @@ def make_tb2rad_lut(self, filepath, normalized=True):
         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/tests/test_reflectance.py

@@ -120,7 +120,7 @@ def test_rsr_integral(self):
             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 @@ def test_rsr_integral(self):
 
         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 @@ def test_reflectance(self):
         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"""