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

Jump to bottom

Add NDVI-scaled hybrid green correction #2280

Merged
merged 19 commits into from
Dec 15, 2022
Merged

Add NDVI-scaled hybrid green correction #2280

Show file tree
Hide file tree
Changes from 15 commits
Commits
Show all changes
19 commits
Select commit Hold shift + click to select a range
848d818
Add NDVI-scaled hybrid green correction.
strandgren Nov 16, 2022
5903d93
Add unit test for NDVIHybridGreen compositor.
strandgren Nov 16, 2022
134dcdc
Refactor code used for green band corrections.
strandgren Nov 17, 2022
9c159d8
Use NDVI weighted hybrid green method for green band correction of FC…
strandgren Nov 17, 2022
0393031
Add missing input data for NDVI hybrid green correction.
strandgren Nov 17, 2022
a73c607
Remove Rayleigh correction for 0.8 micron channel and fix minor issue…
strandgren Nov 17, 2022
17ac671
Replace unit-test with pytest.
strandgren Nov 18, 2022
3d6f2e8
Ensure backwards compatibility.
strandgren Nov 18, 2022
4aeb090
Revert unintended change to Rayleigh correction.
strandgren Nov 18, 2022
caed5fb
Re-introduce satpy/composites/ahi.py in order to ensure compatibility…
strandgren Nov 18, 2022
c23a15e
Fix intendation of composite comments.
strandgren Nov 18, 2022
7b282bb
Fix documentation build errors.
strandgren Nov 21, 2022
34c7731
Fix typos.
strandgren Nov 23, 2022
2f4b419
Harmonize use of white space in lists.
strandgren Nov 23, 2022
0f76c27
Harmonize use of white space in lists.
strandgren Nov 23, 2022
643e252
Fix shebang
strandgren Nov 28, 2022
bf345dd
Clarify docstrings.
strandgren Nov 28, 2022
917250a
Merge branch 'feature-ndvi_corrected_hybrid_green' of https://github.…
strandgren Nov 28, 2022
befb402
Fix typo.
strandgren Dec 6, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
163 changes: 136 additions & 27 deletions satpy/composites/spectral.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,55 +16,164 @@
"""Composite classes for spectral adjustments."""

import logging
import warnings

import dask.array as da

from satpy.composites import GenericCompositor
from satpy.dataset import combine_metadata

LOG = logging.getLogger(__name__)


class GreenCorrector(GenericCompositor):
class SpectralBlender(GenericCompositor):
"""Construct new channel by blending contributions from a set of channels.

This class can be used to compute weighted average of different channels.
Primarily it's used to correct the green band of AHI and FCI in order to
allow for proper true color imagery.

Below is an example used to generate a corrected green channel for AHI using a weighted average from
three channels, with 63% contribution from the native green channel (B02), 29% from the red channel (B03)
and 8% from the near-infrared channel (B04)::

corrected_green:
compositor: !!python/name:satpy.composites.spectral.SpectralBlender
fractions: [0.63, 0.29, 0.08]
prerequisites:
- name: B02
modifiers: [sunz_corrected, rayleigh_corrected]
- name: B03
modifiers: [sunz_corrected, rayleigh_corrected]
- name: B04
modifiers: [sunz_corrected, rayleigh_corrected]
standard_name: toa_bidirectional_reflectance

Other examples can be found in the``ahi.yaml`` composite file in the satpy distribution.
"""

def __init__(self, *args, fractions=(), **kwargs):
"""Set default keyword argument values."""
self.fractions = fractions
super().__init__(*args, **kwargs)

def __call__(self, projectables, optional_datasets=None, **attrs):
"""Blend channels in projectables using the weights in self.fractions."""
if len(self.fractions) != len(projectables):
raise ValueError("fractions and projectables must have the same length.")

projectables = self.match_data_arrays(projectables)
new_channel = sum(fraction * value for fraction, value in zip(self.fractions, projectables))
new_channel.attrs = combine_metadata(*projectables)
return super().__call__((new_channel,), **attrs)


