Add configurable parameters to solar zenith correctors

satpy/composites/__init__.py

@@ -313,17 +313,26 @@ class SunZenithCorrectorBase(CompositeBase):
 
     coszen = WeakValueDictionary()
 
+    def __init__(self, max_sza=95.0, **kwargs):
+        """Collect custom configuration values.
+
+        Args:
+            max_sza (float): Maximum solar zenith angle in degrees that is
+                considered valid and correctable. Default 95.0.
+
+        """
+        self.max_sza = max_sza
+        self.max_sza_cos = np.cos(np.deg2rad(max_sza)) if max_sza is not None else None
+        super(SunZenithCorrectorBase, self).__init__(**kwargs)
+
     def __call__(self, projectables, **info):
         projectables = self.check_areas(projectables)
         vis = projectables[0]
         if vis.attrs.get("sunz_corrected"):
             LOG.debug("Sun zen correction already applied")
             return vis
 
-        if hasattr(vis.attrs["area"], 'name'):
-            area_name = vis.attrs["area"].name
-        else:
-            area_name = 'swath' + str(vis.shape)
+        area_name = hash(vis.attrs['area'])
         key = (vis.attrs["start_time"], area_name)
         tic = time.time()
         LOG.debug("Applying sun zen correction")
@@ -334,46 +343,108 @@ def __call__(self, projectables, **info):
                 LOG.debug("Computing sun zenith angles.")
                 lons, lats = vis.attrs["area"].get_lonlats_dask(CHUNK_SIZE)
 
-                coszen = xr.DataArray(cos_zen(vis.attrs["start_time"],
-                                              lons, lats),
-                                      dims=['y', 'x'],
-                                      coords=[vis['y'], vis['x']])
-                coszen = coszen.where((coszen > 0.035) & (coszen < 1))
+                coszen = xr.DataArray(cos_zen(vis.attrs["start_time"], lons, lats),
+                                      dims=['y', 'x'], coords=[vis['y'], vis['x']])
+                if self.max_sza is not None:
+                    coszen = coszen.where(coszen >= self.max_sza_cos)
                 self.coszen[key] = coszen
         else:
-            coszen = xu.cos(xu.deg2rad(projectables[1]))
+            coszen = np.cos(np.deg2rad(projectables[1]))
             self.coszen[key] = coszen
 
         proj = self._apply_correction(vis, coszen)
         proj.attrs = vis.attrs.copy()
         self.apply_modifier_info(vis, proj)
-        LOG.debug(
-            "Sun-zenith correction applied. Computation time: %5.1f (sec)",
-            time.time() - tic)
+        LOG.debug("Sun-zenith correction applied. Computation time: %5.1f (sec)", time.time() - tic)
         return proj
 
     def _apply_correction(self, proj, coszen):
         raise NotImplementedError("Correction method shall be defined!")
 
 
 class SunZenithCorrector(SunZenithCorrectorBase):
-    """Standard sun zenith correction, 1/cos(sunz)."""
+    """Standard sun zenith correction using ``1 / cos(sunz)``.
+
+    In addition to adjusting the provided reflectances by the cosine of the
+    solar zenith angle, this modifier forces all reflectances beyond a
+    solar zenith angle of ``max_sza`` to 0. It also gradually reduces the
+    amount of correction done between ``correction_limit`` and ``max_sza``. If
+    ``max_sza`` is ``None`` then a constant correction is applied to zenith
+    angles beyond ``correction_limit``.
+
+    To set ``max_sza`` to ``None`` in a YAML configuration file use:
+
+    .. code-block:: yaml
+
+      sunz_corrected:
+        compositor: !!python/name:satpy.composites.SunZenithCorrector
+        max_sza: !!null
+        optional_prerequisites:
+        - solar_zenith_angle
+
+    """
+
+    def __init__(self, correction_limit=88., **kwargs):
+        """Collect custom configuration values.
+
+        Args:
+            correction_limit (float): Maximum solar zenith angle to apply the
+                correction in degrees. Pixels beyond this limit have a
+                constant correction applied. Default 88.
+            max_sza (float): Maximum solar zenith angle in degrees that is
+                considered valid and correctable. Default 95.0.
+
+        """
+        self.correction_limit = correction_limit
+        super(SunZenithCorrector, self).__init__(**kwargs)
 
     def _apply_correction(self, proj, coszen):
         LOG.debug("Apply the standard sun-zenith correction [1/cos(sunz)]")
