diff --git a/docs/examples/irradiance-decomposition/README.rst b/docs/examples/irradiance-decomposition/README.rst new file mode 100644 index 0000000000..95f18c9330 --- /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 0000000000..dbd0406cf6 --- /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 2b89a272d7..f0d1495889 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 fc2afa3321..c5d246310f 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 f6d166d44c..79ce3c3904 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/pvlib/tests/test_irradiance.py b/pvlib/tests/test_irradiance.py index ae3a7ec88a..cffdd23e40 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 diff --git a/readthedocs.yml b/readthedocs.yml index fb2d1374bb..dde255335c 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