class HybridGreen(SpectralBlender):
"""Corrector of the FCI or AHI green band.

The green band in FCI and AHI deliberately misses the chlorophyll peak
in order to focus on aerosol and ash rather than on vegetation. This
affects true colour RGBs, because vegetation looks brown rather than green.
To make vegetation look greener again, this corrector allows
to simulate the green band as a fraction of two or more other channels.
The green band in FCI and AHI (and other bands centered at 0.51 microns) deliberately
misses the chlorophyll spectral reflectance local maximum at 0.55 microns
in order to focus on aerosol and ash rather than on vegetation. This
affects true colour RGBs, because vegetation looks brown rather than green
and barren surface types typically gets a reddish hue.

To be used, the composite takes two or more input channels and a parameter
``fractions`` that should be a list of floats with the same length as the
number of channels.
To correct for this the hybrid green approach proposed by Miller et al. (2016, :doi:`10.1175/BAMS-D-15-00154.2`)
is used. The basic idea is to include some contribution also from the 0.86 micron
channel, which is known for its sensitivity to vegetation. The formula used for this is::

For example, to simulate an FCI corrected green composite, one could use
a combination of 93% from the green band (vis_05) and 7% from the
near-infrared 0.8 µm band (vis_08)::
hybrid_green = (1 - F) * R(0.51) + F * R(0.86)

corrected_green:
compositor: !!python/name:satpy.composites.ahi.GreenCorrector
fractions: [0.93, 0.07]
where F is a constant value. The validation against VIIRS presented in the paper suggests
F = 0.07 as an optimal value to reconstruct a green band optimized for true color imagery
from AHI. This is the default value currently used for the hybrid green correction in Satpy.

For example, to construct the hybrid green for AHI, one should use
a combination of 93% from the green band (B02) and 7% from the
near-infrared 0.85 µm band (B04)::
strandgren marked this conversation as resolved.
Show resolved Hide resolved

hybrid_green:
compositor: !!python/name:satpy.composites.spectral.HybridGreen
fraction: 0.07
prerequisites:
- name: vis_05
- name: B02
modifiers: [sunz_corrected, rayleigh_corrected]
- name: vis_08
- name: B04
modifiers: [sunz_corrected, rayleigh_corrected]
standard_name: toa_bidirectional_reflectance

Other examples can be found in the ``fci.yaml`` and ``ahi.yaml`` composite
files in the satpy distribution.
"""

def __init__(self, *args, fractions=(0.85, 0.15), **kwargs):
def __init__(self, *args, fraction=0.15, **kwargs):
"""Set default keyword argument values."""
# XXX: Should this be 0.93 and 0.07
self.fractions = fractions
super(GreenCorrector, self).__init__(*args, **kwargs)
fractions = (1 - fraction, fraction)
super().__init__(fractions=fractions, *args, **kwargs)


class NDVIHybridGreen(SpectralBlender):
"""Construct a NDVI-weighted hybrid green channel.

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`.

A new green channel using e.g. FCI data and the NDVIHybridGreen compositor can be defined like::

ndvi_hybrid_green:
compositor: !!python/name:satpy.composites.spectral.NDVIHybridGreen
limits: [0.15, 0.05]
prerequisites:
- name: vis_05
modifiers: [sunz_corrected, rayleigh_corrected]
- name: vis_06
modifiers: [sunz_corrected, rayleigh_corrected]
strandgren marked this conversation as resolved.
Show resolved Hide resolved
- name: vis_08
modifiers: [sunz_corrected ]
standard_name: toa_bidirectional_reflectance

In this example, pixels with NDVI=0.0 (default `ndvi_min`) will be a weighted average with 85% contribution from the
Copy link
Member

Choose a reason for hiding this comment

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

Is this correct? I thought that for very green objects we wanted more of the 0.8 channel, not less - as we miss the chlorophyll peak in the FCI/AHI 0.5 micron band we need to boost it.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Well this is the interesting bit with the green band correction. The location of the AHI/FCI green bands has two effects:

  1. Vegetation appears brownish.
  2. Barren appears reddish.

The green color of vegetation can be boosted with a higher contribution from the NIR channel. Similarly, the reddishness of barren is reduced with an increasing contribution from the NIR channel. However, more contribution from NIR is required to get the proper color of barren, compared to vegetation. Hence, I came up with the solution of increasing the contribution from NIR with decreasing NDVI.

Below are a few examples using the default hybrid green method, with decreasing contribution from the NIR channel. As you can see the color of Australia starts to get too red when the color of vegetation start to look realistic. I also included the result of the NDVIHybridGreen (with just qualitative tuning) and VIIRS RGB from NASA worldview as reference.

15% NIR 10% NIR 7% NIR 4% NIR 0% NIR With this PR VIIRS from NASA Worldview
AHI_true_color_hybrid_green_85_15 AHI_true_color_hybrid_green_90_10 AHI_true_color_hybrid_green_93_7 AHI_true_color_hybrid_green_96_4 AHI_true_color_hybrid_green_100_0 AHI_true_color_ndvi_green_15_5 image

Copy link
Member

Choose a reason for hiding this comment

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

Hence, I came up with the solution of increasing the contribution from NIR with decreasing NDVI.

Your solution, does this mean that this is just for FCI so is "new" or is this being applied to AHI by some other group and documented in some paper?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

This is only internal developments in preparation for FCI. We are testing with AHI since we have the same scenario with sub-optimal position of the green band there. Hopefully we can publish it at some point though.

native green vis_05 channel and 15% from the near-infrared vis_08 channel, whereas pixels with an NDVI=1.0 (default
`ndvi_max`) will be a weighted average with 95% contribution from the native green vis_05 channel and 5% from the
near-infrared vis_08 channel. For other values of NDVI (within this range) a linear interpolation will be performed.
"""

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."""
self.ndvi_min = ndvi_min
self.ndvi_max = ndvi_max
self.limits = limits
super().__init__(*args, **kwargs)

def __call__(self, projectables, optional_datasets=None, **attrs):
"""Boost vegetation effect thanks to NIR (0.8µm) band."""
LOG.info('Boosting vegetation on green band')
"""Construct the hybrid green channel weighted by NDVI."""
ndvi_input = self.match_data_arrays([projectables[1], projectables[2]])

projectables = self.match_data_arrays(projectables)
new_green = sum(fraction * value for fraction, value in zip(self.fractions, projectables))
new_green.attrs = combine_metadata(*projectables)
return super(GreenCorrector, self).__call__((new_green,), **attrs)
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]
self.fractions = (1 - fraction, fraction)

return super().__call__([projectables[0], projectables[2]], **attrs)


class GreenCorrector(SpectralBlender):
"""Previous class used to blend channels for green band corrections.

