Merge pull request #2556 from strandgren/feature_sunz_reduction

satpy/etc/composites/ahi.yaml

@@ -107,11 +107,11 @@ composites:
     compositor: !!python/name:satpy.composites.spectral.NDVIHybridGreen
     prerequisites:
     - name: B02
-      modifiers: [sunz_corrected, rayleigh_corrected]
+      modifiers: [sunz_corrected, rayleigh_corrected, sunz_reduced]
     - name: B03
-      modifiers: [sunz_corrected, rayleigh_corrected]
+      modifiers: [sunz_corrected, rayleigh_corrected, sunz_reduced]
     - name: B04
-      modifiers: [sunz_corrected]
+      modifiers: [sunz_corrected, sunz_reduced]
     standard_name: toa_bidirectional_reflectance
 
   airmass:
@@ -275,10 +275,10 @@ composites:
     compositor: !!python/name:satpy.composites.SelfSharpenedRGB
     prerequisites:
     - name: B03
-      modifiers: [sunz_corrected, rayleigh_corrected]
+      modifiers: [sunz_corrected, rayleigh_corrected, sunz_reduced]
     - name: ndvi_hybrid_green
     - name: B01
-      modifiers: [sunz_corrected, rayleigh_corrected]
+      modifiers: [sunz_corrected, rayleigh_corrected, sunz_reduced]
     high_resolution_band: red
     standard_name: true_color
 

satpy/etc/composites/fci.yaml

@@ -13,11 +13,11 @@ composites:
     limits: [0.15, 0.05]
     prerequisites:
       - name: vis_05
-        modifiers: [sunz_corrected, rayleigh_corrected]
+        modifiers: [sunz_corrected, rayleigh_corrected, sunz_reduced]
       - name: vis_06
-        modifiers: [sunz_corrected, rayleigh_corrected]
+        modifiers: [sunz_corrected, rayleigh_corrected, sunz_reduced]
       - name: vis_08
-        modifiers: [sunz_corrected ]
+        modifiers: [sunz_corrected, sunz_reduced ]
     standard_name: toa_bidirectional_reflectance
 
   ndvi_hybrid_green_raw:
@@ -48,10 +48,10 @@ composites:
       of the ndvi_hybrid_green composites for details.
     prerequisites:
       - name: vis_06
-        modifiers: [sunz_corrected, rayleigh_corrected]
+        modifiers: [sunz_corrected, rayleigh_corrected, sunz_reduced]
       - name: ndvi_hybrid_green
       - name: vis_04
-        modifiers: [sunz_corrected, rayleigh_corrected]
+        modifiers: [sunz_corrected, rayleigh_corrected, sunz_reduced]
     standard_name: true_color
 
   true_color_raw_with_corrected_green:

satpy/etc/composites/visir.yaml

@@ -16,6 +16,11 @@ modifiers:
     optional_prerequisites:
       - solar_zenith_angle
 
+  sunz_reduced:
+    modifier: !!python/name:satpy.modifiers.SunZenithReducer
+    optional_prerequisites:
+    - solar_zenith_angle
+
   co2_corrected:
     modifier: !!python/name:satpy.modifiers.CO2Corrector
     prerequisites:

satpy/modifiers/__init__.py

@@ -25,5 +25,6 @@
 from .atmosphere import PSPRayleighReflectance  # noqa: F401
 from .geometry import EffectiveSolarPathLengthCorrector  # noqa: F401
 from .geometry import SunZenithCorrector  # noqa: F401
+from .geometry import SunZenithReducer  # noqa: F401
 from .spectral import NIREmissivePartFromReflectance  # noqa: F401
 from .spectral import NIRReflectance  # noqa: F401

satpy/modifiers/angles.py

