Merge pull request #2554 from strandgren/feature_nonlinear_ndvi_green…

…_correction

satpy/composites/spectral.py

@@ -114,8 +114,8 @@ class NDVIHybridGreen(SpectralBlender):
 
     This green band correction follows the same approach as the HybridGreen compositor, but with a dynamic blend
     factor `f` that depends on the pixel-level Normalized Differece Vegetation Index (NDVI). The higher the NDVI, the
-    smaller the contribution from the nir channel will be, following a liner relationship between the two ranges
-    `[ndvi_min, ndvi_max]` and `limits`.
+    smaller the contribution from the nir channel will be, following a liner (default) or non-linear relationship
+    between the two ranges `[ndvi_min, ndvi_max]` and `limits`.
 
     As an example, a new green channel using e.g. FCI data and the NDVIHybridGreen compositor can be defined like::
 
@@ -124,6 +124,7 @@ class NDVIHybridGreen(SpectralBlender):
         ndvi_min: 0.0
         ndvi_max: 1.0
         limits: [0.15, 0.05]
+        strength: 1.0
         prerequisites:
           - name: vis_05
             modifiers: [sunz_corrected, rayleigh_corrected]
@@ -138,30 +139,70 @@ class NDVIHybridGreen(SpectralBlender):
     pixels with NDVI=1.0 will be a weighted average with 5% contribution from the near-infrared
     vis_08 channel and the remaining 95% from the native green vis_05 channel. For other values of
     NDVI a linear interpolation between these values will be performed.
+
+    A strength larger or smaller than 1.0 will introduce a non-linear relationship between the two ranges
+    `[ndvi_min, ndvi_max]` and `limits`. Hence, a higher strength (> 1.0) will result in a slower transition
+    to higher/lower fractions at the NDVI extremes. Similarly, a lower strength (< 1.0) will result in a
+    faster transition to higher/lower fractions at the NDVI extremes.
     """
 
-    def __init__(self, *args, ndvi_min=0.0, ndvi_max=1.0, limits=(0.15, 0.05), **kwargs):
-        """Initialize class and set the NDVI limits and the corresponding blending fraction limits."""
+    def __init__(self, *args, ndvi_min=0.0, ndvi_max=1.0, limits=(0.15, 0.05), strength=1.0, **kwargs):
+        """Initialize class and set the NDVI limits, blending fraction limits and strength."""
+        if strength <= 0.0:
+            raise ValueError(f"Expected stength greater than 0.0, got {strength}.")
+
         self.ndvi_min = ndvi_min
         self.ndvi_max = ndvi_max
         self.limits = limits
+        self.strength = strength
         super().__init__(*args, **kwargs)
 
     def __call__(self, projectables, optional_datasets=None, **attrs):
         """Construct the hybrid green channel weighted by NDVI."""
+        LOG.info(f"Applying NDVI-weighted hybrid-green correction with limits [{self.limits[0]}, "
+                 f"{self.limits[1]}] and strength {self.strength}.")
+
         ndvi_input = self.match_data_arrays([projectables[1], projectables[2]])
 
         ndvi = (ndvi_input[1] - ndvi_input[0]) / (ndvi_input[1] + ndvi_input[0])
 
         ndvi.data = da.where(ndvi > self.ndvi_min, ndvi, self.ndvi_min)
         ndvi.data = da.where(ndvi < self.ndvi_max, ndvi, self.ndvi_max)
 
-        fraction = (ndvi - self.ndvi_min) / (self.ndvi_max - self.ndvi_min) * (self.limits[1] - self.limits[0]) \
-            + self.limits[0]
+        # Introduce non-linearity to ndvi for non-linear scaling to NIR blend fraction
+        if self.strength != 1.0:  # self._apply_strength() has no effect if strength = 1.0 -> no non-linear behaviour
+            ndvi = self._apply_strength(ndvi)
+
+        # Compute pixel-level NIR blend fractions from ndvi
+        fraction = self._compute_blend_fraction(ndvi)
+
+        # Prepare input as required by parent class (SpectralBlender)
         self.fractions = (1 - fraction, fraction)
 
         return super().__call__([projectables[0], projectables[2]], **attrs)
 
+    def _apply_strength(self, ndvi):
+        """Introduce non-linearity by applying strength factor.
+
+        The method introduces non-linearity to the ndvi for a non-linear scaling from ndvi to blend fraction in
+        `_compute_blend_fraction`. 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 ndvi.
+        """
+        ndvi = ndvi ** self.strength / (ndvi ** self.strength + (1 - ndvi) ** self.strength)
+
+        return ndvi
+
+    def _compute_blend_fraction(self, ndvi):
+        """Compute pixel-level fraction of NIR signal to blend with native green signal.
+
+        This method linearly scales the input ndvi values to pixel-level blend fractions within the range
+        `[limits[0], limits[1]]` following this implementation <https://stats.stackexchange.com/a/281164>.
+        """
+        fraction = (ndvi - self.ndvi_min) / (self.ndvi_max - self.ndvi_min) * (self.limits[1] - self.limits[0]) \
+            + self.limits[0]
+
+        return fraction
+
 
 class GreenCorrector(SpectralBlender):
     """Previous class used to blend channels for green band corrections.

