Add diffuse IAM integration and gallery example

docs/examples/plot_diffuse_aoi_correction.py

@@ -0,0 +1,133 @@
+"""
+Diffuse IAM Calculation
+=======================
+
+Integrating an IAM model across angles to determine the overall reflection
+loss for diffuse irradiance.
+"""
+
+# %%
+# The fraction of light reflected from the front of a module depends on the
+# angle of incidence (AOI) of the light compared to the panel surface.  The
+# steeper the tilt, the larger the reflected fraction is.  The fraction of
+# transmitted light to incident light is called the incident angle modifier
+# (IAM).  Several models exist to calculate the IAM for a given incidence
+# angle (e.g. :py:func:`pvlib.iam.ashrae`, :py:func:`pvlib.iam.martin_ruiz`,
+# :py:func:`pvlib.iam.sapm`, :py:func:`pvlib.iam.physical`).
+# However, evaluating the IAM for diffuse light is
+# not as straightforward because it comes from all directions and therefore
+# has a range of angles of incidence.  Here we show how to integrate the effect
+# of AOI reflection across this AOI range using the process described in [1]_.
+# In particular, we will recreate Figures 3, 4, and 5 in that paper.
+#
+# References
+# ----------
+#  .. [1] B. Marion "Numerical method for angle-of-incidence correction
+#     factors for diffuse radiation incident photovoltaic modules",
+#     Solar Energy, Volume 147, Pages 344-348. 2017.
+#     DOI: 10.1016/j.solener.2017.03.027
+#
+#  .. [2] Duffie, John A. & Beckman, William A. (2013). Solar Engineering
+#     of Thermal Processes.  DOI: 10.1002/9781118671603
+
+
+from pvlib.iam import marion_integrate, physical
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+# %%
+# IAM Model
+# ---------
+#
+# The IAM model used to generate the figures in [1]_ focuses on the air-glass
+# interface.  It uses Snell's, Fresnel's, and Beer's laws to determine the
+# amount of light transmitted through the interface as a function of AOI.
+# The function :py:func:`pvlib.iam.physical` implements this model, except it
+# also includes an exponential term to model attenuation in the glazing layer.
+# To be faithful to Marion's implementation, we will disable this extinction
+# term by setting the attenuation coefficient ``K`` parameter to zero.
+# For more details on this IAM model, see [2]_.
+#
+# Marion generated correction factor profiles for two cases:  a standard
+# uncoated glass with n=1.526 and a glass with anti-reflective (AR) coating
+# with n=1.3.  For convenience, we define a helper function for each case.
+# Comparing them across AOI recreates Figure 3 in [1]_:
+
+
+def calc_uncoated(aoi):
+    return physical(aoi, n=1.526, K=0)
+
+
+def calc_ar_coated(aoi):
+    return physical(aoi, n=1.3, K=0)
+
+
+aoi = np.arange(0, 91)
+cf_uncoated = calc_uncoated(aoi)
+cf_ar_coated = calc_ar_coated(aoi)
+
+plt.plot(aoi, cf_ar_coated, c='b', label='$F_b$, AR coated, n=1.3')
+plt.plot(aoi, cf_uncoated, c='r', label='$F_b$, uncoated, n=1.526')
+plt.xlabel(r'Angle-of-Incidence, AOI $(\degree)$')
+plt.ylabel('Correction Factor')
+plt.legend()
+plt.ylim([0, 1.2])
+plt.grid()
+
+# %%
+# Diffuse sky irradiance (:math:`F_{sky}`)
+# -----------------------------------------
+#
+# Now that we have an AOI model, we use :py:func:`pvlib.iam.marion_integrate`
+# to integrate it across solid angle and determine diffuse irradiance
+# correction factors.  Marion defines three types of diffuse irradiance:
+# sky, horizon, and ground-reflected.  The IAM correction factor is evaluated
+# independently for each type.
+# First we recreate Figure 4 in [1]_, showing the dependence of the sky diffuse
+# correction factor on module tilt.
+
+tilts = np.arange(0, 91, 2.5)
+
+iam_uncoated = marion_integrate(calc_uncoated, tilts, 'sky', N=180)
+iam_ar_coated = marion_integrate(calc_ar_coated, tilts, 'sky', N=180)
+
+plt.plot(tilts, iam_ar_coated, c='b', marker='^',
+         label='$F_{sky}$, AR coated, n=1.3')
+plt.plot(tilts, iam_uncoated, c='r', marker='x',
+         label='$F_{sky}$, uncoated, n=1.526')
+plt.ylim([0.9, 1.0])
+plt.xlabel(r'PV Module Tilt, $\beta (\degree)$')
+plt.ylabel('Correction Factor')
+plt.grid()
+plt.legend()
+plt.show()
+
+
+# %%
+# Diffuse horizon and ground irradiance (:math:`F_{hor}, F_{grd}`)
+# -----------------------------------------------------------------
+#
+# Now we recreate Figure 5 in [1]_, showing the dependence of the correction
+# factors for horizon and ground diffuse irradiance on module tilt.  Note that
+# we use 1800 points instead of 180 for the horizon case to match [1]_.
+
+iam_uncoated_grd = marion_integrate(calc_uncoated, tilts, 'ground', N=180)
+iam_ar_coated_grd = marion_integrate(calc_ar_coated, tilts, 'ground', N=180)
+
+iam_uncoated_hor = marion_integrate(calc_uncoated, tilts, 'horizon', N=1800)
+iam_ar_coated_hor = marion_integrate(calc_ar_coated, tilts, 'horizon', N=1800)
+
+plt.plot(tilts, iam_ar_coated_hor, c='b', marker='^',
+         label='$F_{hor}$, AR coated, n=1.3')
+plt.plot(tilts, iam_uncoated_hor, c='r', marker='x',
+         label='$F_{hor}$, uncoated, n=1.526')
+plt.plot(tilts, iam_ar_coated_grd, c='b', marker='s',
+         label='$F_{grd}$, AR coated, n=1.3')
+plt.plot(tilts, iam_uncoated_grd, c='r', marker='+',
+         label='$F_{grd}$, uncoated, n=1.526')
+plt.xlabel(r'PV Module Tilt, $\beta (\degree)$')
+plt.ylabel('Correction Factor')
+plt.grid()
+plt.legend()
+plt.show()