@@ -235,13 +235,15 @@ def cache_to_zarr_if(
     out old entries. It is up to the user to manage the size of the cache.
 
     """
+
     def _decorator(func: Callable) -> Callable:
         zarr_cacher = ZarrCacheHelper(func,
                                       cache_config_key,
                                       uncacheable_arg_types,
                                       sanitize_args_func)
         wrapper = update_wrapper(zarr_cacher, func)
         return wrapper
+
     return _decorator
 
 
@@ -547,3 +549,44 @@ def _sunzen_corr_cos_ndarray(data: np.ndarray,
     # Force "night" pixels to 0 (where SZA is invalid)
     corr[np.isnan(cos_zen)] = 0
     return data * corr
+
+
+def sunzen_reduction(data: da.Array,
+                     sunz: da.Array,
+                     limit: float = 55.,
+                     max_sza: float = 90.,
+                     strength: float = 1.5) -> da.Array:
+    """Reduced strength of signal at high sun zenith angles."""
+    return da.map_blocks(_sunzen_reduction_ndarray, data, sunz, limit, max_sza, strength,
+                         meta=np.array((), dtype=data.dtype), chunks=data.chunks)
+
+
+def _sunzen_reduction_ndarray(data: np.ndarray,
+                              sunz: np.ndarray,
+                              limit: float,
+                              max_sza: float,
+                              strength: float) -> np.ndarray:
+    # compute reduction factor (0.0 - 1.0) between limit and maz_sza
+    reduction_factor = (sunz - limit) / (max_sza - limit)
+    reduction_factor = reduction_factor.clip(0., 1.)
+
+    # invert the reduction factor such that minimum reduction is done at `limit` and gradually increases towards max_sza
+    with np.errstate(invalid='ignore'):  # we expect space pixels to be invalid
+        reduction_factor = 1. - np.log(reduction_factor + 1) / np.log(2)
+
+    # apply non-linearity to the reduction factor for a non-linear reduction of the signal. This can be used for a
+    # slower or faster transision to higher/lower fractions at the ndvi extremes. If strength equals 1.0, this
+    # operation has no effect on the reduction_factor.
+    reduction_factor = reduction_factor ** strength / (
+                reduction_factor ** strength + (1 - reduction_factor) ** strength)
+
+    # compute final correction term, with no reduction for angles < limit
+    corr = np.where(sunz < limit, 1.0, reduction_factor)
+
+    # force "night" pixels to 0 (where SZA is invalid)
+    corr[np.isnan(sunz)] = 0
+
+    # reduce data signal with correction term
+    res = data * corr
+
+    return res

satpy/modifiers/geometry.py

@@ -24,7 +24,7 @@
 import numpy as np
 
 from satpy.modifiers import ModifierBase
-from satpy.modifiers.angles import sunzen_corr_cos
+from satpy.modifiers.angles import sunzen_corr_cos, sunzen_reduction
 from satpy.utils import atmospheric_path_length_correction
 
 logger = logging.getLogger(__name__)
@@ -159,3 +159,46 @@ def __init__(self, correction_limit=88., **kwargs):
     def _apply_correction(self, proj, coszen):
         logger.debug("Apply the effective solar atmospheric path length correction method by Li and Shibata")
         return atmospheric_path_length_correction(proj, coszen, limit=self.correction_limit, max_sza=self.max_sza)
+
+
+class SunZenithReducer(SunZenithCorrectorBase):
+    """Reduce signal strength at large sun zenith angles.
+
+    Within a given sunz interval [correction_limit, max_sza] the strength of the signal is reduced following the
+    formula:
+
+      res = signal * reduction_factor
+
+    where reduction_factor is a pixel-level value ranging from 0 to 1 within the sunz interval.
+
+    The `strength` parameter can be used for a non-linear reduction within the sunz interval. A strength larger
+    than 1.0 will decelerate the signal reduction towards the sunz interval extremes, whereas a strength
+    smaller than 1.0 will accelerate the signal reduction towards the sunz interval extremes.
+
+    """
+
+    def __init__(self, correction_limit=55., max_sza=90, strength=1.5, **kwargs):
+        """Collect custom configuration values.
+
+        Args:
+            correction_limit (float): Solar zenith angle in degrees where to start the signal reduction. Default 60.
+            max_sza (float): Maximum solar zenith angle in degrees where to apply the signal reduction. Beyond
+                             this solar zenith angle the signal will become zero. Default 90.
+            strength (float): The strength of the non-linear signal reduction. Default 1.5
+
+        """
+        self.correction_limit = correction_limit
+        self.strength = strength
+        super(SunZenithReducer, self).__init__(max_sza=max_sza, **kwargs)
+        if self.max_sza is None:
+            raise ValueError("`max_sza` must be defined when using the SunZenithReducer.")
+
+    def _apply_correction(self, proj, coszen):
+        logger.debug("Apply sun-zenith signal reduction")
+        res = proj.copy()
+        sunz = np.rad2deg(np.arccos(coszen.data))
+        res.data = sunzen_reduction(proj.data, sunz,
+                                    limit=self.correction_limit,
+                                    max_sza=self.max_sza,
+                                    strength=self.strength)
+        return res

satpy/tests/test_modifiers.py

@@ -162,6 +162,38 @@ def test_imcompatible_areas(self, sunz_ds2, sunz_sza):
             comp((sunz_ds2, sunz_sza), test_attr='test')
 
 
+class TestSunZenithReducer:
+    """Test case for the sun zenith reducer."""
+
+    @classmethod
+    def setup_class(cls):
+        """Initialze SunZenithReducer classes that shall be tested."""
+        from satpy.modifiers.geometry import SunZenithReducer
+        cls.default = SunZenithReducer(name='sza_reduction_test_default', modifiers=tuple())
+        cls.custom = SunZenithReducer(name='sza_reduction_test_custom', modifiers=tuple(),
+                                      correction_limit=70, max_sza=95, strength=3.0)
+
+    def test_default_settings(self, sunz_ds1, sunz_sza):
+        """Test default settings with sza data available."""
+        res = self.default((sunz_ds1, sunz_sza), test_attr='test')
+        np.testing.assert_allclose(res.values,
+                                   np.array([[0.00242814, 0.00235669], [0.00245885, 0.00238707]]),
+                                   rtol=1e-5)
+
+    def test_custom_settings(self, sunz_ds1, sunz_sza):
+        """Test custom settings with sza data available."""
+        res = self.custom((sunz_ds1, sunz_sza), test_attr='test')
+        np.testing.assert_allclose(res.values,
+                                   np.array([[0.01041319, 0.01030033], [0.01046164, 0.01034834]]),
+                                   rtol=1e-5)
+
+    def test_invalid_max_sza(self, sunz_ds1, sunz_sza):
+        """Test invalid max_sza with sza data available."""
+        from satpy.modifiers.geometry import SunZenithReducer
+        with pytest.raises(ValueError):
+            SunZenithReducer(name='sza_reduction_test_invalid', modifiers=tuple(), max_sza=None)
+
+
 class TestNIRReflectance(unittest.TestCase):
     """Test NIR reflectance compositor."""