This method has been refactored to make it more generic. The replacement class is 'SpectralBlender' which computes
a weighted average based on N number of channels and N number of corresponding weights/fractions. A new class
called 'HybridGreen' has been created, which performs a correction of green bands centered at 0.51 microns
following Miller et al. (2016, :doi:`10.1175/BAMS-D-15-00154.2`) in order to improve true color imagery.
"""

def __init__(self, *args, fractions=(0.85, 0.15), **kwargs):
"""Set default keyword argument values."""
warnings.warn(
"'GreenCorrector' is deprecated, use 'SpectralBlender' instead, or 'HybridGreen' for hybrid green"
" correction following Miller et al. (2016).", RuntimeWarning)
super().__init__(fractions=fractions, *args, **kwargs)
77 changes: 70 additions & 7 deletions satpy/etc/composites/ahi.yaml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,9 +16,50 @@ modifiers:

composites:
green:
# Deprecated green composite using deprecated class 'GreenCorrector'. Use 'hybrid_green' instead.
compositor: !!python/name:satpy.composites.spectral.GreenCorrector
# FUTURE: Set a wavelength...see what happens. Dependency finding
# probably wouldn't work.
prerequisites:
# should we be using the most corrected or least corrected inputs?
# what happens if something requests more modifiers on top of this?
mraspaud marked this conversation as resolved.
Show resolved Hide resolved
- wavelength: 0.51
modifiers: [sunz_corrected, rayleigh_corrected]
- wavelength: 0.85
modifiers: sunz_corrected]
Copy link
Member

Choose a reason for hiding this comment

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

Something is not being tested. This syntax is not correct.

Copy link
Collaborator Author

@strandgren strandgren Dec 6, 2022 •
edited

Choose a reason for hiding this comment

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

Indeed this syntax leads to an error, with sunz_corrected] getting split into single characters:

[DEBUG: 2022-12-06 16:01:05 : satpy.node] Skipping optional DataQuery(name='solar_zenith_angle', calibration='reflectance'): Unknown dataset DataQuery(name='solar_zenith_angle', calibration='reflectance')
Traceback (most recent call last):
  File "/tcenas/home/strandgren/git/ext/py/satpy/latest/satpy/satpy/scene.py", line 1362, in _update_dependency_tree
    self._dependency_tree.populate_with_keys(needed_datasets, query)
  File "/tcenas/home/strandgren/git/ext/py/satpy/latest/satpy/satpy/dependency_tree.py", line 262, in populate_with_keys
    raise MissingDependencies(unknown_datasets, "Unknown datasets:")
satpy.node.MissingDependencies: Unknown datasets: {DataQuery(name='green'): {DataQuery(wavelength=0.85, modifiers=('s', 'u', 'n', 'z', '_', 'c', 'o', 'r', 'r', 'e', 'c', 't', 'e', 'd', ']')): {DataQuery(wavelength=0.85, modifiers=('s',))}}}

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/tcenas/home/strandgren/git/ext/py/satpy/latest/satpy/satpy/scene.py", line 1364, in _update_dependency_tree
    raise KeyError(str(err))
KeyError: "Unknown datasets: {DataQuery(name='green'): {DataQuery(wavelength=0.85, modifiers=('s', 'u', 'n', 'z', '_', 'c', 'o', 'r', 'r', 'e', 'c', 't', 'e', 'd', ']')): {DataQuery(wavelength=0.85, modifiers=('s',))}}}"

I have fixed this now, but I'm not sure what kind of (unit) test that would detect this kind of syntax? Did I miss something?

And regarding Polar2Grid you can wait with merging this if you want, no hurry from my side. Although, I have been very careful in making sure that it's fully backwards compatible, but I cannot not guarantee that I have not missed anything. For example the ratios of the green correction are not changing with this PR, the AHI true color images should look like before if the same composite names are used.

standard_name: toa_bidirectional_reflectance

green_true_color_reproduction:
# Deprecated green composite using deprecated class 'GreenCorrector'. Use 'reproduced_green' instead.
# JMA True Color Reproduction green band
# http://www.jma.go.jp/jma/jma-eng/satellite/introduction/TCR.html
compositor: !!python/name:satpy.composites.spectral.GreenCorrector
fractions: [0.6321, 0.2928, 0.0751]
prerequisites:
- name: B02
modifiers: [sunz_corrected, rayleigh_corrected]
- name: B03
modifiers: [sunz_corrected, rayleigh_corrected]
- name: B04
modifiers: [sunz_corrected]
standard_name: none

green_nocorr:
# Deprecated green composite using deprecated class 'GreenCorrector'. Use 'hybrid_green_nocorr' instead.
compositor: !!python/name:satpy.composites.spectral.GreenCorrector
# FUTURE: Set a wavelength...see what happens. Dependency finding
# probably wouldn't work.
prerequisites:
# should we be using the most corrected or least corrected inputs?
# what happens if something requests more modifiers on top of this?
- wavelength: 0.51
- wavelength: 0.85
standard_name: toa_reflectance

hybrid_green:
compositor: !!python/name:satpy.composites.spectral.HybridGreen
# FUTURE: Set a wavelength...see what happens. Dependency finding
# probably wouldn't work.
prerequisites:
# should we be using the most corrected or least corrected inputs?
# what happens if something requests more modifiers on top of this?
Expand All @@ -28,10 +69,10 @@ composites:
modifiers: [sunz_corrected]
standard_name: toa_bidirectional_reflectance

green_true_color_reproduction:
reproduced_green:
# JMA True Color Reproduction green band
# http://www.jma.go.jp/jma/jma-eng/satellite/introduction/TCR.html
compositor: !!python/name:satpy.composites.spectral.GreenCorrector
compositor: !!python/name:satpy.composites.spectral.SpectralBlender
fractions: [0.6321, 0.2928, 0.0751]
prerequisites:
- name: B02
Expand All @@ -42,8 +83,8 @@ composites:
modifiers: [sunz_corrected]
standard_name: none

green_nocorr:
compositor: !!python/name:satpy.composites.spectral.GreenCorrector
hybrid_green_nocorr:
compositor: !!python/name:satpy.composites.spectral.HybridGreen
# FUTURE: Set a wavelength...see what happens. Dependency finding
# probably wouldn't work.
prerequisites:
Expand All @@ -53,6 +94,17 @@ composites:
- wavelength: 0.85
standard_name: toa_reflectance

ndvi_hybrid_green:
compositor: !!python/name:satpy.composites.spectral.NDVIHybridGreen
prerequisites:
- name: B02
modifiers: [sunz_corrected, rayleigh_corrected]
- name: B03
modifiers: [sunz_corrected, rayleigh_corrected]
- name: B04
modifiers: [sunz_corrected]
standard_name: toa_bidirectional_reflectance

airmass:
# PDF slides: https://www.eumetsat.int/website/home/News/ConferencesandEvents/DAT_2833302.html
# Under session 2 by Akihiro Shimizu (JMA)
Expand Down Expand Up @@ -204,7 +256,18 @@ composites:
prerequisites:
- name: B03
modifiers: [sunz_corrected, rayleigh_corrected]
- name: green
- name: hybrid_green
- name: B01
modifiers: [sunz_corrected, rayleigh_corrected]
high_resolution_band: red
standard_name: true_color

true_color_ndvi_green:
compositor: !!python/name:satpy.composites.SelfSharpenedRGB
prerequisites:
- name: B03
modifiers: [sunz_corrected, rayleigh_corrected]
- name: ndvi_hybrid_green
- name: B01
modifiers: [sunz_corrected, rayleigh_corrected]
high_resolution_band: red
Expand All @@ -223,7 +286,7 @@ composites:
compositor: !!python/name:satpy.composites.SelfSharpenedRGB
prerequisites:
- name: B03
- name: green_nocorr
- name: hybrid_green_nocorr
- name: B01
high_resolution_band: red
standard_name: true_color
Expand All @@ -235,7 +298,7 @@ composites:
prerequisites:
- name: B03
modifiers: [sunz_corrected, rayleigh_corrected]
- name: green_true_color_reproduction
- name: reproduced_green
- name: B01
modifiers: [sunz_corrected, rayleigh_corrected]
standard_name: true_color_reproduction
Expand Down