docs/sphinx/source/api.rst

@@ -223,6 +223,7 @@ Incident angle modifiers
    iam.martin_ruiz_diffuse
    iam.sapm
    iam.interp
+   iam.marion_integrate
 
 PV temperature models
 ---------------------

docs/sphinx/source/whatsnew/v0.8.0.rst

@@ -8,6 +8,8 @@ API Changes
 
 Enhancements
 ~~~~~~~~~~~~
+* Add :py:func:`pvlib.iam.marion_integrate` to calculate IAM values for
+  diffuse irradiance. (:pull:`984`)
 
 Bug fixes
 ~~~~~~~~~
@@ -26,6 +28,8 @@ Documentation
 * Clarify units for heat loss factors in
   :py:func:`pvlib.temperature.pvsyst_cell` and
   :py:func:`pvlib.temperature.faiman`. (:pull:`960`)
+* Add a gallery example of calculating diffuse IAM using
+  :py:func:`pvlib.iam.marion_integrate`. (:pull:`984`)
 
 Requirements
 ~~~~~~~~~~~~

pvlib/iam.py

@@ -527,3 +527,120 @@ def sapm(aoi, module, upper=None):
         iam = pd.Series(iam, aoi.index)
 
     return iam
+
+
+def marion_integrate(function, surface_tilt, region, N=180):
+    """
+    Integrate an incidence angle modifier (IAM) function over solid angle
+    to determine a diffuse irradiance correction factor using Marion's method.
+
+    Parameters
+    ----------
+    function : callable func(aoi)
+        The IAM function to evaluate across solid angle.  The function must
+        be vectorized and take only one parameter, the angle of incidence in
+        degrees.
+
+    surface_tilt : numeric
+        Surface tilt angles in decimal degrees.
+        The tilt angle is defined as degrees from horizontal
+        (e.g. surface facing up = 0, surface facing horizon = 90).
+
+    region : {'sky', 'horizon', 'ground'}
+        The region to integrate over.  Must be one of:
+
+            * 'sky': radiation from the sky dome
+            * 'horizon': radiation from the region of the sky near the horizon
+            * 'ground': radiation reflected from the ground
+
+        See [1]_ for a detailed description of each class.
+
+    N : int, default 180
+        The number of increments in the zenith integration.
+        If ``region=='horizon'``, a larger number should be used
+        (1800 was used in [1]_).
+
+    Returns
+    -------
+    Fd : numeric
+        AOI diffuse correction factor for the specified region.
+
+    References
+    ----------
+    .. [1] B. Marion "Numerical method for angle-of-incidence correction
+       factors for diffuse radiation incident photovoltaic modules",
+       Solar Energy, Volume 147, Pages 344-348. 2017.
+       DOI: 10.1016/j.solener.2017.03.027
+
+    Examples
+    --------
+    >>> marion_integrate(pvlib.iam.ashrae, 20, 'sky')
+    0.9596085829811408
+
+    >>> from functools import partial
+    >>> f = partial(pvlib.iam.physical, n=1.3)
+    >>> marion_integrate(f, [20, 30], 'sky')
+    array([0.96225034, 0.9653219 ])
+    """
+    beta = np.radians(surface_tilt)
+    if isinstance(beta, pd.Series):
+        # convert Series to np array for broadcasting later
+        beta = beta.values
+    ai = np.pi/N  # angular increment
+
+    phi_range = np.linspace(0, np.pi, N, endpoint=False)
+    psi_range = np.linspace(0, 2*np.pi, 2*N, endpoint=False)
+
+    # the pseudocode in [1] do these checks at the end, but it's
+    # faster to do this criteria check up front instead of later.
+    if region == 'sky':
+        mask = phi_range + ai <= np.pi/2
+    elif region == 'horizon':
+        lo = 89.5 * np.pi/180
+        hi = np.pi/2
+        mask = (lo <= phi_range) & (phi_range + ai <= hi)
+    elif region == 'ground':
+        mask = (phi_range >= np.pi/2)
+    else:
+        raise ValueError('Invalid region: {}'.format(region))
+    phi_range = phi_range[mask]
+
+    # fast Cartesian product of phi and psi
+    angles = np.array(np.meshgrid(phi_range, psi_range)).T.reshape(-1, 2)
+    # index with single-element lists to maintain 2nd dimension so that
+    # these angle arrays broadcast across the beta array
+    phi_1 = angles[:, [0]]
+    psi_1 = angles[:, [1]]
+    phi_2 = phi_1 + ai
+    # psi_2 = psi_1 + ai  # not needed
+    phi_avg = phi_1 + 0.5*ai
+    psi_avg = psi_1 + 0.5*ai
+    term_1 = np.cos(beta) * np.cos(phi_avg)
+    # we are free to choose any gamma between zero and 2pi, so choose zero
+    # and omit it from the cos(psi_avg) term
+    term_2 = np.sin(beta) * np.sin(phi_avg) * np.cos(psi_avg)
+    aoi = np.arccos(term_1 + term_2)
+    # simplify Eq 8, (psi_2 - psi_1) is always ai
+    dAs = ai * (np.cos(phi_1) - np.cos(phi_2))
+    cosaoi_dAs = np.cos(aoi) * dAs
+    # apply the final AOI check, zeroing out non-passing points
+    mask = aoi < np.pi/2
+    cosaoi_dAs = np.where(mask, cosaoi_dAs, 0)
+    numerator = np.sum(function(np.degrees(aoi)) * cosaoi_dAs, axis=0)
+    denominator = np.sum(cosaoi_dAs, axis=0)
+
+    with np.errstate(invalid='ignore'):
+        # in some cases, no points pass the criteria
+        # (e.g. region='ground', surface_tilt=0), so we override the division
+        # by zero to set Fd=0.  Also, preserve nans in beta.
+        Fd = np.where((denominator != 0) | ~np.isfinite(beta),
+                      numerator / denominator,
+                      0)
+
+    # preserve input type
+    if np.isscalar(surface_tilt):
+        Fd = Fd.item()
+    elif isinstance(surface_tilt, pd.Series):
+        Fd = pd.Series(Fd, surface_tilt.index)
+
+    return Fd

