Add diffuse IAM integration and gallery example

[kandersolar]

• Closes Add Bill Marion's diffuse IAM integration method? #983
• I am familiar with the contributing guidelines
• Tests added
• Updates entries to docs/sphinx/source/api.rst for API changes.
• Adds description and name entries in the appropriate "what's new" file in docs/sphinx/source/whatsnew for all changes. Includes link to the GitHub Issue with :issue:`num` or this Pull Request with :pull:`num`. Includes contributor name and/or GitHub username (link with :ghuser:`user`).
• New code is fully documented. Includes numpydoc compliant docstrings, examples, and comments where necessary.
• Pull request is nearly complete and ready for detailed review.
• Maintainer: Appropriate GitHub Labels and Milestone are assigned to the Pull Request and linked Issue.

I know I didn't wait for feedback on #983, I just wanted to implement it for fun and figured as long as the code exists, might as well open a PR. This vectorized implementation is a couple orders of magnitude faster than a straightforward translation of the paper's pseudocode. The gallery example recreates the figures in the paper.

API reference page: https://pvlib-python--984.org.readthedocs.build/en/984/generated/pvlib.iam.marion_integrate.html

Gallery example: https://pvlib-python--984.org.readthedocs.build/en/984/auto_examples/plot_diffuse_aoi_correction.html

[cwhanse]

This addition of a means to calculate the diffuse IAM using the various IAM models is valuable. I would make the marion_integrate function more along the lines of pvlib.iam.martin_ruiz_diffuse in signature.

[cwhanse]

LGTM, thanks @kanderso-nrel

[cwhanse]

@wholmgren @mikofski any comments on this PR? I think it is ready to merge.

[wholmgren]

I focused my review on the API and did not check the implementation.

I'd be in favor of using a scipy integration method over a strict implementation of the reference. Is that feasible?

[kandersolar]

Thanks @wholmgren for this detailed review. I'm happy to make the suggested changes, although I'm not sure about this one:

I'd be in favor of using a scipy integration method over a strict implementation of the reference. Is that feasible?

Here's a quick comparison for the 'sky' case:

from pvlib import tools, irradiance, iam
from scipy import integrate
from functools import partial

def scipy_integrate(function, surface_tilt):

    def integrand(phi, theta, beta, function):
        aoi = irradiance.aoi(beta, 0, phi, theta)
        if aoi >= 90:
            return 0
        area_element = tools.sind(phi)
        return function(aoi) * tools.cosd(aoi) * area_element

    numerator = integrate.dblquad(integrand,
                                  0, 360,  # azimuth
                                  lambda theta: 0, lambda theta: 90,  # zenith
                                  args=(surface_tilt, function))
    denominator = integrate.dblquad(integrand,
                                    0, 360,  # azimuth
                                    lambda theta: 0, lambda theta: 90,  # zenith
                                    args=(surface_tilt, lambda *args: 1))
    Fd = numerator[0] / denominator[0]
    return Fd

iam_model = partial(iam.physical, n=1.3)

%time scipy_iam = scipy_integrate(iam_model, 20)
Wall time: 1.24 s

%time marion_iam = iam.marion_integrate(iam_model, 20, 'sky')
Wall time: 6.98 ms

print(scipy_iam, marion_iam)
0.9622497940399887 0.9622503413733542

So they get the same answer, but this simple scipy version is much slower, presumably for lack of vectorization. I'm not very familiar with scipy so maybe it's possible to nest two integrate.quad_vec calls. I'm also not sure how to scale it across an array of tilts other than a list comprehension or np.vectorize. What do you think?

[wholmgren]

Thanks for exploring this - yikes. Let's stick with the existing implementation.

[mikofski]

Sorry for being late to the game. Would Numpy trapz work? It's vectorised, and I believe it's the preferred way to integrate.

[mikofski]

The other, probably slower option would be SciPy integrate.quad but I've never used it, and I guess I wouldn't water my time, just use trapz I promise you won't regret it

[kandersolar]

Sorry for being late to the game. Would Numpy trapz work? It's vectorised, and I believe it's the preferred way to integrate.

Is the advantage that we would use trapezoidal integration instead of rectangular integration? I think it would only really be able to replace the two summation lines, so it wouldn't really simplify the code. OK with me if the change is for accuracy, not simplicity. It'll be a little weird to handle the dependency of the solid angle element on zenith; I think the best way is to just lump it into the values to be integrated like it is now.

[mikofski]

Your call, I haven't had a chance to look at the code yet. But in the 2d infinite sheds PR , which is just a bunch more of Bill Marion's formulas I use trapz, I've just always used it, going back to the MATLAB days. I do think trapezoid rule is preferable to rectangle, not only is it more accurate but it avoids the whole left side or right side problem

[kandersolar]

Trying out np.trapz on the test cases, the difference from the current implementation is on the order of 10^-4 at worst, sometimes 10^-5 or -6. Although it's not a perfect comparison because I think the integration bounds would need updating -- with an array of length N, np.sum will sum N intervals while np.trapz integrates only (N-1). Anyway with so small a difference for the test cases, I think I'd prefer np.sum to keep to the reference because the summation is a little nuanced. But I won't stand in the way if we think summing is for chumps and the bikeshed should be shaped like a trapezoid :)

[mikofski]

I'll try to finish reviewing this today, but if you all are already happy, don't let me hold you up.

[mikofski]

LGTM! I left some comments, read them or not, it's all good!

[mikofski]

@kanderso-nrel is this ready to go? There are some "conda-linux" tests that are still not passing. Have you looked into them? Do you know if Will and Cliff's reviews are complete? Shall I merge this if everyone is okay with it?

[kandersolar]

The test failures are unrelated to this PR:

• iotools/test_psm3.py::test_get_psm3_tmy
  • Seems to be an issue on the API side, see NSRDB API returning bad data for tmy-2017 NREL/developer.nrel.gov#159
• iotools/test_tmy.py::test_read_tmy3_remote
  • Hadn't noticed this one before. Looks like the RREDC website has changed and I'm not sure where to find those files anymore. Possibly relevant: https://nsrdb.nrel.gov/data-sets/archives.html

I hesitate to speak for Will and Cliff, but I'd guess the changes since their reviews would be fairly uncontroversial. OK with me to merge.

[wholmgren]

Thanks @kanderso-nrel. Looks good to me. Nice job on the tests.

[mikofski]

Okay @wholmgren, @cwhanse, @kanderso-nrel - I will merge this tomorrow morning unless I hear otherwise. Is squash best?

[wholmgren]

Yes always squash

[mikofski]

Thanks @kanderso-nrel and everyone for reviewing. This is merged now in master.