Repository: pvlib/pvlib-python
Pull request #1179 (2021-02-26 04:41:51)

add Boland DF estimation
<!-- Thank you for your contribution! The following items must be addressed before the code can be merged. Please don't hesitate to ask for help if you're unsure of how to accomplish any of the items. Feel free to remove checklist items that are not relevant to your change. -->

 - [x] referenced in Google group: https://groups.google.com/g/pvlib-python/c/TsDIiU19AUU/m/Bwi9595PAwAJ
 - [x] I am familiar with the [contributing guidelines](https://pvlib-python.readthedocs.io/en/latest/contributing.html)
 - [x] Tests added
 - [x] Updates entries to [`docs/sphinx/source/api.rst`](https://github.com/pvlib/pvlib-python/blob/master/docs/sphinx/source/api.rst) for API changes.
 - [x] Adds description and name entries in the appropriate "what's new" file in [`docs/sphinx/source/whatsnew`](https://github.com/pvlib/pvlib-python/tree/master/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` ``).
 - [x] New code is fully documented. Includes [numpydoc](https://numpydoc.readthedocs.io/en/latest/format.html) compliant docstrings, examples, and comments where necessary.
 - [x] Pull request is nearly complete and ready for detailed review.
 - [x] Maintainer: Appropriate GitHub Labels and Milestone are assigned to the Pull Request and linked Issue.

<!-- Brief description of the problem and proposed solution (if not already fully described in the issue linked to above): -->
`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.

mikofski: Hi @cwhanse, interested in reviewing the gallery example and added Boland DNI estimation function? Thanks
mikofski: Hi @cwhanse sorry for the very long delay in responding. I have accepted your changes, added discussions for your suggestions, and responded to your comments. Let me know if you have any additional feedback. Thanks!
mikofski: BTW: I think the readthedocs failure is unrelated, I'm trying to figure it out - the bifacial gallery examples are failing with this error:
>ImportError: cannot import name 'geos_geometrycollection_from_py' from 'shapely.geometry.collection' (/home/docs/checkouts/readthedocs.org/user_builds/pvlib-python/envs/1179/lib/python3.7/site-packages/shapely/geometry/collection.py)
mikofski: looks like we do not have reproducible readthedocs builds - see this: https://docs.readthedocs.io/en/stable/guides/reproducible-builds.html
mikofski: @kanderso-nrel I've figured out why readthedocs is going bonkers. For some reason the version is being read as pvlib-0.9.0a5.dev58+g<commit-SHA> EG: the last commit (9715084):
![image](https://user-images.githubusercontent.com/1385621/221510781-3059b047-3bc0-4bcb-9c68-28e7eeb83c1e.png)

This only started happening after I merged with main 9d9421e, but oddly, that was the last commit that _did_ have the correct version: pvlib-0.9.5.dev16+g9d9421e
![image](https://user-images.githubusercontent.com/1385621/221515109-9e2dd36f-5b1a-4bc2-a7b8-ed39909e9581.png)

And readthedocs _did_ build. This has something to do with the pep517/518 changes, `pyproject.toml` and `importlib_metadata` but none of these files are affected as far as I can tell.

When the version is misread as v0.9.0, that causes pvfactors to downgrade, and that causes shapely>v2 to be installed leading to the build failure.

I don't understand how the new build system works, all I can hope is that it all works out when/if this gets merged.
mikofski: You can view a rendering of the docs here: https://pvlib-python--1179.org.readthedocs.build/en/1179/

I guess the readthedocs failures has something to do with [`setuptools_scm`](https://github.com/pypa/setuptools_scm). I had to delete the `pvlib.egg-info` file in my local working copy to get it to return the correct version:
```python
from importlib_metadata import version
import pvlib
version('pvlib')
# '0.9.5.dev34+g1ae35cd'
```

Before I purged the egg-info file, it was returning 0.9.0!
kanderso-nrel: @mikofski I think the problem was that RTD does a shallow clone by default and didn't give `setuptools_scm` enough info to determine the correct version string for this PR (which has a long/strange pedigree now?  not sure...)

Anyway to get around this problem we can do "unshallow" clones instead (see [docs](https://docs.readthedocs.io/en/stable/build-customization.html#unshallow-git-clone)).  I hope you don't mind that I went ahead and pushed a fix directly to this branch!
AdamRJensen: I'm out back country skiing but I'd like to review on Monday if it isn't merged by then 😎
mikofski: @kanderso-nrel thanks for tips. I think I've addressed everything. @AdamRJensen looking forward to your comments.
mikofski: @pvlib/pvlib-maintainer ready to merge

Patch:
diff --git a/docs/examples/irradiance-decomposition/README.rst b/docs/examples/irradiance-decomposition/README.rst
new file mode 100644
index 0000000..95f18c9
--- /dev/null
+++ b/docs/examples/irradiance-decomposition/README.rst
@@ -0,0 +1,3 @@
+Irradiance Decomposition
+------------------------
+
diff --git a/docs/examples/irradiance-decomposition/plot_diffuse_fraction.py b/docs/examples/irradiance-decomposition/plot_diffuse_fraction.py
new file mode 100644
index 0000000..dbd0406
--- /dev/null
+++ b/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
diff --git a/docs/sphinx/source/reference/irradiance/decomposition.rst b/docs/sphinx/source/reference/irradiance/decomposition.rst
index 2b89a27..f0d1495 100644
--- a/docs/sphinx/source/reference/irradiance/decomposition.rst
+++ b/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
diff --git a/docs/sphinx/source/whatsnew/v0.9.5.rst b/docs/sphinx/source/whatsnew/v0.9.5.rst
index fc2afa3..c5d2463 100644
--- a/docs/sphinx/source/whatsnew/v0.9.5.rst
+++ b/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`)
 
diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py
index f6d166d..79ce3c3 100644
--- a/pvlib/irradiance.py
+++ b/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):
     '''
diff --git a/readthedocs.yml b/readthedocs.yml
index fb2d137..dde2553 100644
--- a/readthedocs.yml
+++ b/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
\ No newline at end of file
+   # only use the packages specified in setup.py
+   system_packages: false
+
+   install:
+      - method: pip
+        path: .
+        extra_requirements:
+          - doc


Test patch:
diff --git a/pvlib/tests/test_irradiance.py b/pvlib/tests/test_irradiance.py
index ae3a7ec..cffdd23 100644
--- a/pvlib/tests/test_irradiance.py
+++ b/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