Refactor code for Sun zenith angle corrections and change to effective_solar_pathlength_corrected instead of sunz_corrected for built-in RGB recipes

[strandgren]

This PR aims to clarify the code used for the Sun zenith angle corrections in satpy, including improving documentation and log messages to make it more transparent to the user. Furthermore, the code is prepared to eventually change the default behavior of the SunZenithCorrector/sunz_corrected modifier to remove the default reduction at high angles in order to compute the true reflectance. This change is planned in Satpy v1.0 and can be done by changing the new configuration parameter introduced in _config.py here to "use_legacy_sunz_correction": False. Finally, the built-in recipes are modified to improve imagery and comply with the recommendation from the WMO RGB workshop in Norrkoeping in April 2025.

Background
In satpy we have two corrections for the Solar zenith angle/Solar path length available, the standard SunZenithCorrector which applies the 1/cos(sunz) correction and the EffectiveSolarPathLengthCorrector which applies the parameterization proposed by Li and Shibata (2006).

There is also a default reduction of the correction starting at a given Sun zenith angle (88 degrees). This is intended for (RGB) imagery to avoid overcorrection at very large angles. However, this is not desirable for quantitative or scientific use of the data, where the standard 1/cos(sunz) should be used as is in order to compute the true reflectance. This functionality for capping/reduction has also been implemented in the EffectiveSolarPathLengthCorrector (duplicate code), but as shown in #3096 this reduction doesn't have any added value for the Li and Shibata parameterization, which works best as is, since it already accounts for the overcorrection at higher angles. For this reason, it was also agreed at the WMO RGB workshop in Spring 2025, that the Li and Shibata parameterization should be the recommended method for normalizing the data for the Solar zenith angle (without any capping or reduction).

Furthermore, some datasets come with the 1/cos(sunz) already applied, in which case the sunz_corrected modifier is attached to the dataset. However, this is not accurate since the default satpy sunz_corrected modifier is not only the simple 1/cos(sunz), but also includes the reduction of the correction at higher Solar zenith angels. Therefore it would also be good if the default sunz_corrected would rather be the simple 1/cos(sunz) corrected, which would then be in-line with these datasets and also allow users in general to compute the true reflectance - this becomes even more relevant with #3292 where we clearly separate between "reflectance" data normalized by the Solar zenith angel and not.

---

Changes in this PR
In response to the points above this PR includes the following changes:

Python code:

• Remove the capping/reduction of the EffectiveSolarPathLengthCorrector modifier. The parameters for this are still kept as valid input to avoid run-time errors if provided, but have no effect. A warning is issues if provided.
• Move the underlying code for EffectiveSolarPathLengthCorrector to the same place as the other Solar zenith angle correction methods (./satpy/satpy/modifiers/angles.py)
• Refactor SunZenithCorrector and underlying methods to make the functionality more clear.
• Change the default behavior of SunZenithCorrector to compute the true reflectance without any reduction of the the correction. Kept as is for now, to be changed for Satpy v1.0.
• Add the configuration parameter use_legacy_sunz_correction which can be used to control whether to use the current legacy behavior of SunZenithCorrector (with reduction) or not (simple 1/cos(sunz) correction) if no reduction limits are explicitly defined in the modifier yaml configuration. For now this is set to True and should be changed to False in Satpy v1.0. Corresponding RST documentation has been added.
• Clarify documentation for SunZenithCorrector and EffectiveSolarPathLengthCorrector, also highlighting the differences between the two and clarify that the former can be used to compute the true reflectance, whereas the latter is tailored for imagery.
• Modify methods to be dask-compatible and remove wrapping using map_blocks
• Add more log-messages and warnings depending on the use of the two corrections.
• Extend test coverage, especially since some parts were previously not tested

yaml-recipes