satpy/etc/composites/ahi.yaml

@@ -105,6 +105,8 @@ composites:
 
   ndvi_hybrid_green:
     compositor: !!python/name:satpy.composites.spectral.NDVIHybridGreen
+    limits: [0.15, 0.05]
+    strength: 3.0
     prerequisites:
     - name: B02
       modifiers: [sunz_corrected, rayleigh_corrected, sunz_reduced]

satpy/etc/composites/ami.yaml

@@ -69,6 +69,7 @@ composites:
       currently implemented are experimental and may change in future versions of Satpy.
     compositor: !!python/name:satpy.composites.spectral.NDVIHybridGreen
     limits: [0.15, 0.05]
+    strength: 3.0
     prerequisites:
       - name: VI005
         modifiers: [sunz_corrected, rayleigh_corrected]
@@ -83,6 +84,7 @@ composites:
       Alternative to ndvi_hybrid_green, but without solar zenith or rayleigh correction.
     compositor: !!python/name:satpy.composites.spectral.NDVIHybridGreen
     limits: [0.15, 0.05]
+    strength: 3.0
     prerequisites:
       - name: VI005
       - name: VI006

satpy/etc/composites/fci.yaml

@@ -7,10 +7,11 @@ composites:
       The FCI green band at 0.51 µm deliberately misses the chlorophyll band, such that
       the signal comes from aerosols and ash rather than vegetation. An effect
       is that vegetation in a true colour RGB looks rather brown than green and barren rather red. Mixing in
-      some part of the NIR 0.8 channel reduced this effect.  Note that the fractions
+      some part of the NIR 0.8 channel reduced this effect.  Note that the fractions and non-linear strength
       currently implemented are experimental and may change in future versions of Satpy.
     compositor: !!python/name:satpy.composites.spectral.NDVIHybridGreen
     limits: [0.15, 0.05]
+    strength: 3.0
     prerequisites:
       - name: vis_05
         modifiers: [sunz_corrected, rayleigh_corrected, sunz_reduced]
@@ -25,6 +26,7 @@ composites:
       Alternative to ndvi_hybrid_green, but without solar zenith or rayleigh correction.
     compositor: !!python/name:satpy.composites.spectral.NDVIHybridGreen
     limits: [0.15, 0.05]
+    strength: 3.0
     prerequisites:
       - name: vis_05
       - name: vis_06

satpy/tests/compositor_tests/test_spectral.py

@@ -14,7 +14,6 @@
 # You should have received a copy of the GNU General Public License along with
 # satpy.  If not, see <http://www.gnu.org/licenses/>.
 """Tests for spectral correction compositors."""
-import warnings
 
 import dask.array as da
 import numpy as np
@@ -66,18 +65,37 @@ def test_hybrid_green(self):
         data = res.compute()
         np.testing.assert_allclose(data, 0.23)
 
