add Boland DF estimation

docs/examples/plot_diffuse_fraction.py

@@ -0,0 +1,186 @@
+"""
+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). PV systems are often tilted to optimize performance, so the
+# entire diffuse sky dome may not be visible to the PV surface. Determining the
+# total irradiance incident on the plane of the array requires transposing the
+# diffuse component, but irradiance sensors such as pyranometers typically only
+# measure GHI. Therefore pvlib provides several correlations to estimate the
+# diffuse fraction of the GHI, that can be used to resolve the diffuse and
+# direct components.
+
+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, temperature=greensboro.DryBulb)
+
+# %%
+# DISC
+# ----
+#
+# DISC :py:meth:`~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:meth:`~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:meth:`~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:meth:`~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'})
+
+# %%
+# Combine everything together.
+
+dni_renames = {
+    'DNI': 'TMY', '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': 'TMY', '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)
+
+# %%
+# 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()
+
+# %%
+# 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()
+
+# %%
+# 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
+# ----------
+# The Erbs and Boland are correlations with only kt, which is derived from the
+# horizontal component of the extra-terrestrial irradiance. Therefore at low
+# sun elevation (zenith ~ 90-deg), especially near sunset, this causes kt to
+# explode as the denominator approaches zero. This is controlled in pvlib by
+# setting ``min_cos_zenith`` and ``max_clearness_index`` which each have
+# reasonable defaults, but there are still concerning spikes at sunset for Jan.
+# 5th & 7th, April 4th, 5th, & 7th, and July 6th & 7th. The DISC & DIRINT
+# methods differ from Erbs and Boland be including airmass, which seems to
+# reduce DNI spikes over 1000[W/m^2], but still have errors at other times.
+#
+# Another difference is that DISC & DIRINT return DNI whereas Erbs & Boland
+# calculate the diffuse fraction which is then used to derive DNI from GHI and
+# the solar zenith, which exacerbates errors at low sun elevation due to the
+# relation: DNI = GHI*(1 - DF)/cos(zenith).

docs/sphinx/source/api.rst

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

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

@@ -96,6 +96,10 @@ Enhancements
 * :py:meth:`~pvlib.pvsystem.PVSystem.get_ac` is added to calculate AC power
   from DC power. Use parameter ``model`` to specify which inverter model to use.
   (:pull:`1147`, :issue:`998`, :pull:`1150`)
+* :py:meth:`~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

@@ -2186,6 +2186,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 ratio of GHI to extraterrestrial irradiance or clearness index, kt.
+
+    .. math::
+
+        \mathit{DF} = \frac{1}{1 + \exp\left(-5 + 8.6 k_t\right)}
+
+    where :math:`k_t` is the clearness index.
+
+    Parameters
+    ----------
+    ghi: numeric
+        Global horizontal irradiance in W/m^2.
+    zenith: numeric
+        True (not refraction-corrected) zenith angles in decimal degrees.
+    datetime_or_doy : int, float, array, pd.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] 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`
+    .. [2] 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`
+
+    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

@@ -636,6 +636,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