Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Daskify NIR reflectance calculations #59

Merged
merged 13 commits into from
Apr 9, 2019
Merged

Daskify NIR reflectance calculations #59

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
d94cb88
Use dask instead of Numpy where suitable
pnuu Nov 28, 2018
4b88bdb
Use np.testing.assert_allclose instead of assertAlmostEqual
pnuu Nov 29, 2018
a695588
Revert back to use numpy operations
pnuu Nov 29, 2018
71ff135
Use np.asanyarray instead of np.array
pnuu Nov 29, 2018
45e38b1
Remove references to masked array attributes
pnuu Nov 29, 2018
654e3e7
Change documentation to use plain Numpy arrays instead of masked arrays
pnuu Nov 29, 2018
01933a7
Remove obsolete commented-out lines
pnuu Nov 29, 2018
b977c21
Don't rely on dask to be available
pnuu Nov 29, 2018
92b743c
Compute the results if the input is not dask array
pnuu Nov 29, 2018
b22d176
Fix da.log and da.exp to np.log and da.log
pnuu Nov 30, 2018
500dc7f
Convert the result to masked array if also the input is masked
pnuu Nov 30, 2018
781e477
Add tests for checking output for masked and dask arrays stay the same
pnuu Nov 30, 2018
86086df
Ensure that radiance from planck() isn't masked
pnuu Nov 30, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
26 changes: 13 additions & 13 deletions doc/37_reflectance.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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'
Expand All @@ -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'
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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:
Expand All @@ -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]


Expand All @@ -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']
8 changes: 5 additions & 3 deletions doc/usage.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
1 change: 1 addition & 0 deletions pyspectral/blackbody.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down
40 changes: 30 additions & 10 deletions pyspectral/near_infrared_reflectance.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does np.asanyarray convert dask arrays to numpy arrays?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It does:

In [1]: import dask.array as da                                                                                                                                      

In [2]: a = da.arange(5., chunks=2)                                                                                                                                  

In [3]: import numpy as np                                                                                                                                           

In [4]: type(np.asanyarray(a))                                                                                                                                       
Out[4]: numpy.ndarray

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@djhoese what do you suggest?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As mentioned on slack, I'd prefer if everything could be kept as dask arrays for as long as possible and if things have to be computed they should be done in a map_blocks or delayed function. As @pnuu mentioned on slack converting to numpy here made the code faster and use less memory in his test cases. I'm fairly certain this is because the array is used later on in a non-dask-friendly way and the dask array is actually being computed multiple times. I think I pointed them out on slack during the PCW.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As this part of pyspectral wasn't dask friendly before, could we perhaps merge this PR now and make sure things works and then work further to improve the daskfriendliness?
@mraspaud @djhoese @pnuu ?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm fine with merging. Then we could also merge pytroll/satpy#529 so this PR would be usable via SatpY.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So this is working now ?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think there hasn't been anything that didn't work since I last touched this, just that it might not be quite optimal in performance and how things should be done with dask.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fine. +0

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fine, I will just try do some system testing off line first to be sure. Any suggested recipes/tests @pnuu ?


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(
Expand All @@ -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 "
Expand All @@ -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)
Expand All @@ -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
2 changes: 1 addition & 1 deletion pyspectral/radiance_tb_conversion.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down
32 changes: 24 additions & 8 deletions pyspectral/tests/test_reflectance.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -116,7 +116,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
Expand All @@ -128,7 +128,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"""
Expand Down Expand Up @@ -158,28 +158,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"""
Expand Down