add Boland DF estimation

docs/examples/plot_diffuse_fraction.py

@@ -0,0 +1,213 @@
+"""
+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 numpy as np
+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.
+solpos = get_solarposition(
+    greensboro.index, latitude=metadata['latitude'],
+    longitude=metadata['longitude'], altitude=metadata['altitude'],
+    pressure=greensboro.Pressure*100,  # convert from millibar to Pa
+    temperature=greensboro.DryBulb)
+
+# %%
+# 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)
+out_disc = out_disc.rename(columns={'dni': 'dni_disc'})
+out_disc['dhi_disc'] = (
+    greensboro.GHI
+    - out_disc.dni_disc*np.cos(np.radians(solpos.apparent_zenith)))
+
+# %%
+# 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)
+dhi_dirint = (
+    greensboro.GHI
+    - dni_dirint*np.cos(np.radians(solpos.apparent_zenith)))
+out_dirint = pd.DataFrame(
+    {'dni_dirint': dni_dirint, 'dhi_dirint': dhi_dirint},
+    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
+
+JAN6AM, JAN6PM = '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[JAN6AM:JAN6PM].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[JAN6AM:JAN6PM].plot(ax=ax[1])
+ax[1].grid(which="both")
+ax[1].set_ylabel('DHI $[W/m^2]$')
+ghi_kt[JAN6AM:JAN6PM].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 ...
+
+APR6AM, APR6PM = '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[APR6AM:APR6PM].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[APR6AM:APR6PM].plot(ax=ax[1])
+ax[1].grid(which="both")
+ax[1].set_ylabel('DHI $[W/m^2]$')
+ghi_kt[APR6AM:APR6PM].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.
+
+JUL6AM, JUL6PM = '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[JUL6AM:JUL6PM].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[JUL6AM:JUL6PM].plot(ax=ax[1])
+ax[1].grid(which="both")
+ax[1].set_ylabel('DHI $[W/m^2]$')
+ghi_kt[JUL6AM:JUL6PM].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 with only kt, which is derived
+# from 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`)
 
 Bug fixes
 ~~~~~~~~~

pvlib/irradiance.py

@@ -2265,6 +2265,94 @@ def erbs(ghi, zenith, datetime_or_doy, min_cos_zenith=0.065, max_zenith=87):
     return data
 
 
+def boland(ghi, zenith, datetime_or_doy, 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 extraterrestrial
+    irradiance.
+
+    .. math::
+
+        \mathit{DF} = \frac{1}{1 + \exp\left(-5 + 8.6 k_t\right)}
+
+
+    Parameters
+    ----------
+    ghi: numeric
+        Global horizontal irradiance. [W/m^2]
+    zenith: numeric
+        True (not refraction-corrected) zenith angles in decimal degrees.
+    datetime_or_doy : int, float, numpy.ndarray, pandas.DatetimeIndex
+        Day of year or array of days of year e.g.
+        pd.DatetimeIndex.dayofyear, or pd.DatetimeIndex.
+    min_cos_zenith : numeric, default 0.065
+        Minimum value of cos(zenith) to allow when calculating global
+        clearness index `kt`. 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
+    """
+
+    dni_extra = get_extra_radiation(datetime_or_doy)
+
+    kt = clearness_index(ghi, zenith, dni_extra, min_cos_zenith=min_cos_zenith,
+                         max_clearness_index=1)
+
+    # Boland equation
+    df = 1.0 / (1.0 + np.exp(-5.0 + 8.6 * kt))
+    # NOTE: [1] 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(zenith)
+    bad_values = (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],
+         [103.735879, 40.958822,  0.405724],
+         [776.006568, 235.782716, 0.718133],
+         [845.794317, 167.055199, 0.768214]]),
+        columns=['dni', 'dhi', 'kt'], index=index)
+
+    out = irradiance.boland(ghi, zenith, index)
+
+    assert_frame_equal(np.round(out, 0), np.round(expected, 0))
+
+
 def test_erbs_min_cos_zenith_max_zenith():
     # map out behavior under difficult conditions with various
     # limiting kwargs settings

readthedocs.yml

@@ -1,7 +1,26 @@
+# .readthedocs.yaml
+# Read the Docs configuration file
+# See https://docs.readthedocs.io/en/stable/config-file/v2.html for details
+
+# Required
+version: 2
+
+build:
+   os: ubuntu-20.04
+   tools:
+     python: "3.7"
+   jobs:
+     # fetch the full history so that setuptools_scm can determine the
+     # correct version string for long PRs with many commits
+     post_checkout:
+       - git fetch --unshallow
+
 python:
-    version: 3
-    # only use the packages specified in setup.py
-    use_system_site_packages: false
-    pip_install: true
-    extra_requirements:
-        - doc
+   # only use the packages specified in setup.py
+   system_packages: false
+
+   install:
+      - method: pip
+        path: .
+        extra_requirements:
+          - doc