-        return sunzen_corr_cos(proj, coszen)
+        return sunzen_corr_cos(proj, coszen, limit=self.correction_limit, max_sza=self.max_sza)
 
 
 class EffectiveSolarPathLengthCorrector(SunZenithCorrectorBase):
     """Special sun zenith correction with the method proposed by Li and Shibata.
 
     (2006): https://doi.org/10.1175/JAS3682.1
 
+    In addition to adjusting the provided reflectances by the cosine of the
+    solar zenith angle, this modifier forces all reflectances beyond a
+    solar zenith angle of `max_sza` to 0 to reduce noise in the final data.
+    It also gradually reduces the amount of correction done between
+    ``correction_limit`` and ``max_sza``. If ``max_sza`` is ``None`` then a
+    constant correction is applied to zenith angles beyond
+    ``correction_limit``.
+
+    To set ``max_sza`` to ``None`` in a YAML configuration file use:
+
+    .. code-block:: yaml
+
+      effective_solar_pathlength_corrected:
+        compositor: !!python/name:satpy.composites.EffectiveSolarPathLengthCorrector
+        max_sza: !!null
+        optional_prerequisites:
+        - solar_zenith_angle
+
     """
 
+    def __init__(self, correction_limit=88., **kwargs):
+        """Collect custom configuration values.
+
+        Args:
+            correction_limit (float): Maximum solar zenith angle to apply the
+                correction in degrees. Pixels beyond this limit have a
+                constant correction applied. Default 88.
+            max_sza (float): Maximum solar zenith angle in degrees that is
+                considered valid and correctable. Default 95.0.
+
+        """
+        self.correction_limit = correction_limit
+        super(EffectiveSolarPathLengthCorrector, self).__init__(**kwargs)
+
     def _apply_correction(self, proj, coszen):
         LOG.debug("Apply the effective solar atmospheric path length correction method by Li and Shibata")
-        return atmospheric_path_length_correction(proj, coszen)
+        return atmospheric_path_length_correction(proj, coszen, limit=self.correction_limit, max_sza=self.max_sza)
 
 
 class PSPRayleighReflectance(CompositeBase):
@@ -829,7 +900,7 @@ def __call__(self, projectables, **kwargs):
         lim_low = np.cos(np.deg2rad(self.lim_low))
         lim_high = np.cos(np.deg2rad(self.lim_high))
         try:
-            coszen = xu.cos(xu.deg2rad(projectables[2]))
+            coszen = np.cos(np.deg2rad(projectables[2]))
         except IndexError:
             from pyorbital.astronomy import cos_zen
             LOG.debug("Computing sun zenith angles.")
@@ -1163,7 +1234,7 @@ def __call__(self, datasets, optional_datasets=None, **info):
             r, g, b = datasets[:3]
 
         # combine the masks
-        mask = ~(da.isnull(r.data) | da.isnull(g.data) | da.isnull(b.data))
+        mask = ~(r.isnull() | g.isnull() | b.isnull())
         r = r.where(mask)
         g = g.where(mask)
         b = b.where(mask)
@@ -1202,8 +1273,6 @@ def four_element_average_dask(d):
         """Average every 4 elements (2x2) in a 2D array"""
         def _mean4(data, offset=(0, 0), block_id=None):
             rows, cols = data.shape
-            rows2, cols2 = data.shape
-            pad = []
             # we assume that the chunks except the first ones are aligned
             if block_id[0] == 0:
                 row_offset = offset[0] % 2

satpy/composites/viirs.py

@@ -160,7 +160,7 @@ def get_angles(self, vis):
         lons, lats = vis.attrs['area'].get_lonlats_dask(
             chunks=vis.data.chunks)
         suna = get_alt_az(vis.attrs['start_time'], lons, lats)[1]
