Daskify NIR reflectance calculations

pyspectral/blackbody.py

@@ -23,6 +23,7 @@
 
 """Planck radiation equation"""
 import numpy as np
+import dask.array as da
 
 import logging
 LOG = logging.getLogger(__name__)
@@ -41,25 +42,21 @@ def blackbody_rad2temp(wavelength, radiance):
 
     """
 
-    mask = False
     if np.isscalar(radiance):
-        rad = np.array([radiance, ], dtype='float64')
+        rad = da.array([radiance, ], dtype='float64')
     else:
-        rad = np.array(radiance, dtype='float64')
-        if np.ma.is_masked(radiance):
-            mask = radiance.mask
+        rad = da.array(radiance, dtype='float64')
 
-    rad = np.ma.masked_array(rad, mask=mask)
-    rad = np.ma.masked_less_equal(rad, 0)
+    rad = da.where(rad < 0, np.nan, rad)
 
     if np.isscalar(wavelength):
-        wvl = np.array([wavelength, ], dtype='float64')
+        wvl = da.array([wavelength, ], dtype='float64')
     else:
-        wvl = np.array(wavelength, dtype='float64')
+        wvl = da.array(wavelength, dtype='float64')
 
     const1 = H_PLANCK * C_SPEED / K_BOLTZMANN
     const2 = 2 * H_PLANCK * C_SPEED**2
-    res = const1 / (wvl * np.log(np.divide(const2, rad * wvl**5) + 1.0))
+    res = const1 / (wvl * da.log(da.divide(const2, rad * wvl**5) + 1.0))
 
     shape = rad.shape
     resshape = res.shape
@@ -68,33 +65,33 @@ def blackbody_rad2temp(wavelength, radiance):
         if rad.shape[0] == 1:
             return res[0]
         else:
-            return res[::].reshape(shape)
+            return da.reshape(res[::], shape)
     else:
         if rad.shape[0] == 1:
             return res[0, :]
         else:
             if len(shape) == 1:
-                return np.reshape(res, (shape[0], resshape[1]))
+                return da.reshape(res, (shape[0], resshape[1]))
             else:
-                return np.reshape(res, (shape[0], shape[1], resshape[1]))
+                return da.reshape(res, (shape[0], shape[1], resshape[1]))
 
 
 def blackbody_wn_rad2temp(wavenumber, radiance):
     """Derive brightness temperatures from radiance using the Planck 
     function. Wavenumber space"""
 
     if np.isscalar(radiance):
-        rad = np.array([radiance, ], dtype='float64')
+        rad = da.array([radiance, ], dtype='float64')
     else:
-        rad = np.array(radiance, dtype='float64')
+        rad = da.array(radiance, dtype='float64')
     if np.isscalar(wavenumber):
-        wavnum = np.array([wavenumber, ], dtype='float64')
+        wavnum = da.array([wavenumber, ], dtype='float64')
     else:
-        wavnum = np.array(wavenumber, dtype='float64')
+        wavnum = da.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 / da.log(da.divide(const2 * wavnum**3, rad) + 1.0)
 
     shape = rad.shape
     resshape = res.shape
@@ -103,15 +100,15 @@ def blackbody_wn_rad2temp(wavenumber, radiance):
         if rad.shape[0] == 1:
             return res[0]
         else:
-            return res[::].reshape(shape)
+            return da.reshape(res[::], shape)
     else:
         if rad.shape[0] == 1:
             return res[0, :]
         else:
             if len(shape) == 1:
-                return np.reshape(res, (shape[0], resshape[1]))
+                return da.reshape(res, (shape[0], resshape[1]))
             else:
-                return np.reshape(res, (shape[0], shape[1], resshape[1]))
+                return da.reshape(res, (shape[0], shape[1], resshape[1]))
 
 
 def planck(wave, temp, wavelength=True):
@@ -139,15 +136,15 @@ def planck(wave, temp, wavelength=True):
             units[(wavelength == True) - 1]))
 
     if np.isscalar(temp):
-        temperature = np.array([temp, ], dtype='float64')
+        temperature = da.array([temp, ], dtype='float64')
     else:
-        temperature = np.array(temp, dtype='float64')
+        temperature = da.array(temp, dtype='float64')
 
     shape = temperature.shape
     if np.isscalar(wave):
-        wln = np.array([wave, ], dtype='float64')
+        wln = da.array([wave, ], dtype='float64')
     else:
-        wln = np.array(wave, dtype='float64')
+        wln = da.array(wave, dtype='float64')
 
     if wavelength:
         const = 2 * H_PLANCK * C_SPEED ** 2
@@ -157,30 +154,31 @@ 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)
+    arg2 = da.reshape(da.where(da.greater(da.absolute(temperature), EPSILON),
+                               da.array(1. / temperature), np.nan),
+                      (-1, 1))
     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()))
+    # LOG.debug("Max and min - arg2: %s  %s", str(arg2.max()), str(arg2.min()))
     try:
-        exp_arg = np.multiply(arg1.astype('float32'), arg2.astype('float32'))
+        # exp_arg = da.multiply(arg1.astype('float32'), arg2.astype('float32'))
+        exp_arg = da.multiply(arg1, arg2)
     except MemoryError:
         LOG.warning(("Dimensions used in numpy.multiply probably reached "
                      "limit!\n"
                      "Make sure the Radiance<->Tb table has been created "
                      "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:
+    # LOG.debug("Max and min before exp: %s  %s", str(exp_arg.max()),
+    #           str(exp_arg.min()))
+    if da.min(exp_arg) < 0:
         LOG.warning("Something is fishy: \n" +
                     "\tDenominator might be zero or negative in radiance derivation:")
-        dubious = np.where(exp_arg < 0)[0]
+        dubious = da.where(exp_arg < 0)[0]
         LOG.warning(
-            "Number of items having dubious values: " + str(dubious.shape[0]))
+            "Number of items having dubious values: %d", dubious.shape[0])
 
-    denom = np.exp(exp_arg) - 1
+    denom = da.exp(exp_arg) - 1
     rad = nom / denom
     radshape = rad.shape
     if wln.shape[0] == 1:
@@ -193,9 +191,9 @@ def planck(wave, temp, wavelength=True):
             return rad[0, :]
         else:
             if len(shape) == 1:
-                return np.reshape(rad, (shape[0], radshape[1]))
+                return da.reshape(rad, (shape[0], radshape[1]))
             else:
-                return np.reshape(rad, (shape[0], shape[1], radshape[1]))
+                return da.reshape(rad, (shape[0], shape[1], radshape[1]))
 
 
 def blackbody_wn(wavenumber, temp):

pyspectral/near_infrared_reflectance.py

@@ -28,6 +28,7 @@
 
 import os
 import numpy as np
+import dask.array as da
 from pyspectral.solar import (SolarIrradianceSpectrum,
                               TOTAL_IRRADIANCE_SPECTRUM_2000ASTM)
 from pyspectral.utils import BANDNAMES, get_bandname_from_wavelength
@@ -239,9 +240,11 @@ 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 = da.where(sunzmask, 88.0, sun_zenith)
+        # sunz = np.ma.masked_outside(sun_zenith, 0.0, 88.0)
+        # sunzmask = sunz.mask
+        # sunz = sunz.filled(88.)
 
         mu0 = np.cos(np.deg2rad(sunz))
         # mu0 = np.where(np.less(mu0, 0.1), 0.1, mu0)
@@ -262,10 +265,10 @@ 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)
+        da.logical_or(sunzmask, mask, out=mask)
+        da.logical_or(mask, np.isnan(tb_nir), out=mask)
 
-        self._r3x = np.ma.masked_where(mask, data)
+        self._r3x = da.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

pyspectral/radiance_tb_conversion.py

@@ -26,6 +26,7 @@
 """
 
 import numpy as np
+import dask.array as da
 from pyspectral.blackbody import blackbody, blackbody_wn
 from pyspectral.utils import BANDNAMES
 from pyspectral.utils import get_bandname_from_wavelength
@@ -228,7 +229,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)  #.compressed())
 
     @staticmethod
     def read_tb2rad_lut(filepath):
@@ -308,7 +309,7 @@ def radiance2tb(self, rad):
         beta = SEVIRI[self.bandname][self.platform_name][2]
 
         tb_ = c_2 * vc_ / \
-            (alpha * np.log(c_1 * vc_ ** 3 / rad + 1)) - beta / alpha
+            (alpha * da.log(c_1 * vc_ ** 3 / rad + 1)) - beta / alpha
 
         return tb_
 
@@ -338,7 +339,7 @@ def tb2radiance(self, tb_, **kwargs):
         beta = SEVIRI[self.bandname][self.platform_name][2]
 
         radiance = c_1 * vc_ ** 3 / \
-            (np.exp(c_2 * vc_ / (alpha * tb_ + beta)) - 1)
+            (da.exp(c_2 * vc_ / (alpha * tb_ + beta)) - 1)
 
         unit = 'W/m^2 sr^-1 (m^-1)^-1'
         scale = 1.0