From 38bda954fc2a7b80bf394ea6bf71a4616bfb332f Mon Sep 17 00:00:00 2001 From: Johan Strandgren Date: Wed, 23 Aug 2023 14:28:32 +0200 Subject: [PATCH 1/6] Implement non-linearity term for NDVI-weighted hybrid-green correction when converting NDVI to blend fraction. --- satpy/composites/spectral.py | 27 ++++++++++++++++--- satpy/tests/compositor_tests/test_spectral.py | 20 ++++++++++---- 2 files changed, 38 insertions(+), 9 deletions(-) diff --git a/satpy/composites/spectral.py b/satpy/composites/spectral.py index c0ccaff64f..5e6e03c148 100644 --- a/satpy/composites/spectral.py +++ b/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,17 +139,29 @@ 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]) @@ -156,6 +169,12 @@ def __call__(self, projectables, optional_datasets=None, **attrs): ndvi.data = da.where(ndvi > self.ndvi_min, ndvi, self.ndvi_min) ndvi.data = da.where(ndvi < self.ndvi_max, ndvi, self.ndvi_max) + # Apply non-linearity to the ndvi for a non-linear conversion from ndvi to 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) + + # Compute blending fraction from ndvi fraction = (ndvi - self.ndvi_min) / (self.ndvi_max - self.ndvi_min) * (self.limits[1] - self.limits[0]) \ + self.limits[0] self.fractions = (1 - fraction, fraction) diff --git a/satpy/tests/compositor_tests/test_spectral.py b/satpy/tests/compositor_tests/test_spectral.py index 467adf119b..03e51a5043 100644 --- a/satpy/tests/compositor_tests/test_spectral.py +++ b/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 . """Tests for spectral correction compositors.""" -import warnings import dask.array as da import numpy as np @@ -78,6 +77,7 @@ def test_ndvi_hybrid_green(self): 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,12 +86,22 @@ 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) + # Test invalid strength + with pytest.raises(ValueError): + _ = NDVIHybridGreen('ndvi_hybrid_green', strength=0.0, prerequisites=(0.51, 0.65, 0.85), + standard_name='toa_bidirectional_reflectance') + + # Test non-linear strength + 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_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') + 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) From 478a3b5ef6e81bbfb8489cca5df50ef78698670e Mon Sep 17 00:00:00 2001 From: Johan Strandgren Date: Wed, 23 Aug 2023 14:32:47 +0200 Subject: [PATCH 2/6] Modify default strength for FCI green band correction to be in line with EUMETSAT recipe used for first public images. --- satpy/etc/composites/fci.yaml | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/satpy/etc/composites/fci.yaml b/satpy/etc/composites/fci.yaml index 193415656b..30860fbb28 100644 --- a/satpy/etc/composites/fci.yaml +++ b/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] @@ -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 From 6e410f4aa8483b79cdd5900e5e64a93cbcff8882 Mon Sep 17 00:00:00 2001 From: Johan Strandgren Date: Thu, 24 Aug 2023 08:33:12 +0200 Subject: [PATCH 3/6] Add non-linearity term to AHI NDVI green correction. --- satpy/etc/composites/ahi.yaml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/satpy/etc/composites/ahi.yaml b/satpy/etc/composites/ahi.yaml index 9e36bf7d7f..9f771b712c 100644 --- a/satpy/etc/composites/ahi.yaml +++ b/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] From ed1c4a171e703c10dc8d9fa1e4e9260a2bde3abd Mon Sep 17 00:00:00 2001 From: Johan Strandgren Date: Mon, 18 Sep 2023 16:03:28 +0200 Subject: [PATCH 4/6] Add non-linear strength for AMI ndvi_hybrid_green band. --- satpy/etc/composites/ami.yaml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/satpy/etc/composites/ami.yaml b/satpy/etc/composites/ami.yaml index 5a8608d795..b0d943c5ca 100644 --- a/satpy/etc/composites/ami.yaml +++ b/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 From 049a2a3690e07166aa301d68e9d7e9ee569e0d2d Mon Sep 17 00:00:00 2001 From: Johan Strandgren Date: Mon, 18 Sep 2023 16:49:49 +0200 Subject: [PATCH 5/6] Refactor code for applying non-linearity and computing blend fractions. --- satpy/composites/spectral.py | 34 ++++++++++++++++++++++++++++------ 1 file changed, 28 insertions(+), 6 deletions(-) diff --git a/satpy/composites/spectral.py b/satpy/composites/spectral.py index 5e6e03c148..7d05a000d6 100644 --- a/satpy/composites/spectral.py +++ b/satpy/composites/spectral.py @@ -169,17 +169,39 @@ def __call__(self, projectables, optional_datasets=None, **attrs): ndvi.data = da.where(ndvi > self.ndvi_min, ndvi, self.ndvi_min) ndvi.data = da.where(ndvi < self.ndvi_max, ndvi, self.ndvi_max) - # Apply non-linearity to the ndvi for a non-linear conversion from ndvi to 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. + # 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) - # Compute blending fraction from ndvi + 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 . + """ fraction = (ndvi - self.ndvi_min) / (self.ndvi_max - self.ndvi_min) * (self.limits[1] - self.limits[0]) \ + self.limits[0] - self.fractions = (1 - fraction, fraction) - return super().__call__([projectables[0], projectables[2]], **attrs) + return fraction class GreenCorrector(SpectralBlender): From 3767f61a4b2545843232a97c2b7c02dd04448c33 Mon Sep 17 00:00:00 2001 From: Johan Strandgren Date: Mon, 18 Sep 2023 17:03:03 +0200 Subject: [PATCH 6/6] Refactor unit tests for NDVI-weighted hybrid green correction. --- satpy/tests/compositor_tests/test_spectral.py | 46 +++++++++++-------- 1 file changed, 27 insertions(+), 19 deletions(-) diff --git a/satpy/tests/compositor_tests/test_spectral.py b/satpy/tests/compositor_tests/test_spectral.py index 03e51a5043..b6460c911b 100644 --- a/satpy/tests/compositor_tests/test_spectral.py +++ b/satpy/tests/compositor_tests/test_spectral.py @@ -65,8 +65,24 @@ 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), @@ -74,6 +90,8 @@ def test_ndvi_hybrid_green(self): 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') @@ -86,26 +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) - # Test invalid strength - with pytest.raises(ValueError): - _ = NDVIHybridGreen('ndvi_hybrid_green', strength=0.0, prerequisites=(0.51, 0.65, 0.85), - standard_name='toa_bidirectional_reflectance') - - # Test non-linear strength + 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_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) + 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')