-        suna = xu.rad2deg(suna)
+        suna = np.rad2deg(suna)
         sunz = sun_zenith_angle(vis.attrs['start_time'], lons, lats)
         sata, satel = get_observer_look(
             vis.attrs['satellite_longitude'],

satpy/etc/composites/visir.yaml

@@ -4,10 +4,12 @@ modifiers:
   sunz_corrected:
     compositor: !!python/name:satpy.composites.SunZenithCorrector
     optional_prerequisites:
-    - satellite_zenith_angle
+    - solar_zenith_angle
 
   effective_solar_pathlength_corrected:
     compositor: !!python/name:satpy.composites.EffectiveSolarPathLengthCorrector
+    optional_prerequisites:
+      - solar_zenith_angle
 
   co2_corrected:
     compositor: !!python/name:satpy.composites.CO2Corrector

satpy/tests/compositor_tests/__init__.py

@@ -223,6 +223,58 @@ def test_self_sharpened_basic(self):
         np.testing.assert_allclose(res[2], np.array([[4, 4], [4, 4]], dtype=np.float64))
 
 
+class TestSunZenithCorrector(unittest.TestCase):
+    def setUp(self):
+        """Create test data."""
+        from pyresample.geometry import AreaDefinition
+        area = AreaDefinition('test', 'test', 'test',
+                              {'proj': 'merc'}, 2, 2,
+                              (-2000, -2000, 2000, 2000))
+        attrs = {'area': area,
+                 'start_time': datetime(2018, 1, 1, 18),
+                 'modifiers': tuple(),
+                 'name': 'test_vis'}
+        ds1 = xr.DataArray(da.ones((2, 2), chunks=2, dtype=np.float64),
+                           attrs=attrs, dims=('y', 'x'),
+                           coords={'y': [0, 1], 'x': [0, 1]})
+        self.ds1 = ds1
+        self.sza = xr.DataArray(
+            np.rad2deg(np.arccos(da.from_array([[0.0149581333, 0.0146694376], [0.0150812684, 0.0147925727]],
+                                               chunks=2))),
+            attrs={'area': area},
+            dims=('y', 'x'),
+            coords={'y': [0, 1], 'x': [0, 1]},
+        )
+
+    def test_basic_default_not_provided(self):
+        """Test default limits when SZA isn't provided."""
+        from satpy.composites import SunZenithCorrector
+        comp = SunZenithCorrector(name='sza_test', modifiers=tuple())
+        res = comp((self.ds1,), test_attr='test')
+        np.testing.assert_allclose(res.values, np.array([[22.401667, 22.31777], [22.437503, 22.353533]]))
+
+    def test_basic_lims_not_provided(self):
+        """Test custom limits when SZA isn't provided."""
+        from satpy.composites import SunZenithCorrector
+        comp = SunZenithCorrector(name='sza_test', modifiers=tuple(), correction_limit=90)
+        res = comp((self.ds1,), test_attr='test')
+        np.testing.assert_allclose(res.values, np.array([[66.853262, 68.168939], [66.30742, 67.601493]]))
+
+    def test_basic_default_provided(self):
+        """Test default limits when SZA is provided."""
+        from satpy.composites import SunZenithCorrector
+        comp = SunZenithCorrector(name='sza_test', modifiers=tuple())
+        res = comp((self.ds1, self.sza), test_attr='test')
+        np.testing.assert_allclose(res.values, np.array([[22.401667, 22.31777], [22.437503, 22.353533]]))
+
+    def test_basic_lims_provided(self):
+        """Test custom limits when SZA is provided."""
+        from satpy.composites import SunZenithCorrector
+        comp = SunZenithCorrector(name='sza_test', modifiers=tuple(), correction_limit=90)
+        res = comp((self.ds1, self.sza), test_attr='test')
+        np.testing.assert_allclose(res.values, np.array([[66.853262, 68.168939], [66.30742, 67.601493]]))
+
+
 class TestDayNightCompositor(unittest.TestCase):
     """Test DayNightCompositor."""
 
@@ -464,6 +516,7 @@ def suite():
     mysuite.addTests(test_viirs.suite())
     mysuite.addTest(loader.loadTestsFromTestCase(TestCheckArea))
     mysuite.addTest(loader.loadTestsFromTestCase(TestRatioSharpenedCompositors))
+    mysuite.addTest(loader.loadTestsFromTestCase(TestSunZenithCorrector))
     mysuite.addTest(loader.loadTestsFromTestCase(TestDayNightCompositor))
     mysuite.addTest(loader.loadTestsFromTestCase(TestFillingCompositor))
     mysuite.addTest(loader.loadTestsFromTestCase(TestSandwichCompositor))

satpy/utils.py

@@ -31,6 +31,7 @@
 import os
 import re
 
+import numpy as np
 import xarray.ufuncs as xu
 
 try:
@@ -208,56 +209,82 @@ def proj_units_to_meters(proj_str):
 
 
 def _get_sunz_corr_li_and_shibata(cos_zen):
+    return 24.35 / (2. * cos_zen + np.sqrt(498.5225 * cos_zen**2 + 1))
 
-    return 24.35 / (2. * cos_zen +
-                    xu.sqrt(498.5225 * cos_zen**2 + 1))
 
-
-def sunzen_corr_cos(data, cos_zen, limit=88.):
+def sunzen_corr_cos(data, cos_zen, limit=88., max_sza=95.):
     """Perform Sun zenith angle correction.
 
     The correction is based on the provided cosine of the zenith
-    angle (*cos_zen*).  The correction is limited
-    to *limit* degrees (default: 88.0 degrees).  For larger zenith
-    angles, the correction is the same as at the *limit*.  Both *data*
-    and *cos_zen* are given as 2-dimensional Numpy arrays or Numpy
-    MaskedArrays, and they should have equal shapes.
+    angle (``cos_zen``).  The correction is limited
+    to ``limit`` degrees (default: 88.0 degrees).  For larger zenith
+    angles, the correction is the same as at the ``limit`` if ``max_sza``
+    is `None`. The default behavior is to gradually reduce the correction
+    past ``limit`` degrees up to ``max_sza`` where the correction becomes
+    0. Both ``data`` and ``cos_zen`` should be 2D arrays of the same shape.
 
     """
 
     # Convert the zenith angle limit to cosine of zenith angle
-    limit = xu.cos(xu.deg2rad(limit))
+    limit_rad = np.deg2rad(limit)
+    limit_cos = np.cos(limit_rad)
+    max_sza_rad = np.deg2rad(max_sza) if max_sza is not None else max_sza
 
     # Cosine correction
     corr = 1. / cos_zen
-    # Use constant value (the limit) for larger zenith
-    # angles
-    corr = corr.where(cos_zen > limit).fillna(1 / limit)
+    if max_sza is not None:
+        # gradually fall off for larger zenith angle
+        grad_factor = (np.arccos(cos_zen) - limit_rad) / (max_sza_rad - limit_rad)
+        # invert the factor so maximum correction is done at `limit` and falls off later
+        grad_factor = 1. - np.log(grad_factor + 1) / np.log(2)
+        # make sure we don't make anything negative
+        grad_factor = grad_factor.clip(0.)
+    else:
+        # Use constant value (the limit) for larger zenith angles
+        grad_factor = 1.
+    corr = corr.where(cos_zen > limit_cos, grad_factor / limit_cos)
+    # Force "night" pixels to 0 (where SZA is invalid)
+    corr = corr.where(cos_zen.notnull(), 0)
 
     return data * corr
 
 
-def atmospheric_path_length_correction(data, cos_zen, limit=88.):
+def atmospheric_path_length_correction(data, cos_zen, limit=88., max_sza=95.):
     """Perform Sun zenith angle correction.
 
     This function uses the correction method proposed by
     Li and Shibata (2006): https://doi.org/10.1175/JAS3682.1
 
-    The correction is limited to *limit* degrees (default: 88.0 degrees). For
-    larger zenith angles, the correction is the same as at the *limit*. Both
-    *data* and *cos_zen* are given as 2-dimensional Numpy arrays or Numpy
-    MaskedArrays, and they should have equal shapes.
+    The correction is limited to ``limit`` degrees (default: 88.0 degrees). For
+    larger zenith angles, the correction is the same as at the ``limit`` if
+    ``max_sza`` is `None`. The default behavior is to gradually reduce the
+    correction past ``limit`` degrees up to ``max_sza`` where the correction
+    becomes 0. Both ``data`` and ``cos_zen`` should be 2D arrays of the same
+    shape.
 
     """
-
     # Convert the zenith angle limit to cosine of zenith angle
-    limit = xu.cos(xu.radians(limit))
+    limit_rad = np.deg2rad(limit)
+    limit_cos = np.cos(limit_rad)
+    max_sza_rad = np.deg2rad(max_sza) if max_sza is not None else max_sza
 
     # Cosine correction
     corr = _get_sunz_corr_li_and_shibata(cos_zen)
-    # Use constant value (the limit) for larger zenith
-    # angles
-    corr_lim = _get_sunz_corr_li_and_shibata(limit)
-    corr = corr.where(cos_zen > limit).fillna(corr_lim)
+    # Use constant value (the limit) for larger zenith angles
+    corr_lim = _get_sunz_corr_li_and_shibata(limit_cos)
+
+    if max_sza is not None:
+        # gradually fall off for larger zenith angle
+        grad_factor = (np.arccos(cos_zen) - limit_rad) / (max_sza_rad - limit_rad)
+        # invert the factor so maximum correction is done at `limit` and falls off later
+        grad_factor = 1. - np.log(grad_factor + 1) / np.log(2)
+        # make sure we don't make anything negative
+        grad_factor = grad_factor.clip(0.)
+    else:
+        # Use constant value (the limit) for larger zenith angles
+        grad_factor = 1.
+    corr = corr.where(cos_zen > limit_cos, grad_factor * corr_lim)
+    # Force "night" pixels to 0 (where SZA is invalid)
+    corr = corr.where(cos_zen.notnull(), 0)
 
     return data * corr

setup.py

@@ -31,7 +31,7 @@
 
 from setuptools import find_packages, setup
 
-requires = ['numpy >=1.12', 'pillow', 'pyresample >=1.10.3', 'trollsift',
+requires = ['numpy >=1.13', 'pillow', 'pyresample >=1.10.3', 'trollsift',
             'trollimage >=1.5.1', 'pykdtree', 'six', 'pyyaml', 'xarray >=0.10.1',
             'dask[array] >=0.17.1']