-    def test_ndvi_hybrid_green(self):
-        """Test NDVI-scaled hybrid green correction of 'green' band."""
+    def test_green_corrector(self):
+        """Test the deprecated class for green corrections."""
+        comp = GreenCorrector('blended_channel', fractions=(0.85, 0.15), prerequisites=(0.51, 0.85),
+                              standard_name='toa_bidirectional_reflectance')
+        res = comp((self.c01, self.c03))
+        assert isinstance(res, xr.DataArray)
+        assert isinstance(res.data, da.Array)
+        assert res.attrs['name'] == 'blended_channel'
+        assert res.attrs['standard_name'] == 'toa_bidirectional_reflectance'
+        data = res.compute()
+        np.testing.assert_allclose(data, 0.23)
+
+
+class TestNdviHybridGreenCompositor:
+    """Test NDVI-weighted hybrid green correction of green band."""
+
+    def setup_method(self):
+        """Initialize channels."""
         self.c01 = xr.DataArray(da.from_array([[0.25, 0.30], [0.20, 0.30]], chunks=25),
                                 dims=('y', 'x'), attrs={'name': 'C02'})
         self.c02 = xr.DataArray(da.from_array([[0.25, 0.30], [0.25, 0.35]], chunks=25),
                                 dims=('y', 'x'), attrs={'name': 'C03'})
         self.c03 = xr.DataArray(da.from_array([[0.35, 0.35], [0.28, 0.65]], chunks=25),
                                 dims=('y', 'x'), attrs={'name': 'C04'})
 
+    def test_ndvi_hybrid_green(self):
+        """Test General functionality with linear scaling from ndvi to blend fraction."""
         comp = NDVIHybridGreen('ndvi_hybrid_green', limits=(0.15, 0.05), prerequisites=(0.51, 0.65, 0.85),
                                standard_name='toa_bidirectional_reflectance')
 
+        # Test General functionality with linear strength (=1.0)
         res = comp((self.c01, self.c02, self.c03))
         assert isinstance(res, xr.DataArray)
         assert isinstance(res.data, da.Array)
@@ -86,16 +104,16 @@ def test_ndvi_hybrid_green(self):
         data = res.values
         np.testing.assert_array_almost_equal(data, np.array([[0.2633, 0.3071], [0.2115, 0.3420]]), decimal=4)
 
-    def test_green_corrector(self):
-        """Test the deprecated class for green corrections."""
-        with warnings.catch_warnings():
-            warnings.filterwarnings("ignore", category=UserWarning, message=r'.*deprecated.*')
-            comp = GreenCorrector('blended_channel', fractions=(0.85, 0.15), prerequisites=(0.51, 0.85),
-                                  standard_name='toa_bidirectional_reflectance')
-        res = comp((self.c01, self.c03))
-        assert isinstance(res, xr.DataArray)
-        assert isinstance(res.data, da.Array)
-        assert res.attrs['name'] == 'blended_channel'
-        assert res.attrs['standard_name'] == 'toa_bidirectional_reflectance'
-        data = res.compute()
-        np.testing.assert_allclose(data, 0.23)
+    def test_nonliniear_scaling(self):
+        """Test non-linear scaling using `strength` term."""
+        comp = NDVIHybridGreen('ndvi_hybrid_green', limits=(0.15, 0.05), strength=2.0, prerequisites=(0.51, 0.65, 0.85),
+                               standard_name='toa_bidirectional_reflectance')
+
+        res = comp((self.c01, self.c02, self.c03))
+        np.testing.assert_array_almost_equal(res.values, np.array([[0.2646, 0.3075], [0.2120, 0.3471]]), decimal=4)
+
+    def test_invalid_strength(self):
+        """Test using invalid `strength` term for non-linear scaling."""
+        with pytest.raises(ValueError):
+            _ = NDVIHybridGreen('ndvi_hybrid_green', strength=0.0, prerequisites=(0.51, 0.65, 0.85),
+                                standard_name='toa_bidirectional_reflectance')