Skip to content

Refactor code for solar zenith angle correction modifiers #2955

Description

@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:

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:

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

Metadata

Metadata

Assignees

Labels

Type

No type

Projects

No projects

Milestone

No milestone

Relationships

None yet

Development

No branches or pull requests

Issue actions