add Boland DF estimation (#1179)

* add Boland DF estimation

* update what's new, add to docs

* add example

* finish example

* trailing whitespace

* respond to comments

- put Badescu textbook before paper
- revise wording defining kt

* respond to comments

- remove redundant kt definition
- use intersphinx to link to pandas and numpy types

* Update docs/examples/plot_diffuse_fraction.py

responding to comments, simplify wording L12-15 to omit "systems often tilted to optimize performance..."

Co-authored-by: Cliff Hansen <[email protected]>

* Update docs/examples/plot_diffuse_fraction.py

respond to comments, reword decomposition example intro

Co-authored-by: Cliff Hansen <[email protected]>

* respond to comments

- add intro before functions in example

* respond to comments in diffuse fraction example

- fix use :py:func: to xref functions in docs, instead of :meth:
- add comment to explain conversion of mbars to Pa
- use Gsc not E0 in kt comparison plot

* respond to comments

- add section to conclusions to warn users that TMY3 and NSRDB are also models,
 therefore not to assume differences are errors,
 and link to TMY3 & NSRDB documentation
- use implicit targets to link to DISC, DIRINT, Erbs, & Boland sections
- reverse change of scaled GHI by E0 to Gsc (aka: solar constant)
- use math directive for DNI formula

* respond to comments

- revise & condense conclusions, don't refer to differences as errors
 without operational data
- add header before plots section, explain why comparing normalized GHI
- add subheaders for seasonal plots

* fix stickler

* update readthedocs.yml to v2

* oops, extra requires is doc, no s

* use requirements to pin requirements

* install pvlib again?

* Revert "update readthedocs.yml to v2"

This reverts commit ccf40a8.

* update RTD config to fix shallow clone issue

* Update pvlib/irradiance.py

Use Sphinx math role to render k_t in docstring for min cosine zenith

Co-authored-by: Cliff Hansen <[email protected]>

* use solar_zenith in irradiance.boland docstring since #1403

* update reference to modeling diffuse fraction by J. Boland, available online https://www.researchgate.net/publication/229873508_Modelling_the_diffuse_fraction_of_global_solar_radiation_on_a_horizontal_surface

* let user pass boland coeffs as kwargs

- add a_coeff, b_coeff to docstring params
- update equation to match Boland paper, use A, B coeffs
- update coeffs to match Boland paper too (8.6-> 8.645 and 5 -> 5.3
- add note to explain different coeffs for different time intervals
- give coeffs for 15-min & 1-hr
- update zenith to solar_zenith everywhere
- update test expected to match output with new coeffs

* move plot diffuse fraction example to subdir

- create irradiance-decomposition subdir in gallery, add readme.rst
- use complete_irradiance() to get DHI using closure eqn's

* shift solar position calc to hour center

- in plot diffuse fraction examples
- add note to explain b/c TMY timestampe at hour end
- reset index to align with TMY3
- fix Jan, July, Apr timestamp names didn't match actual times

* don't need numpy in plot diffuse fraction examples

* minor edits in plot diffuse fraction example

- chnage df to be specific to disc or dirint
- change date names to Jan-04, Jan-07, etc. instead of _AM, _PM etc.

* Apply suggestions from code review by Adam

* use lower case coeffs in equation
* periods in docstrings
* some wordsmithing
* replace parameter types for datetime_or_doy with numeric
* consistent def for kt is ratio of GHI to horizontal extraterrestrial-irradiance

Co-authored-by: Adam R. Jensen <[email protected]>

* fix stickler - long lines in plot SF example

---------

Co-authored-by: Cliff Hansen <[email protected]>
Co-authored-by: Kevin Anderson <[email protected]>
Co-authored-by: Adam R. Jensen <[email protected]>

Author: 4 people

docs/examples/irradiance-decomposition/README.rst

@@ -0,0 +1,3 @@
+Irradiance Decomposition
+------------------------
+

docs/examples/irradiance-decomposition/plot_diffuse_fraction.py

@@ -0,0 +1,219 @@
+"""
+Diffuse Fraction Estimation
+===========================
+
+Comparison of diffuse fraction estimation methods used to derive direct and
+diffuse components from measured global horizontal irradiance.
+"""
+
+# %%
+# This example demonstrates how to use diffuse fraction estimation methods to
+# obtain direct and diffuse components from measured global horizontal
+# irradiance (GHI). Irradiance sensors such as pyranometers typically only
+# measure GHI. pvlib provides several functions that can be used to separate
+# GHI into the diffuse and direct components. The separate components are
+# needed to estimate the total irradiance on a tilted surface.
+
+import pathlib
+from matplotlib import pyplot as plt
+import pandas as pd
+from pvlib.iotools import read_tmy3
+from pvlib.solarposition import get_solarposition
+from pvlib import irradiance
+import pvlib
+
+# For this example we use the Greensboro, North Carolina, TMY3 file which is
+# in the pvlib data directory. TMY3 are made from the median months from years
+# of data measured from 1990 to 2010. Therefore we change the timestamps to a
+# common year, 1990.
+DATA_DIR = pathlib.Path(pvlib.__file__).parent / 'data'
+greensboro, metadata = read_tmy3(DATA_DIR / '723170TYA.CSV', coerce_year=1990)
+
+# Many of the diffuse fraction estimation methods require the "true" zenith, so
+# we calculate the solar positions for the 1990 at Greensboro, NC.
+# NOTE: TMY3 files timestamps indicate the end of the hour, so shift indices
+# back 30-minutes to calculate solar position at center of the interval
+solpos = get_solarposition(
+    greensboro.index.shift(freq="-30T"), latitude=metadata['latitude'],
+    longitude=metadata['longitude'], altitude=metadata['altitude'],
+    pressure=greensboro.Pressure*100,  # convert from millibar to Pa
+    temperature=greensboro.DryBulb)
+solpos.index = greensboro.index  # reset index to end of the hour
+
+# %%
+# pvlib Decomposition Functions
+# -----------------------------
+# Methods for separating DHI into diffuse and direct components include:
+# `DISC`_, `DIRINT`_, `Erbs`_, and `Boland`_.
+
+# %%
+# DISC
+# ----
+#
+# DISC :py:func:`~pvlib.irradiance.disc` is an empirical correlation developed
+# at SERI (now NREL) in 1987. The direct normal irradiance (DNI) is related to
+# clearness index (kt) by two polynomials split at kt = 0.6, then combined with
+# an exponential relation with airmass.
+
+out_disc = irradiance.disc(
+    greensboro.GHI, solpos.zenith, greensboro.index, greensboro.Pressure*100)
+# use "complete sum" AKA "closure" equations: DHI = GHI - DNI * cos(zenith)
+df_disc = irradiance.complete_irradiance(
+    solar_zenith=solpos.apparent_zenith, ghi=greensboro.GHI, dni=out_disc.dni,
+    dhi=None)
+out_disc = out_disc.rename(columns={'dni': 'dni_disc'})
+out_disc['dhi_disc'] = df_disc.dhi
+
+# %%
+# DIRINT
+# ------
+#
+# DIRINT :py:func:`~pvlib.irradiance.dirint` is a modification of DISC
+# developed by Richard Perez and Pierre Ineichen in 1992.
+
+dni_dirint = irradiance.dirint(
+    greensboro.GHI, solpos.zenith, greensboro.index, greensboro.Pressure*100,
+    temp_dew=greensboro.DewPoint)
+# use "complete sum" AKA "closure" equation: DHI = GHI - DNI * cos(zenith)
+df_dirint = irradiance.complete_irradiance(
+    solar_zenith=solpos.apparent_zenith, ghi=greensboro.GHI, dni=dni_dirint,
+    dhi=None)
+out_dirint = pd.DataFrame(
+    {'dni_dirint': dni_dirint, 'dhi_dirint': df_dirint.dhi},
+    index=greensboro.index)
+
+# %%
+# Erbs
+# ----
+#
+# The Erbs method, :py:func:`~pvlib.irradiance.erbs` developed by Daryl Gregory
+# Erbs at the University of Wisconsin in 1982 is a piecewise correlation that
+# splits kt into 3 regions: linear for kt <= 0.22, a 4th order polynomial
+# between 0.22 < kt <= 0.8, and a horizontal line for kt > 0.8.
+
+out_erbs = irradiance.erbs(greensboro.GHI, solpos.zenith, greensboro.index)
+out_erbs = out_erbs.rename(columns={'dni': 'dni_erbs', 'dhi': 'dhi_erbs'})
+
+# %%
+# Boland
+# ------
+#
+# The Boland method, :py:func:`~pvlib.irradiance.boland` is a single logistic
+# exponential correlation that is continuously differentiable and bounded
+# between zero and one.
+
+out_boland = irradiance.boland(greensboro.GHI, solpos.zenith, greensboro.index)
+out_boland = out_boland.rename(
+    columns={'dni': 'dni_boland', 'dhi': 'dhi_boland'})
+
+# %%
+# Comparison Plots
+# ----------------
+# In the plots below we compare the four decomposition models to the TMY3 file
+# for Greensboro, North Carolina. We also compare the clearness index, kt, with
+# GHI normalized by a reference irradiance, E0 = 1000 [W/m^2], to highlight
+# spikes caused when cosine of zenith approaches zero, particularly at sunset.
+#
+# First we combine the dataframes for the decomposition models and the TMY3
+# file together to make plotting easier.
+
+dni_renames = {
+    'DNI': 'TMY3', 'dni_disc': 'DISC', 'dni_dirint': 'DIRINT',
+    'dni_erbs': 'Erbs', 'dni_boland': 'Boland'}
+dni = [
+    greensboro.DNI, out_disc.dni_disc, out_dirint.dni_dirint,
+    out_erbs.dni_erbs, out_boland.dni_boland]
+dni = pd.concat(dni, axis=1).rename(columns=dni_renames)
+dhi_renames = {
+    'DHI': 'TMY3', 'dhi_disc': 'DISC', 'dhi_dirint': 'DIRINT',
+    'dhi_erbs': 'Erbs', 'dhi_boland': 'Boland'}
+dhi = [
+    greensboro.DHI, out_disc.dhi_disc, out_dirint.dhi_dirint,
+    out_erbs.dhi_erbs, out_boland.dhi_boland]
+dhi = pd.concat(dhi, axis=1).rename(columns=dhi_renames)
+ghi_kt = pd.concat([greensboro.GHI/1000.0, out_erbs.kt], axis=1)
+
+# %%
+# Winter
+# ++++++
+# Finally, let's plot them for a few winter days and compare
+
+JAN04, JAN07 = '1990-01-04 00:00:00-05:00', '1990-01-07 23:59:59-05:00'
+f, ax = plt.subplots(3, 1, figsize=(8, 10), sharex=True)
+dni[JAN04:JAN07].plot(ax=ax[0])
+ax[0].grid(which="both")
+ax[0].set_ylabel('DNI $[W/m^2]$')
+ax[0].set_title('Comparison of Diffuse Fraction Estimation Methods')
+dhi[JAN04:JAN07].plot(ax=ax[1])
+ax[1].grid(which="both")
+ax[1].set_ylabel('DHI $[W/m^2]$')
+ghi_kt[JAN04:JAN07].plot(ax=ax[2])
+ax[2].grid(which='both')
+ax[2].set_ylabel(r'$\frac{GHI}{E0}, k_t$')
+f.tight_layout()
+
+# %%
+# Spring
+# ++++++
+# And a few spring days ...
+
+APR04, APR07 = '1990-04-04 00:00:00-05:00', '1990-04-07 23:59:59-05:00'
+f, ax = plt.subplots(3, 1, figsize=(8, 10), sharex=True)
+dni[APR04:APR07].plot(ax=ax[0])
+ax[0].grid(which="both")
+ax[0].set_ylabel('DNI $[W/m^2]$')
+ax[0].set_title('Comparison of Diffuse Fraction Estimation Methods')
+dhi[APR04:APR07].plot(ax=ax[1])
+ax[1].grid(which="both")
+ax[1].set_ylabel('DHI $[W/m^2]$')
+ghi_kt[APR04:APR07].plot(ax=ax[2])
+ax[2].grid(which='both')
+ax[2].set_ylabel(r'$\frac{GHI}{E0}, k_t$')
+f.tight_layout()
+
+# %%
+# Summer
+# ++++++
+# And few summer days to finish off the seasons.
+
+JUL04, JUL07 = '1990-07-04 00:00:00-05:00', '1990-07-07 23:59:59-05:00'
+f, ax = plt.subplots(3, 1, figsize=(8, 10), sharex=True)
+dni[JUL04:JUL07].plot(ax=ax[0])
+ax[0].grid(which="both")
+ax[0].set_ylabel('DNI $[W/m^2]$')
+ax[0].set_title('Comparison of Diffuse Fraction Estimation Methods')
+dhi[JUL04:JUL07].plot(ax=ax[1])
+ax[1].grid(which="both")
+ax[1].set_ylabel('DHI $[W/m^2]$')
+ghi_kt[JUL04:JUL07].plot(ax=ax[2])
+ax[2].grid(which='both')
+ax[2].set_ylabel(r'$\frac{GHI}{E0}, k_t$')
+f.tight_layout()
+
+# %%
+# Conclusion
+# ----------
+# This example compares several decomposition models to a TMY3 file for
+# Greensboro, North Carolina. However, DNI and DHI in TMY3 files are themselves
+# the output of models (either METSTAT or SUNY), and so differences between
+# *e.g.* DISC output and the TMY3 file shouldn't be regarded as errors, and
+# it's not a reasonable expectation to assume that the four models should
+# reproduce the TMY3 values. We refer those interested to the `TMY3`_ and
+# `NSRDB`_ user manuals.
+#
+# The Erbs and Boland models are correlations only based on the clearness index
+# kt, which is the ratio of GHI to the the horizontal component of the
+# extra-terrestrial irradiance. At low sun elevation (zenith near 90 degrees),
+# especially near sunset, kt can explode because the denominator
+# (extra-terrestrial irradiance) approaches zero. In pvlib this behavior is
+# moderated by ``min_cos_zenith`` and ``max_clearness_index`` which each have
+# reasonable defaults. Even so, near sunset there are still spikes in kt and
+# DNI from Erbs and Boland for Jan. 5th & 7th, April 4th, 5th, & 7th, and July
+# 6th & 7th.
+#
+# By contrast, the DISC and DIRINT methods estimate DNI first by means of
+# correlations, which include additional variables such as airmass. These
+# methods seem to reduce DNI spikes over 1000 [W/m^2].
+#
+# .. _TMY3: https://www.nrel.gov/docs/fy08osti/43156.pdf
+# .. _NSRDB: https://www.nrel.gov/docs/fy12osti/54824.pdf

docs/sphinx/source/reference/irradiance/decomposition.rst

@@ -12,5 +12,6 @@ DNI estimation models
    irradiance.dirint
    irradiance.dirindex
    irradiance.erbs
+   irradiance.boland
    irradiance.campbell_norman
    irradiance.gti_dirint

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

@@ -20,6 +20,10 @@ Enhancements
   :py:func:`pvlib.snow.loss_townsend` (:issue:`1636`, :pull:`1653`)
 * Added optional ``n_ar`` parameter to :py:func:`pvlib.iam.physical` to
   support an anti-reflective coating. (:issue:`1501`, :pull:`1616`)
+* :py:func:`~pvlib.irradiance.boland` is another diffuse fraction, DF,
+  estimation method similar to Erbs but uses a single logistic exponential
+  correlation between DF and clearness index, kt, that is continuously
+  differentiable and bounded between zero and one. (:pull:`1179`)
 * Add ``model='gueymard2003'``, the airmass model used for REST and REST2,
   to :py:func:`~pvlib.atmosphere.get_relative_airmass`. (:pull:`1655`)
 

pvlib/irradiance.py

@@ -2265,6 +2265,110 @@ def erbs(ghi, zenith, datetime_or_doy, min_cos_zenith=0.065, max_zenith=87):
     return data
 
 
+def boland(ghi, solar_zenith, datetime_or_doy, a_coeff=8.645, b_coeff=0.613,
+           min_cos_zenith=0.065, max_zenith=87):
+    r"""
+    Estimate DNI and DHI from GHI using the Boland clearness index model.
+
+    The Boland model [1]_, [2]_ estimates the diffuse fraction, DF, from global
+    horizontal irradiance, GHI, through an empirical relationship between DF
+    and the clearness index, :math:`k_t`, the ratio of GHI to horizontal
+    extraterrestrial irradiance.
+
+    .. math::
+
+        \mathit{DF} = \frac{1}{1 + \exp\left(a \left(k_t - b\right)\right)}
+
+
+    Parameters
+    ----------
+    ghi: numeric
+        Global horizontal irradiance. [W/m^2]
+    solar_zenith: numeric
+        True (not refraction-corrected) zenith angles in decimal degrees.
+    datetime_or_doy : numeric, pandas.DatetimeIndex
+        Day of year or array of days of year e.g.
+        pd.DatetimeIndex.dayofyear, or pd.DatetimeIndex.
+    a_coeff : float, default 8.645
+        Logistic curve fit coefficient.
+    b_coeff : float, default 0.613
+        Logistic curve fit coefficient.
+    min_cos_zenith : numeric, default 0.065
+        Minimum value of cos(zenith) to allow when calculating global
+        clearness index :math:`k_t`. Equivalent to zenith = 86.273 degrees.
+    max_zenith : numeric, default 87
+        Maximum value of zenith to allow in DNI calculation. DNI will be
+        set to 0 for times with zenith values greater than `max_zenith`.
+
+    Returns
+    -------
+    data : OrderedDict or DataFrame
+        Contains the following keys/columns:
+
+            * ``dni``: the modeled direct normal irradiance in W/m^2.
+            * ``dhi``: the modeled diffuse horizontal irradiance in
+              W/m^2.
+            * ``kt``: Ratio of global to extraterrestrial irradiance
+              on a horizontal plane.
+
+    References
+    ----------
+    .. [1] J. Boland, B. Ridley (2008) Models of Diffuse Solar Fraction. In:
+       Badescu V. (eds) Modeling Solar Radiation at the Earth’s Surface.
+       Springer, Berlin, Heidelberg. :doi:`10.1007/978-3-540-77455-6_8`
+    .. [2] John Boland, Lynne Scott, and Mark Luther, Modelling the diffuse
+       fraction of global solar radiation on a horizontal surface,
+       Environmetrics 12(2), pp 103-116, 2001,
+       :doi:`10.1002/1099-095X(200103)12:2%3C103::AID-ENV447%3E3.0.CO;2-2`
+
+    See also
+    --------
+    dirint
+    disc
+    erbs
+
+    Notes
+    -----
+    Boland diffuse fraction differs from other decomposition algorithms by use
+    of a logistic function to fit the entire range of clearness index,
+    :math:`k_t`. Parameters ``a_coeff`` and ``b_coeff`` are reported in [2]_
+    for different time intervals:
+
+    * 15-minute: ``a = 8.645`` and ``b = 0.613``
+    * 1-hour:  ``a = 7.997`` and ``b = 0.586``
+    """
+
+    dni_extra = get_extra_radiation(datetime_or_doy)
+
+    kt = clearness_index(
+        ghi, solar_zenith, dni_extra, min_cos_zenith=min_cos_zenith,
+        max_clearness_index=1)
+
+    # Boland equation
+    df = 1.0 / (1.0 + np.exp(a_coeff * (kt - b_coeff)))
+    # NOTE: [2] has different coefficients, for different time intervals
+    # 15-min: df = 1 / (1 + exp(8.645 * (kt - 0.613)))
+    # 1-hour: df = 1 / (1 + exp(7.997 * (kt - 0.586)))
+
+    dhi = df * ghi
+
+    dni = (ghi - dhi) / tools.cosd(solar_zenith)
+    bad_values = (solar_zenith > max_zenith) | (ghi < 0) | (dni < 0)
+    dni = np.where(bad_values, 0, dni)
+    # ensure that closure relationship remains valid
+    dhi = np.where(bad_values, ghi, dhi)
+
+    data = OrderedDict()
+    data['dni'] = dni
+    data['dhi'] = dhi
+    data['kt'] = kt
+
+    if isinstance(datetime_or_doy, pd.DatetimeIndex):
+        data = pd.DataFrame(data, index=datetime_or_doy)
+
+    return data
+
+
 def campbell_norman(zenith, transmittance, pressure=101325.0,
                     dni_extra=1367.0):
     '''

pvlib/tests/test_irradiance.py

@@ -804,6 +804,22 @@ def test_erbs():
     assert_frame_equal(np.round(out, 0), np.round(expected, 0))
 
 
+def test_boland():
+    index = pd.DatetimeIndex(['20190101']*3 + ['20190620'])
+    ghi = pd.Series([0, 50, 1000, 1000], index=index)
+    zenith = pd.Series([120, 85, 10, 10], index=index)
+    expected = pd.DataFrame(np.array(
+        [[0.0,        0.0,        0.0],
+         [81.9448546, 42.8580353, 0.405723511],
+         [723.764990, 287.230626, 0.718132729],
+         [805.020419, 207.209650, 0.768214312]]),
+        columns=['dni', 'dhi', 'kt'], index=index)
+
+    out = irradiance.boland(ghi, zenith, index)
+
+    assert np.allclose(out, expected)
+
+
 def test_erbs_min_cos_zenith_max_zenith():
     # map out behavior under difficult conditions with various
     # limiting kwargs settings