Refactor code for solar zenith angle correction modifiers

[strandgren]

Feature Request

Is your feature request related to a problem? Please describe.
Satpy includes two methods for solar zenith angle correction of VIS/NIR data:

• standard 1/cos_zen correction: https://github.com/pytroll/satpy/blob/main/satpy/modifiers/angles.py#L534
• effective path length correction following Li and Shibata: https://github.com/pytroll/satpy/blob/main/satpy/utils.py#L270

The overall code in these methods is in principle identical with only minor differences, namely the computation of the correction corr:

• 1/cos_zen vs 24.35 / (2. * cos_zen + np.sqrt(498.5225 * cos_zen**2 + 1))

. The same applies to the corresponding modifier classes:

• https://github.com/pytroll/satpy/blob/main/satpy/modifiers/geometry.py#L77
• https://github.com/pytroll/satpy/blob/main/satpy/modifiers/geometry.py#L120

My impression is that having two modifiers with separate code that do very similar things can lead to confusion regarding the actual differences. Hence I would propose some refactoring work

Describe the solution you'd like

Alternative 1:
As a first alternative I propose to do a "light" refactoring of the the actual code for the sunz correction and merge the methods used for the sunz correction and instead pass a method keyword, e.g. standard or li_shibata.

It could look something like this:

def _sunzen_corr_ndarray(data: np.ndarray,
                             cos_zen: np.ndarray,
                             method: str,
                             limit: float,
                             max_sza: Optional[float],r) -> np.ndarray:
    # Convert the zenith angle limit to cosine of zenith angle
    limit_rad = np.deg2rad(limit)
    limit_cos = np.cos(limit_rad)
    max_sza_rad = np.deg2rad(max_sza) if max_sza is not None else max_sza

    # Cosine correction
    if method == "standard":
        corr = (1. / cos_zen).astype(data.dtype, copy=False):
        corr_lim = (1. / limit_cos).astype(data.dtype, copy=False):
    elif method == "li_shibata":
        corr = get_sunz_corr_li_and_shibata(cos_zen)
        corr_lim = get_sunz_corr_li_and_shibata(limit_cos)
    else:
        raise error

    if max_sza is not None:
        # gradually fall off for larger zenith angle
        grad_factor = (np.arccos(cos_zen) - limit_rad) / (max_sza_rad - limit_rad)
        # invert the factor so maximum correction is done at `limit` and falls off later
        with np.errstate(invalid="ignore"):  # we expect space pixels to be invalid
            grad_factor = 1. - np.log(grad_factor + 1) / np.log(2)
        # make sure we don't make anything negative
        grad_factor = grad_factor.clip(0.)
    else:
        # Use constant value (the limit) for larger zenith angles
        grad_factor = 1.
    corr = np.where(
        cos_zen > limit_cos,
        corr,
        (grad_factor * corr_lim).astype(data.dtype, copy=False)
    )
    # Force "night" pixels to 0 (where SZA is invalid)
    corr[np.isnan(cos_zen)] = 0
    return data * corr

Alternative 2
The second alternative, and my preference, is to do the above but also remove the EffectiveSolarPathLengthCorrector modifier class completely and instead pass the method keyword already in SunZenithCorrector. This would in my opinion be the cleanest and most transparent solution, but we'd need to have a deprecation warning for a few of releases.

Describe any changes to existing user workflow

Alternative 1: None.

Alternative 2:
modifiers using the EffectiveSolarPathLengthCorrector class would have to change from

  effective_solar_pathlength_corrected:
    modifier: !!python/name:satpy.modifiers.EffectiveSolarPathLengthCorrector
    optional_prerequisites:
      - solar_zenith_angle

to

  effective_solar_pathlength_corrected:
    modifier: !!python/name:satpy.modifiers.SunZenithCorrector
    method: "li_shibata"
    optional_prerequisites:
      - solar_zenith_angle

Additional context
Furthermore the reduction functionality of both sunz correction methods has similarities with the reduction applied in the SunZenithReducer modifier here: https://github.com/pytroll/satpy/blob/main/satpy/modifiers/angles.py#L576. Hence, I would also consider refactoring this part of the code and perhaps integrate the SunZenithReducer modifier in the refactored SunZenithCorrector class (TBD depending on final readability and transparency of the code).

[strandgren]

I would be happy to work on this at some point, but it won't be immediately. Hence, I'm happy to get feedback from others before I would start.

[strandgren]

For some context and future reference, this is how the two sunz correction methods work, with and without the default reduction above 88 degrees:

I think the sharp transition at 88 degrees is the reason why we sometimes see a bright line along the terminator. Hence, we could also consider making the transition smoother, for example using the strength functionality in SunZenithReducer, which would anyway available if including that modifier in the refactoring work.

[mraspaud]

I like the idea of the refactoring! I'm just wondering how that will affect the DataIDs we use in satpy, where we now describe the modifications to a band with a list of strings...

[strandgren]

@mraspaud you mean when we do something like this: scn.load(['B03'], modifiers=('sunz_corrected', 'rayleigh_corrected'))?

in that case I don't think anything would change since my understanding is that these strings point to predefined modifiers in e.g. etc/composites/visir.yaml? And these names would not change, only the actual modifier class used for the pre-defined effective_solar_pathlength_corrected modifier (as shown in the issue description). Or do I miss(understand) something else?

[mraspaud]

@strandgren you didn't missunderstand, it was exactly what I was wondering. Thanks for refreshing my memory, I see there won't be any problem on that side then.

[djhoese]

Seems reasonable to me. Option 2 makes sense.