• Change pre-configured (RGB) imagery recipes to use the effective_solar_pathlength_corrected modifier instead of sunz_corrected for the solar zenith angle correction. As demonstrated in Standardize the correction for atmospheric path length for solar channels #3096 this leads to better imagery and was also suggested as recommendation at WMO RGB workshop in spring 2025.
• Change pre-configured Rayleigh correction modifiers for (RGB) imagery recipes to use the effective_solar_pathlength_corrected modifier instead of sunz_corrected for the solar zenith angle correction of the red band used for reduced correction over clouds. This is done since sunz_corrected is changing to become the simple 1/cos(sunz) correction which would not work well for imagery purposes since the reflectance of the red band would become 0 for angles larger than 90 degrees and thus break the Rayleigh correction. The Rayleigh correction reduction over clouds would also not work well close to 90 degrees, since also clear-sky becomes very bright at high sun zenith angles with the 1/cos(sunz) correction.
• Remove sunz_reduced modifier for FCI true_color, cloud_phase and cloud_type RGBs to be in-line with other sensors and recommendation from RGB workshop. Additional true_color_sunz_reduced composite is added to still support the generation of this RGB with a smoother transition into deep-space. There may be some further changes/revisions to this since I plan some further optimization of the imagery with focus on the twilight area and blending with night-time imagery in a later PR once this one is approved and merged.

---

Special cases
Some readers provide the reflectance data with the 1/cos(sunz) correction already applied, in which case the sunz_corrected modifier is already attached to the dataset. This is the case for the following readers:

• meris_nc_sen3
• msi_safe_l2a
• olci_l2
• viirs_sdr
• clavrx

For meris_nc_sen3, olci_l2 and clavrx there are no pre-configured RGB recipes, so the changes in this PR are not relevant for these data. For msi_safe_l2a the data have the dedicated modifier esa_sunz_corrected attached to them, which is also used for the corresponding RGB recipes. Hence, the changes in the PR are not relevant for these data either.

The only relevant one is VIIRS (SDR). Since the VIIRS SDR data come with the Solar zenith angle correction already applied, the dataset gets the sunz_corrected modifier attached. However, VIIRS L1B data come without the 1/cos(sunz) correction and therefore no sunz_corrected modifier is attached to the dataset. For this reason, we can use the EffectiveSolarPathLengthCorrector modifier for the recipes, but we have to name it sunz_corrected, in order for satpy to understand that this is equivalent to the sunz correction already in the SDR data in the sense that no further correction should be applied.

A better and more generic solution would be to modify the dependency tree to know that if the user requests a dataset with the effective_solar_pathlength_corrected modifier, but only finds a DataID with the modifier sunz_corrected, it would be regarded as equivalent and a substitute for the effective_solar_pathlength_corrected modifier. However, such a solution is outside the scope of this PR.

After the monthly meeting on June 2, it was agreed to use the regular effective_solar_pathlength_corrected modifier also for VIIRS data. Currently this won't work with VIIRS SDR data, since satpy will fail to find the "base" version of the requested without modifiers. It was agreed that this should be solved by adding the unnormalized_reflectance calibration level to the viirs_sdr reader once #3292 is merged, see details in #3411. We then just need to make sure that the VIIRS composites specify the expected calibration level in the composite and/or modifier recipes.

---

Changes for users and backwards incompatibility

• For the built-in satpy recipes the changes are generally minor and to the better with smoother imagery close the the terminator. For example the bright line along the terminator is no longer there which was caused by the very sharp peak of the sunz_corrected modifier at 88 degrees. See Standardize the correction for atmospheric path length for solar channels #3096 and comments in the PR for examples.
• For FCI users, the removal of the sunz_reduced modifier leads to significantly more data being retained/visible close to the terminator. see examples below.
• The reduction of the effective_solar_pathlength_corrected modifier has been removed as it's not needed, given that the Li and Shibata parameterization already deals with the over-correction at high solar zenith angels. If the reduction parameters are still provided they will be accepted but ignored and a warning is raised.
• VIIRS composites using the effective_solar_pathlength_corrected modifier will not be possible with VIIRS SDR data until Add unnormalized_reflectance to VIIRS SDR (viirs_sdr) reader #3411 has been fixed.
• The default reduction of the sunz_corrected has been disabled (but still with support to apply it if needed). Hence, the sunz_corrected modifier will compute the true reflectance if applied with it's new default configuration. This will impact local composites which use the upstream sunz_corrected modifiers defined in satpy (except for VIIRS which is treated differently given the differences between the L1B and SDR data). The result will be very bright imagery close to 90 degrees and no data beyond 90 degrees. A warning has been added when using these modifiers to inform users on this change and how to fix it.
• Similarly, if users have defined there own local sunz_corrected modifier with default configuration (i.e. no correction_limit or max_sza), the imagery will change close to the terminator (again becoming very bright close to 90 degrees and black beyond).

