Keep or exclude total diffuse irradiance in the output when `return_components=True`

[cbcrespo]

As part of my GSoC 2026 project (#2750) I am working on diffuse irradiance models in pvlib to expand support for the return_components parameter, which is currently only implemented in some functions. When return_components=True, these functions return a dictionary containing the individual diffuse irradiance components (isotropic, circumsolar, horizon), but also their sum (total diffuse irradiance).

Several people (@AdamRJensen, @adriesse) have raised questions on whether it makes sense to keep the total or if it should be removed. Including the total is convenient and avoids requiring users to recompute it themselves if necessary. Excluding the total avoids redundant outputs and makes the meaning of return_components more explicit: return the components only.

The purpose of this issue is therefore to open a more focused discussion on this point and hopefully reach a consensus.

For reference, irradiance.get_total_irradiance and irradiance.poa_components currently also return several values that could be considered redundant: poa_global, poa_direct, poa_diffuse, poa_sky_diffuse, and poa_ground_diffuse. But in this case, there is no parameter controlling which values are returned.

[markcampanelli]

I'm going to make a bold opening move here with goal to "keep it simple stupid (KISS)" and suggest--

• Remove entirely the return_components input argument in all relevant functions in irradiance.py.
• Always return a dictionary with numpy.ndarray fields (possibly this is a 0-dimensional numpy.ndarray for the "scalar" case)--
  • A poa_sky_diffuse_total field.
  • A poa_sky_diffuse_components field, which is--
    • None if this function doesn't provide any components, OR
    • A sub-dictionary with the named POA diffuse components, when calculated by the function.

NOTE: Further investigate here to ensure "ease of use" by expected downstream consumers such as SAPM- vs. Marion-based effective irradiance computations for diffuse-component IAMs. Also, I mentioned elsewhere to consider we should try to indicate when, say, the circumsolar or horizon components are considered "concentrated", because that may affect how downstream consumers compute with the provided values. Including this info here seems to me like a stretch, but considering the downstream consumers NOW might shed some more light on if this would be worth the additional complexity.

YES, I am suggesting dropping the output of pandas objects and removing the input flag that changes the format of the output (KISS), but the functions could still accept other things than numpy.ndarray such as pandas.Series (AFAIK the timestamps are entirely optional for all current implementations). This follows the Robustness Principle, which says "be conservative in what you send, be liberal in what you accept". (Note, I do not view this as a hard and fast rule, but it seems applicable here given how the existing highly heterogeneous interfaces have evolved.)

NOTE: Arguments against always returning the components (when they are otherwise available from the algorithm) might include:

• They are somehow not always computed when computing the total.
• They hog up a lot of memory in some situations.
• ???

EXAMPLES:

perez (or perez_driesse) would return--

{
  "poa_sky_diffuse_total": numpy.ndarray(...),
  "poa_sky_diffuse_components": {
    "isotropic": numpy.ndarray(...),
    "circumsolar": numpy.ndarray(...),
     "horizon": numpy.ndarray(...),
  }
}

isotropic would return--

{
  "poa_sky_diffuse_total": numpy.ndarray(...),
  "poa_sky_diffuse_components": {
    "isotropic": numpy.ndarray(...),  # This would presumably equal poa_sky_diffuse_total.
  }
}

klucher would return--

{
  "poa_sky_diffuse_total": numpy.ndarray(...),
  "poa_sky_diffuse_components": None
}

Finally, because I hate naming things...

poa_sky_diffuse could simply be the name of the dictionary returned, with this string removed in the field names:

poa_sky_diffuse = {
    "total": numpy.ndarray(...),
    "components": {
      "isotropic": numpy.ndarray(...),
      "circumsolar": numpy.ndarray(...),
       "horizon": numpy.ndarray(...),
    }
}
return poa_sky_diffuse

Alternatively, a more explicit nested dictionary that stutters less:

return {
  "poa_sky_diffuse": {
    "total": numpy.ndarray(...),
    "components": {
      "isotropic": numpy.ndarray(...),
      "circumsolar": numpy.ndarray(...),
       "horizon": numpy.ndarray(...),
    }
  }
}

[wholmgren]

These functions are heavily used in PV modeling workflows (likely more so for perez, haydavies, and isotropic). We should be extremely cautious about making breaking API changes to these functions. So I would not remove keys that already exist.

The "always return a dict" approach has merit too but I am skeptical that it's worth the pain to users.

[AdamRJensen]

As Will pointed out, changing the existing functions is hard due to extensive usage and this being a breaking change. I personally am less hesitant in breaking existing functionality to be consistent, whereas Will is kinder to our users.

@cbcrespo For your future work, I guess it is more consistent to keep poa_sky_diffuse, since the existing functions are doing this (even though I don't agree with the practice, but that is now the standard).

[adriesse]

As Will pointed out, changing the existing functions is hard due to extensive usage and this being a breaking change. I personally am less hesitant in breaking existing functionality to be consistent, whereas Will is kinder to our users.

Claims about users' wishes are always hard to discuss because we don't have hard data. In this case, if someone is using the option return_components=True they must be using the names of those components to access their values. If they want the components, it appears highly unlikely that they are also accessing the "component" called poa_sky_diffuse; therefore, although removing poa_sky_diffuse would technically be an API change, it seems very unlikely to cause havoc in existing workflows.

[markcampanelli]

I note that my main workflow uses both the total POA irradiance (for cell temperature calculations) and the POA components (for effective irradiance calculations). I'm sure using the total POA irradiance (as opposed to effective POA irradiance) for cell temperature could be argued to be significantly wrong, but I view it as being conservative, and I think different models are either not explicit here or explicitly make different choices. In any event, I like having both outputs when available (esp. since components appear to commonly be computed as a byproduct of total), and the total is the invariant across all implementations.

Interestingly, the current Perez implementation appears to compute the total without separating the components first, then goes back and REcomputes the components when requested by return_components == True.

pvlib-python/pvlib/irradiance.py

Line 1211 in 0f7a205

# Calculate the different components

As an aside, in my application's default test case using a SAT, including the Marion IAM corrections on the diffuse components lowers generation by roughly 0.4%.

[ramaroesilva]

From @markcampanelli

Interestingly, the current Perez implementation appears to compute the total without separating the components first, then goes back and REcomputes the components when requested by return_components == True.

Similar redundancies have been reported by @adriesse in #2775 for reindl (link). As I mention there, this is implicitly resolved if it is decided that sky_diffuse is not provided with the components (not necessarily suggestion this solution, just stating).

[markcampanelli]

@ramaroesilva I have been focusing on pvlib.irradiance.perez right now (which I happen to be translating into Rust). I should be more careful in suggesting that "simply" summing the components to get the total is the solution here, because there are some edge cases to work around (in Perez).

There is a first edge case when relative airmass is NaN which appears to be handled in the existing algorithm.

A second edge case is when DNI and DHI are both zero in the computation of epsilon, which is apparently not currently handled. The offender appears to be in the computation of epsilon in the factor (dhi + dni) / dhi. Two questions/points:

1. Why is this not instead computed as 1 + dni / dhi ?
2. Physically, it seems to make sense to me that DNI -> 0^+ faster than DHI -> 0^+, so that it makes sense to define dni / dhi = 0 when both DNI and DHI are zero, and the whole factor should thus be unity. I personally find the NaNs returned in this edge case to be rather disruptive in my existing computational workflows. Perhaps @adriesse can comment on this, as he has spent a lot of time with this model?

Thus, one must presumably consider these edge cases appropriately when computing the components first, then the sum for the total should (hopefully) follow without issue.

[adriesse]

A second edge case is when DNI and DHI are both zero in the computation of epsilon, which is apparently not currently handled. The offender appears to be in the computation of epsilon in the factor (dhi + dni) / dhi. Two questions/points:

1. Why is this not instead computed as 1 + dni / dhi ?

I think this is done without thinking about it much, but I like it when the math in the code resembles the publications.

2. Physically, it seems to make sense to me that DNI -> 0^+ faster than DHI -> 0^+, so that it makes sense to define dni / dhi = 0 when both DNI and DHI are zero, and the whole factor should thus be unity. I personally find the NaNs returned in this edge case to be rather disruptive in my existing computational workflows. Perhaps @adriesse can comment on this, as he has spent a lot of time with this model?

I think we can safely say that in real life beam irradiance will disappear before diffuse, so it would be logical to set epsilon to 1 when both measurements are zero. This should only happen when the sun is well below the horizon. The calculation of zeta (replacement for epsilon) in the Perez-Driesse model handles this case explicitly.

But, we are digressing from the main question here...

[markcampanelli]

But, we are digressing from the main question here...

Apologizes if this continues to be out of scope, but I view this GSoC project as a good opportunity to tidy up things (esp. w.r.t. consistency of inputs and results) that might otherwise be considered a breaking change or disruptive bug fix.

@adriesse I very much appreciate your clarification here. I did notice (after writing the above) that in practice perez_driesse handles the potential 0/0 case more gracefully (via the dhi== 0 check in _calc_zeta, but perhaps dni==0 should be checked instead also?). Please know that my questions here are not asked with me "knowing the answer already". For example, the alternate form 1 + dni / dhi is what I would have had my calculus students do in order to try to reason about the limit as dni and dhi both go to zero from the positive side. In addition, this form would appear to be numerically optimal because it eliminates one variable reference. However, I have been surprised by the numerical _in_stability of such rewrites in the past! Hence my honest question.

I spotted the NaN epsilon issue with a very innocuous TMY file from NSRDB. FWIW, here's one of several offending time steps:

Year                     2008.000
Month                       2.000
Day                        19.000
Hour                       17.000
Minute                     30.000
temp_air                   10.000
Alpha                       1.270
Aerosol Optical Depth       0.130
Asymmetry                   0.640
dhi_clear                   0.000
dni_clear                   0.000
ghi_clear                   0.000
Cloud Fill Flag             7.000
Cloud Type                  7.000
temp_dew                    4.700
dhi                         0.000
dni                         0.000
Fill Flag                   7.000
ghi                         0.000
Ozone                       0.286
relative_humidity          69.730
solar_zenith               89.520
SSA                         0.940
albedo                      0.210
pressure                  931.000
precipitable_water          1.800
wind_direction            252.000
wind_speed                  1.600
Name: 2008-02-19 17:30:00-08:00, dtype: float64

Chasing down the issue now suggests another question: Is pvlib responsible for checking the (apparent) assumption that whenever DHI is zero, then DNI should likewise be zero? (And which weather providers might be guaranteeing this?)

Similar "quality" questions arise when, say, a negative diffuse horizontal component is so big that the total diffuse is negative. (Note the line sky_diffuse = np.maximum(dhi * (term1 + term2 + term3), 0) here.) Does a source reference indicate what to do here with this non-physical total irradiance, or must this be the fault of a non-physical weather file (as opposed to the algorithm)? (The docstring and later single comment # Set values of components to 0 when sky_diffuse is 0 leaves a bit to be desired here, IMHO.)

[markcampanelli]

I put my thoughts into action, if anyone would care to take a look at #2808.

[adriesse]

Chasing down the issue now suggests another question: Is pvlib responsible for checking the (apparent) assumption that whenever DHI is zero, then DNI should likewise be zero? (And which weather providers might be guaranteeing this?)

pvlib is free and has no responsibilities! More seriously though, negative solar irradiance doesn't exist and I don't know of any model in pvlib that is validated for negative irradiance.

But back to the main question: Keep or exclude total diffuse irradiance in the output when return_components?