pvlib/tests/test_iam.py

@@ -218,3 +218,61 @@ def test_sapm_limits():
 
     module_parameters = {'B0': -5, 'B1': 0, 'B2': 0, 'B3': 0, 'B4': 0, 'B5': 0}
     assert _iam.sapm(1, module_parameters) == 0
+
+
+@pytest.mark.parametrize('region,N,expected', [
+    ('sky', 180, 0.9596085829811408),
+    ('horizon', 1800, 0.8329070417832541),
+    ('ground', 180, 0.719823559106309)
+])
+def test_marion_integrate_scalar(region, N, expected):
+    actual = _iam.marion_integrate(_iam.ashrae, 20, region, N)
+    assert_allclose(actual, expected)
+
+    with np.errstate(invalid='ignore'):
+        actual = _iam.marion_integrate(_iam.ashrae, np.nan, region, N)
+    expected = np.nan
+    assert_allclose(actual, expected)
+
+
+@pytest.mark.parametrize('region,N,expected', [
+    ('sky', 180, [0.9523611991069362, 0.9596085829811408, 0.9619811198105501]),
+    ('horizon', 1800, [0.0, 0.8329070417832541, 0.8987287652347437]),
+    ('ground', 180, [0.0, 0.719823559106309, 0.8186360238536674])
+])
+def test_marion_integrate_list(region, N, expected):
+    actual = _iam.marion_integrate(_iam.ashrae, [0, 20, 30], region, N)
+    assert_allclose(actual, expected)
+
+    with np.errstate(invalid='ignore'):
+        actual = _iam.marion_integrate(_iam.ashrae, [0, 20, np.nan], region, N)
+    assert_allclose(actual, [expected[0], expected[1], np.nan])
+
+
+@pytest.mark.parametrize('region,N,expected', [
+    ('sky', 180, [0.9523611991069362, 0.9596085829811408, 0.9619811198105501]),
+    ('horizon', 1800, [0.0, 0.8329070417832541, 0.8987287652347437]),
+    ('ground', 180, [0.0, 0.719823559106309, 0.8186360238536674])
+])
+def test_marion_integrate_series(region, N, expected):
+    idx = pd.date_range('2019-01-01', periods=3, freq='h')
+    tilt = pd.Series([0, 20, 30], index=idx)
+    expected = pd.Series(expected, index=idx)
+    actual = _iam.marion_integrate(_iam.ashrae, tilt, region, N)
+    assert_series_equal(actual, expected)
+
+    tilt.iloc[1] = np.nan
+    expected.iloc[1] = np.nan
+    with np.errstate(invalid='ignore'):
+        actual = _iam.marion_integrate(_iam.ashrae, tilt, region, N)
+    assert_allclose(actual, expected)
+
+
+def test_marion_integrate_ground_flat():
+    iam = _iam.marion_integrate(_iam.ashrae, 0, 'horizon', 1800)
+    assert_allclose(iam, 0)
+
+
+def test_marion_integrate_invalid():
+    with pytest.raises(ValueError):
+        _iam.marion_integrate(_iam.ashrae, 0, 'bad', 180)