Pvlib-python: Bug in extraterrestrial radiation calculation

[wholmgren]

The ET radiation calculator needs the "day angle" of the Earth's position in its orbit around the Sun. Our extraradiation function and the matlab pvl_extraradiation calculate the angle as:

B = (2. * pi / 365.) * doy

Assuming that doy is a value between 1 and 365, I believe it should be

B = (2. * pi / 365.) * (doy - 1)

The disc and pvl_disc functions use the (doy-1) formula.

[wholmgren]

Here are some plots to show the difference.

doy

doy - 1

doy

doy - 1

Some quick pandas.describe() statistics for the difference between the spencer method and the pyephem methods:

doy

mean       0.003336
std        1.260789
min       -1.855202
25%       -1.184669
50%       -0.065613
75%        1.233035
max        1.913360

doy - 1:

mean       0.003264
std        0.935147
min       -1.390897
25%       -0.870512
50%       -0.055248
75%        0.963799
max        1.356251

doy - 1 wins, as expected.

[cwhanse]

Will,

Thanks catching this.

What do you think we should do with a leap year?

B = 2 pi / max(DOY) * (DOY – 1) perhaps?

Cliff

[wholmgren]

I'd say that we should follow the literature that these methods come from, which, if I recall correctly, means ignoring leap years. Also, I think that max(DOY) would only work if DOY is an array with a full year of data.

On a related note, I wrote some code to pull the Earth-Sun distance from pvlib pythons implementation of the NREL SPA. It's about 50x slower than the Spencer method but 10x faster than using PyEphem. I'll post it later today.

[dacoex]

A general question:
When such kind of flaws are discovered, the issue and its resolution may change results for a certain location and data set.

For the point of testing but als tracebaility of results obtained with different versions, would it make sense to have a accuracy & uncertainty log file?
We could have a few BSRN locations which we use as reference sites along with test data sets.
For each version, we record in a CSV, HDF5 or alike the results for a typical calculations chain.

With this, it will be possible to demonstrate:

• Improvement of the library over time not only in code quality & new functionality but also accuracy
• Level of inaccuracy of older output, e.g. when version 0.2 was used to calculate the results 1 yr. ago and now a re-run shows different values.

Just as an example, it is common in industry for important projects to fix a certain PVSyst version and have several versions of this program in parallel to avoid inconsistencies on the software level when comparing different configurations of components.

[adriesse]

If one method takes into account leap years and the other doesn't, the differences will look different from year to year. How do the stats and plots look for a four-year period with a leap year?

[wholmgren]

@adriesse there's not much difference. Here's 2015 and 2016 at hourly resolution. Note that the wiggles in the difference plot are due to the fact that the extraradiation function converts a DatetimeIndex into integer days.

In [36]: spencer_pyephem['2015'].std()
Out[36]: 1.1428934417419712

In [37]: spencer_pyephem['2016'].std()
Out[37]: 1.2049695065475967

In [38]: spencer_doy1_pyephem['2015'].std()
Out[38]: 0.8636022510479551

In [39]: spencer_doy1_pyephem['2016'].std()
Out[39]: 0.8812626386839624

Some perspective... The differences between the spencer(doy), spencer(doy-1), and pyephem methods shown above are on the order of 1 to 2 W/m^2. The solar "constant" varies by a few W/m^2 over the 11 year solar cycle. The accepted value of the solar constant has also been revised downward by a few percent over the last 30 years. Changing pvlib's default solar constant would have a 3-5x larger impact on model results than changing doy to doy - 1. If you want to read more, the solar constant wikipedia page and CU's LASP TIM page are reasonable places to start.

So, while I think we should change doy to doy - 1 for the sake of mathematical consistency and literature consistency, I do not think that it makes a meaningful difference for PV modeling. For example, PR #215 demonstrates that the PV power test results are consistent to better than 1 part in 1000.

@dacoex I'm not sure what you're asking for that's not already covered by the whats new documentation or the unit tests. Feel free to make a pull request with a concrete suggestion.

[wholmgren]

I'm closing this via PR #215, though feel free to continue to discuss here or in another issue if you feel otherwise.

[Schmelzersolar]

While looking at the code for the calculations of extraterrestrial radiation get_extra_radiation() and crosschecking it with literature, I noticed that for both methods 'asce' and 'spencer' the solarposition._calculate_simple_day_angle(doy) function is used. This function basically just calculates (2. * np.pi / 365.) * (doy- 1), which is correct for method='spencer'.

But I think for method='asce' it should just be (2. * np.pi * doy / 365.) - without "- 1" - as described on page 9 in "J. A. Duffie and W. A. Beckman, "Solar Engineering of Thermal Processes, 3rd Edition" J. Wiley and Sons, New York (2006)"

--> Before this bug fix it was wrong for method='spencer', now it's wrong for method='asce'?

[wholmgren]

@Schmelzersolar thanks for the report. I'll reopen the issue. Pull request welcome.

[wholmgren]

@cwhanse can you provide any guidance on this?

[cwhanse]

@wholmgren yes but not today, I think we'll need to inspect the source reference (ACSE publication) and I don't have a copy handy.

[mikofski]

@wholmgren Nice catch! I believe this function already exists in solarposition._calculate_simple_day_angle():

def _calculate_simple_day_angle(dayofyear):
    """
    Calculates the day angle for the Earth's orbit around the Sun.

    Parameters
    ----------
    dayofyear : numeric

    Returns
    -------
    day_angle : numeric
    """
    return (2. * np.pi / 365.) * (dayofyear - 1)

It's used in the equation_of_time_spencer71, equation_of_time_pvcdrom, declination_cooper69, and declination_spencer71. Perhaps extraterrestrial radiation can use the existing function?

[cwhanse]

@Schmelzersolar is correct. To faithfully implement the extraterrestrial radiation model in Spencer 1971, the day angle is (2. * np.pi / 365.) * (dayofyear - 1) where dayofyear runs from 1 to 365 (366 in a leap year). However, in the ASCE report the day angle is (2. * np.pi / 365.) * (dayofyear). We should reference the ASCE document:

ASCE (2005), THE ASCE STANDARDIZED REFERENCE EVAPOTRANSPIRATION EQUATION, Environmental and Water Resources Institute of the American Civil Engineers, Ed. R. G. Allen et al.

[cwhanse]

A reminder of the original motivation for pvlib - to provide reference implementations of models from literature.