Following the monthly meeting on June 2, this PR has been changed such that there are no changes to the default values in SunZenithCorrector or corresponding pre-configured sunz_corrected modifiers in satpy at this point. Instead we will change this in v1.0 (#3412) and a warning is therefore now raised about this upcoming change of the default values of SunZenithCorrector are used. A new configuration parameters use_legacy_sunz_correction has been added to control whether to use the legacy correction (with reduction) or not. Currently it's set to True, but should be changed to False in Satpy v1.0

---

• Closes Refactor code for solar zenith angle correction modifiers #2955, Standardize the correction for atmospheric path length for solar channels #3096
• Tests added
• Fully documented

---

AI was used to support this PR, mainly for refactoring tests and writing doc-strings

[djhoese]

In that case, we could change all composites in Satpy to use EffectiveSolarPathLengthCorrector instead (already proposed in #3096). However, that would lead to backward incompatibility and change of behavior for users having defined SunZenithCorrector with default configuration in the local configurations (local configurations specifying correction_limit and max_sza would still work as before).

I think in general I'm not against this PR "in spirit". The one thing I had trouble with when trying to switch to effective path length as described in your point above is that for some readers the / cos(SZA) is already applied as part of the file creation. For example, VIIRS SDRs. So for most instrument-specific composite definitions we could define RGBs with a the effective path length modifier as part of the prerequisites and users could request the "sunz_corrected" if they wanted to "do science".

However, I think this means that for VIIRS composites we'd have to redefine the sunz_corrected modifier to be the effective path length and use it for all RGBs. This way VIIRS SDRs don't get any extra correction and the RGBs can be made, but then for VIIRS L1Bs (where no / cos(SZA) is pre-applied) they would get the effective solar pathlength correction and then RGBs could be made. This does mean RGBs would look different between VIIRS SDR and L1B inputs though.

[simon-sat]

You could always reverse the cos(sza) and then apply the effective path length...

[djhoese]

🤔 I'd prefer not to.

[strandgren]

@djhoese I understand the concern, but if I understand the comment correctly I don't think the proposed changes would be a problem and no significant change to the current behavior.

Now the default sunz_corrected modified already reduces the correction between 88-95 degrees in order to avoid over-correction and get better looking imagery. Hence, the standard sunz_corrected does not compute the real reflectance at high sunz angles and tus dehaves differently compared to the correction already applied in the VIIRS SDR data. This reduction is a rather crude alternative to the Li and Shibata parameterization, which I tried to demonstrate in #3096. What this means:

1. If you derive RGBs from VIIRS SRD data with the pure 1/cos(sunz) correction already applied in the data, they will already look different compared to RGBs derived from VIIRS L1B data where the /cos(sunz) correction including default reduction is applied in satpy.
2. If preferred, you will still be able to preserve the current default behavior of sunz_corrected by defining it like this:

  sunz_corrected_legacy:
    modifier: !!python/name:satpy.modifiers.SunZenithCorrector
    correction_limit: 88
    max_sza: 95
    prerequisites:
    - name: solar_zenith_angle
      resolution: 742

Does that make sense and sound ok for you you, or did I miss something else?

I'll try to implement my suggested changes today such that it's more clear what I propose, and then we can see i the review process if this is acceptable or if I should revert something.

[djhoese]

What you say makes sense but I can't tell if you're saying we would still need to do what I mentioned. That is, VIIRS composites would need to use "sunz_corrected" for RGBs but either have to:

1. Stick with regular sunz correction with or without caps.
2. Switch to path length to match the rest of satpy but still call the modifier "sunz_corrected" so that VIIRS SDRs could still be used for RGBs.

[strandgren]

@djhoese Thanks a lot for the review and your comments! I have addressed them all now

[djhoese]

Is it possible to move this check to where the map_blocks call is so the error happens immediately rather than at compute time?

[djhoese]

Oh but also, these debug messages (logger.debug) will show up for every chunk that gets processed 🤔 What's the best way to do this...

[djhoese]

I mean...would it be terrible to move the logging and error of the functions up to just before the map_blocks call? Oh and there's the debug message for the else: case below.

[strandgren]

uups.. lesson learned not to put log messages in map_blocks, quite some spamming in the debug output indeed.

I moved the log messages and the ValueError to the parent method where I now also define the method and pass that as an argument to the map_blocks call - I think this looks cleaner overall. I still had to keep a check for a possible ValuesError in _sunzen_corr_cos_ndarray to satisfy mypy, but given how sunzen_corr_cos is defined, we would never reach that ValuesError.

I also moved the log message in _atmospheric_path_length_correction_ndarray for the same reason.

[djhoese]

Sorry this is turning into such a headache, but Claude points out that the dtype in the map_blocks call is actually wrong. Rather than summarize what it said, here is the exact content:

_atmospheric_path_length_correction_ndarray (and _sunzen_reduction_ndarray) can return a
 dtype other than the one declared to map_blocks. _sunzen_corr_cos_ndarray casts
 (angles.py:628), but the other two do not. map_blocks now declares dtype=data.dtype /
 meta=np.array((), dtype=data.dtype), so if cos_zen is float64 while data is float32 the dask
 array claims float32 and produces float64 chunks. Verified locally:

 declared float32, actual float64

 That combination is reachable: get_cos_sza() casts lons/lats to the data dtype, but the
 SZA-provided branch does coszen = np.cos(np.deg2rad(projectables[1])) with whatever dtype the
 reader gave. The tests never hit it because call_sunz_modifier casts both arrays to the same
 dtype. One-line fix in each: corr = corr.astype(data.dtype, copy=False) before the multiply, and
 ideally a test with mismatched data/SZA dtypes.

[strandgren]

good catch. I just added one more level of parameterization to test different dtype combinations for data_arr and sunz_sza making sure that the dtype of data_arr is preserved . As expected this lead to test failures for _atmospheric_path_length_correction_ndarray and _sunzen_reduction_ndarray, but after adding corr = corr.astype(data.dtype, copy=False) to these methods, all tests now pass.

I did keep the dtype=dtype in the map_blocks call, which I think should be the case now when the dtype issue is fixed.

[djhoese]

You should be able to check this in the test and is something we have in tests. You can check the .dtype of the dask array returned and the .dtype of the computed numpy array created from that dask array. If they are the same and are np.float32 then 🎉

[strandgren]

I tried to satisfy codescene, but for the remaining complaints I don't really agree:

• sunzen_corr_cos - I see no way of reducing the "complexity" without making the code less readable
• SunZenithCorrector._apply_correction Since this concerns checks only needed for backwards compatibility towards v1.0 and will be removed later, I see no big need changing this (since again I think it's clear and readable now)
• call_sunz_modifier - One argument too much, please accept this😄

[djhoese]

I think sunzen_corr_cos could be refactored into "standard" being at the top then an else that calls a helper function to determine some other method...but just don't do it. I agree with you on the code being less readable and I agree with you on the others.

Now if other things are all cleaned up maybe I should get to reviewing the overall idea again. What is the last thing we decided on VIIRS? Wait on the unnormalized_reflectance calibration level?

[djhoese]

You mentioned in the other thread that this couldn't be moved. What if as a first try to make mypy happy you:

1. Moved this to the caller.
2. Change max_sza here to not be optional.
3. In the method if statements set correction_limit to 0.0 and max_sza to 0.0 if they are None. Oh! Or set them to NaN.

Is that going too far? Either way I think correction_limits type annotation here is wrong since it says float and no optional.