Feature request for extending the Rayleigh reduction functionality

[strandgren]

Background

The built-in method (reduce_rayleigh_highzenith)](https://github.com/pytroll/pyspectral/blob/main/pyspectral/rayleigh.py#L287) offers a way to linearly reduce the Rayleigh correction at high sun zenith angles in order to compensate for the fact that RTM models and thus the Rayleigh correction computation breaks down at large angles. This reduction Rayleigh correction can be defined/tuned using three parameters:

• thres_zen: the lower sun zenith angle limit, i.e. where to start reducing the Rayleigh correction.
• maxzen: the upper sun zenith angle limit, typically where to the Rayleigh correction is completely removed but can be ambiguous (see below).
• strength: the aggressiveness of the of the Rayleigh correction reduction between thres_zen and maxzen, i.e. the steepness of the linear Rayleigh reduction function that starts decreasing at thres_zen.

Using strength=1.0 will result in full Rayleigh correction at thres_zen and zero Rayleigh correction at maxzen (and beyond), with a linear decrease in between. A strength below 1.0 will reduce the slope of the linear decrease and thus shift the point of zero Rayleigh correction to a larger angle.

New feature proposal

To make the transition from full to no Rayleigh correction smoother I suggest to use a smoothstep function for the Rayleigh reduction rather than the current linear function. This could better capture the non-linear behavior at large angles. Below is an example of a smoothstep function where the strength is kept and still controls the steepness of the reduction between thres_zen and maxzen:

Note that thres_zen and maxzen always define the points of full and zero Rayleigh correction, respectively, independent of strength, as opposed the the current linear scaling function. Furthermore, the the smoothstep function used here simplifies to a linear function if strength=1.0, hence preserving the current functionality and backwards compatibility.

Below are two examples of an AHI true color RGB with the current linear reduction function (left, current main) and a slight non-linear reduction function (right, proposed feature):

main: thresh_zen=70, maxzen=105, strength=1.0	Proposed feature: thresh_zen=70, maxzen=105, strength=1.5

Although the differences may not strike as significant, the non-linear reduction function clearly results in more vivid colors in eastern Australia and Monglia as a result of slower reduction in the Rayleigh correction. This is a single example and there might be more optimized tuning settings for the current linear implementation. However, this is also true for the non-linear example and the the proposed non-linear function including the tunable steepness parameter strength in combination with different combinations of thres_zen and maxzen should offer even more extensive tuning and optimization possibilities of the Rayleigh correction reduction.

If this is considered a useful feature, I would be happy to make a PR since I have most of the the code ready.

Open question

The main question is how to integrate this into the current main with minimum changes for other users and full backwards compatibility.

Actually I would be in favor of removing the current functionality of the strength parameter which simply changes the steepness of the linear reduction function and increases the sun zenith angle at with the Rayleigh correction is completely removed. For any value of strength less than 1.0 it's not transparent to the user at which sun zenith angle the Rayleigh correction reaches zero. I would rather be in favor of keeping the linear strength constant (=1.0) and let the user modify maxzen accordingly instead, which would give the user full transparency since maxzen would always define the point where the Rayleigh correction is completely removed. The figure below shows how the linear reduction function for constant thres_zen/maxzen and variable strength can be replicated simply using variable maxzen and constant strength (=1.0)):

Suggested solution

Given the fact that:

1. The current strength is only valid for the range [0.0, 1.0], and
2. The new proposed non-linear strength is only valid for the range [1.0, inf)

I suggest the following:

Short-term implementation (following PR):

• For strength < 1.0 compute the Rayleigh reduction following the current linear implementation.
• For strength >= 1.0 compute the Rayleigh reduction following the non-linear smoothstep function proposed above.
• Add a deprecation warning that strength < 1.0 will loose support and that users shall modify maxzen instead to weaken the "strength" of the linear decrease.

Long-term implementation (a couple of releases after suggested PR is released):

• For strength < 1.0: Remove deprecation warning and throw error that this is an invalid strength.

One could also implement the proposed non-linear function independently, but given the fact that any change to current strength can be achieved by increasing maxzen, I think the above would be the most clean solution. But if preferred, one could also skip the proposed deprecation warning and long-term implementation and simply have different behavior for strength=[0.0, 1.0) (current functionality with linear scaling) and strength=[1.0, inf) (new proposed non-linear scaling).

---

Before making a PR I would be happy to hear any feedback on this, especially from e.g. @simonrp84 , @mherbertson and @yukaribbba who have implemented and/or shown interest in the current Rayleigh reduction functionality :)

[mherbertson]

This is very interesting! As this follows a non linear approach to correction I think it would have fixed issues I ran into when trying to get the overall balance between strength, thresh and max zen. I've now moved onto the JMA method of blending the uncorrected image so unfortunately wont be able to provide much technical assistance.