Skip to content

Commit bbc863d

Browse files
mikofskicwhansekandersolarAdamRJensen
authored
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]>
1 parent bee329b commit bbc863d

File tree

7 files changed

+372
-6
lines changed

7 files changed

+372
-6
lines changed

‎docs/examples/irradiance-decomposition/README.rst‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,3 @@
1+
Irradiance Decomposition
2+
------------------------
3+

‎docs/examples/irradiance-decomposition/plot_diffuse_fraction.py‎

Lines changed: 219 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,219 @@
1+
"""
2+
Diffuse Fraction Estimation
3+
===========================
4+
5+
Comparison of diffuse fraction estimation methods used to derive direct and
6+
diffuse components from measured global horizontal irradiance.
7+
"""
8+
9+
# %%
10+
# This example demonstrates how to use diffuse fraction estimation methods to
11+
# obtain direct and diffuse components from measured global horizontal
12+
# irradiance (GHI). Irradiance sensors such as pyranometers typically only
13+
# measure GHI. pvlib provides several functions that can be used to separate
14+
# GHI into the diffuse and direct components. The separate components are
15+
# needed to estimate the total irradiance on a tilted surface.
16+
17+
import pathlib
18+
from matplotlib import pyplot as plt
19+
import pandas as pd
20+
from pvlib.iotools import read_tmy3
21+
from pvlib.solarposition import get_solarposition
22+
from pvlib import irradiance
23+
import pvlib
24+
25+
# For this example we use the Greensboro, North Carolina, TMY3 file which is
26+
# in the pvlib data directory. TMY3 are made from the median months from years
27+
# of data measured from 1990 to 2010. Therefore we change the timestamps to a
28+
# common year, 1990.
29+
DATA_DIR = pathlib.Path(pvlib.__file__).parent / 'data'
30+
greensboro, metadata = read_tmy3(DATA_DIR / '723170TYA.CSV', coerce_year=1990)
31+
32+
# Many of the diffuse fraction estimation methods require the "true" zenith, so
33+
# we calculate the solar positions for the 1990 at Greensboro, NC.
34+
# NOTE: TMY3 files timestamps indicate the end of the hour, so shift indices
35+
# back 30-minutes to calculate solar position at center of the interval
36+
solpos = get_solarposition(
37+
greensboro.index.shift(freq="-30T"), latitude=metadata['latitude'],
38+
longitude=metadata['longitude'], altitude=metadata['altitude'],
39+
pressure=greensboro.Pressure*100, # convert from millibar to Pa
40+
temperature=greensboro.DryBulb)
41+
solpos.index = greensboro.index # reset index to end of the hour
42+
43+
# %%
44+
# pvlib Decomposition Functions
45+
# -----------------------------
46+
# Methods for separating DHI into diffuse and direct components include:
47+
# `DISC`_, `DIRINT`_, `Erbs`_, and `Boland`_.
48+
49+
# %%
50+
# DISC
51+
# ----
52+
#
53+
# DISC :py:func:`~pvlib.irradiance.disc` is an empirical correlation developed
54+
# at SERI (now NREL) in 1987. The direct normal irradiance (DNI) is related to
55+
# clearness index (kt) by two polynomials split at kt = 0.6, then combined with
56+
# an exponential relation with airmass.
57+
58+
out_disc = irradiance.disc(
59+
greensboro.GHI, solpos.zenith, greensboro.index, greensboro.Pressure*100)
60+
# use "complete sum" AKA "closure" equations: DHI = GHI - DNI * cos(zenith)
61+
df_disc = irradiance.complete_irradiance(
62+
solar_zenith=solpos.apparent_zenith, ghi=greensboro.GHI, dni=out_disc.dni,
63+
dhi=None)
64+
out_disc = out_disc.rename(columns={'dni': 'dni_disc'})
65+
out_disc['dhi_disc'] = df_disc.dhi
66+
67+
# %%
68+
# DIRINT
69+
# ------
70+
#
71+
# DIRINT :py:func:`~pvlib.irradiance.dirint` is a modification of DISC
72+
# developed by Richard Perez and Pierre Ineichen in 1992.
73+
74+
dni_dirint = irradiance.dirint(
75+
greensboro.GHI, solpos.zenith, greensboro.index, greensboro.Pressure*100,
76+
temp_dew=greensboro.DewPoint)
77+
# use "complete sum" AKA "closure" equation: DHI = GHI - DNI * cos(zenith)
78+
df_dirint = irradiance.complete_irradiance(
79+
solar_zenith=solpos.apparent_zenith, ghi=greensboro.GHI, dni=dni_dirint,
80+
dhi=None)
81+
out_dirint = pd.DataFrame(
82+
{'dni_dirint': dni_dirint, 'dhi_dirint': df_dirint.dhi},
83+
index=greensboro.index)
84+
85+
# %%
86+
# Erbs
87+
# ----
88+
#
89+
# The Erbs method, :py:func:`~pvlib.irradiance.erbs` developed by Daryl Gregory
90+
# Erbs at the University of Wisconsin in 1982 is a piecewise correlation that
91+
# splits kt into 3 regions: linear for kt <= 0.22, a 4th order polynomial
92+
# between 0.22 < kt <= 0.8, and a horizontal line for kt > 0.8.
93+
94+
out_erbs = irradiance.erbs(greensboro.GHI, solpos.zenith, greensboro.index)
95+
out_erbs = out_erbs.rename(columns={'dni': 'dni_erbs', 'dhi': 'dhi_erbs'})
96+
97+
# %%
98+
# Boland
99+
# ------
100+
#
101+
# The Boland method, :py:func:`~pvlib.irradiance.boland` is a single logistic
102+
# exponential correlation that is continuously differentiable and bounded
103+
# between zero and one.
104+
105+
out_boland = irradiance.boland(greensboro.GHI, solpos.zenith, greensboro.index)
106+
out_boland = out_boland.rename(
107+
columns={'dni': 'dni_boland', 'dhi': 'dhi_boland'})
108+
109+
# %%
110+
# Comparison Plots
111+
# ----------------
112+
# In the plots below we compare the four decomposition models to the TMY3 file
113+
# for Greensboro, North Carolina. We also compare the clearness index, kt, with
114+
# GHI normalized by a reference irradiance, E0 = 1000 [W/m^2], to highlight
115+
# spikes caused when cosine of zenith approaches zero, particularly at sunset.
116+
#
117+
# First we combine the dataframes for the decomposition models and the TMY3
118+
# file together to make plotting easier.
119+
120+
dni_renames = {
121+
'DNI': 'TMY3', 'dni_disc': 'DISC', 'dni_dirint': 'DIRINT',
122+
'dni_erbs': 'Erbs', 'dni_boland': 'Boland'}
123+
dni = [
124+
greensboro.DNI, out_disc.dni_disc, out_dirint.dni_dirint,
125+
out_erbs.dni_erbs, out_boland.dni_boland]
126+
dni = pd.concat(dni, axis=1).rename(columns=dni_renames)
127+
dhi_renames = {
128+
'DHI': 'TMY3', 'dhi_disc': 'DISC', 'dhi_dirint': 'DIRINT',
129+
'dhi_erbs': 'Erbs', 'dhi_boland': 'Boland'}
130+
dhi = [
131+
greensboro.DHI, out_disc.dhi_disc, out_dirint.dhi_dirint,
132+
out_erbs.dhi_erbs, out_boland.dhi_boland]
133+
dhi = pd.concat(dhi, axis=1).rename(columns=dhi_renames)
134+
ghi_kt = pd.concat([greensboro.GHI/1000.0, out_erbs.kt], axis=1)
135+
136+
# %%
137+
# Winter
138+
# ++++++
139+
# Finally, let's plot them for a few winter days and compare
140+
141+
JAN04, JAN07 = '1990-01-04 00:00:00-05:00', '1990-01-07 23:59:59-05:00'
142+
f, ax = plt.subplots(3, 1, figsize=(8, 10), sharex=True)
143+
dni[JAN04:JAN07].plot(ax=ax[0])
144+
ax[0].grid(which="both")
145+
ax[0].set_ylabel('DNI $[W/m^2]$')
146+
ax[0].set_title('Comparison of Diffuse Fraction Estimation Methods')
147+
dhi[JAN04:JAN07].plot(ax=ax[1])
148+
ax[1].grid(which="both")
149+
ax[1].set_ylabel('DHI $[W/m^2]$')
150+
ghi_kt[JAN04:JAN07].plot(ax=ax[2])
151+
ax[2].grid(which='both')
152+
ax[2].set_ylabel(r'$\frac{GHI}{E0}, k_t$')
153+
f.tight_layout()
154+
155+
# %%
156+
# Spring
157+
# ++++++
158+
# And a few spring days ...
159+
160+
APR04, APR07 = '1990-04-04 00:00:00-05:00', '1990-04-07 23:59:59-05:00'
161+
f, ax = plt.subplots(3, 1, figsize=(8, 10), sharex=True)
162+
dni[APR04:APR07].plot(ax=ax[0])
163+
ax[0].grid(which="both")
164+
ax[0].set_ylabel('DNI $[W/m^2]$')
165+
ax[0].set_title('Comparison of Diffuse Fraction Estimation Methods')
166+
dhi[APR04:APR07].plot(ax=ax[1])
167+
ax[1].grid(which="both")
168+
ax[1].set_ylabel('DHI $[W/m^2]$')
169+
ghi_kt[APR04:APR07].plot(ax=ax[2])
170+
ax[2].grid(which='both')
171+
ax[2].set_ylabel(r'$\frac{GHI}{E0}, k_t$')
172+
f.tight_layout()
173+
174+
# %%
175+
# Summer
176+
# ++++++
177+
# And few summer days to finish off the seasons.
178+
179+
JUL04, JUL07 = '1990-07-04 00:00:00-05:00', '1990-07-07 23:59:59-05:00'
180+
f, ax = plt.subplots(3, 1, figsize=(8, 10), sharex=True)
181+
dni[JUL04:JUL07].plot(ax=ax[0])
182+
ax[0].grid(which="both")
183+
ax[0].set_ylabel('DNI $[W/m^2]$')
184+
ax[0].set_title('Comparison of Diffuse Fraction Estimation Methods')
185+
dhi[JUL04:JUL07].plot(ax=ax[1])
186+
ax[1].grid(which="both")
187+
ax[1].set_ylabel('DHI $[W/m^2]$')
188+
ghi_kt[JUL04:JUL07].plot(ax=ax[2])
189+
ax[2].grid(which='both')
190+
ax[2].set_ylabel(r'$\frac{GHI}{E0}, k_t$')
191+
f.tight_layout()
192+
193+
# %%
194+
# Conclusion
195+
# ----------
196+
# This example compares several decomposition models to a TMY3 file for
197+
# Greensboro, North Carolina. However, DNI and DHI in TMY3 files are themselves
198+
# the output of models (either METSTAT or SUNY), and so differences between
199+
# *e.g.* DISC output and the TMY3 file shouldn't be regarded as errors, and
200+
# it's not a reasonable expectation to assume that the four models should
201+
# reproduce the TMY3 values. We refer those interested to the `TMY3`_ and
202+
# `NSRDB`_ user manuals.
203+
#
204+
# The Erbs and Boland models are correlations only based on the clearness index
205+
# kt, which is the ratio of GHI to the the horizontal component of the
206+
# extra-terrestrial irradiance. At low sun elevation (zenith near 90 degrees),
207+
# especially near sunset, kt can explode because the denominator
208+
# (extra-terrestrial irradiance) approaches zero. In pvlib this behavior is
209+
# moderated by ``min_cos_zenith`` and ``max_clearness_index`` which each have
210+
# reasonable defaults. Even so, near sunset there are still spikes in kt and
211+
# DNI from Erbs and Boland for Jan. 5th & 7th, April 4th, 5th, & 7th, and July
212+
# 6th & 7th.
213+
#
214+
# By contrast, the DISC and DIRINT methods estimate DNI first by means of
215+
# correlations, which include additional variables such as airmass. These
216+
# methods seem to reduce DNI spikes over 1000 [W/m^2].
217+
#
218+
# .. _TMY3: https://www.nrel.gov/docs/fy08osti/43156.pdf
219+
# .. _NSRDB: https://www.nrel.gov/docs/fy12osti/54824.pdf

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

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,5 +12,6 @@ DNI estimation models
1212
irradiance.dirint
1313
irradiance.dirindex
1414
irradiance.erbs
15+
irradiance.boland
1516
irradiance.campbell_norman
1617
irradiance.gti_dirint

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

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,10 @@ Enhancements
2020
:py:func:`pvlib.snow.loss_townsend` (:issue:`1636`, :pull:`1653`)
2121
* Added optional ``n_ar`` parameter to :py:func:`pvlib.iam.physical` to
2222
support an anti-reflective coating. (:issue:`1501`, :pull:`1616`)
23+
* :py:func:`~pvlib.irradiance.boland` is another diffuse fraction, DF,
24+
estimation method similar to Erbs but uses a single logistic exponential
25+
correlation between DF and clearness index, kt, that is continuously
26+
differentiable and bounded between zero and one. (:pull:`1179`)
2327
* Add ``model='gueymard2003'``, the airmass model used for REST and REST2,
2428
to :py:func:`~pvlib.atmosphere.get_relative_airmass`. (:pull:`1655`)
2529

‎pvlib/irradiance.py‎

Lines changed: 104 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2265,6 +2265,110 @@ def erbs(ghi, zenith, datetime_or_doy, min_cos_zenith=0.065, max_zenith=87):
22652265
return data
22662266

22672267

2268+
def boland(ghi, solar_zenith, datetime_or_doy, a_coeff=8.645, b_coeff=0.613,
2269+
min_cos_zenith=0.065, max_zenith=87):
2270+
r"""
2271+
Estimate DNI and DHI from GHI using the Boland clearness index model.
2272+
2273+
The Boland model [1]_, [2]_ estimates the diffuse fraction, DF, from global
2274+
horizontal irradiance, GHI, through an empirical relationship between DF
2275+
and the clearness index, :math:`k_t`, the ratio of GHI to horizontal
2276+
extraterrestrial irradiance.
2277+
2278+
.. math::
2279+
2280+
\mathit{DF} = \frac{1}{1 + \exp\left(a \left(k_t - b\right)\right)}
2281+
2282+
2283+
Parameters
2284+
----------
2285+
ghi: numeric
2286+
Global horizontal irradiance. [W/m^2]
2287+
solar_zenith: numeric
2288+
True (not refraction-corrected) zenith angles in decimal degrees.
2289+
datetime_or_doy : numeric, pandas.DatetimeIndex
2290+
Day of year or array of days of year e.g.
2291+
pd.DatetimeIndex.dayofyear, or pd.DatetimeIndex.
2292+
a_coeff : float, default 8.645
2293+
Logistic curve fit coefficient.
2294+
b_coeff : float, default 0.613
2295+
Logistic curve fit coefficient.
2296+
min_cos_zenith : numeric, default 0.065
2297+
Minimum value of cos(zenith) to allow when calculating global
2298+
clearness index :math:`k_t`. Equivalent to zenith = 86.273 degrees.
2299+
max_zenith : numeric, default 87
2300+
Maximum value of zenith to allow in DNI calculation. DNI will be
2301+
set to 0 for times with zenith values greater than `max_zenith`.
2302+
2303+
Returns
2304+
-------
2305+
data : OrderedDict or DataFrame
2306+
Contains the following keys/columns:
2307+
2308+
* ``dni``: the modeled direct normal irradiance in W/m^2.
2309+
* ``dhi``: the modeled diffuse horizontal irradiance in
2310+
W/m^2.
2311+
* ``kt``: Ratio of global to extraterrestrial irradiance
2312+
on a horizontal plane.
2313+
2314+
References
2315+
----------
2316+
.. [1] J. Boland, B. Ridley (2008) Models of Diffuse Solar Fraction. In:
2317+
Badescu V. (eds) Modeling Solar Radiation at the Earth’s Surface.
2318+
Springer, Berlin, Heidelberg. :doi:`10.1007/978-3-540-77455-6_8`
2319+
.. [2] John Boland, Lynne Scott, and Mark Luther, Modelling the diffuse
2320+
fraction of global solar radiation on a horizontal surface,
2321+
Environmetrics 12(2), pp 103-116, 2001,
2322+
:doi:`10.1002/1099-095X(200103)12:2%3C103::AID-ENV447%3E3.0.CO;2-2`
2323+
2324+
See also
2325+
--------
2326+
dirint
2327+
disc
2328+
erbs
2329+
2330+
Notes
2331+
-----
2332+
Boland diffuse fraction differs from other decomposition algorithms by use
2333+
of a logistic function to fit the entire range of clearness index,
2334+
:math:`k_t`. Parameters ``a_coeff`` and ``b_coeff`` are reported in [2]_
2335+
for different time intervals:
2336+
2337+
* 15-minute: ``a = 8.645`` and ``b = 0.613``
2338+
* 1-hour: ``a = 7.997`` and ``b = 0.586``
2339+
"""
2340+
2341+
dni_extra = get_extra_radiation(datetime_or_doy)
2342+
2343+
kt = clearness_index(
2344+
ghi, solar_zenith, dni_extra, min_cos_zenith=min_cos_zenith,
2345+
max_clearness_index=1)
2346+
2347+
# Boland equation
2348+
df = 1.0 / (1.0 + np.exp(a_coeff * (kt - b_coeff)))
2349+
# NOTE: [2] has different coefficients, for different time intervals
2350+
# 15-min: df = 1 / (1 + exp(8.645 * (kt - 0.613)))
2351+
# 1-hour: df = 1 / (1 + exp(7.997 * (kt - 0.586)))
2352+
2353+
dhi = df * ghi
2354+
2355+
dni = (ghi - dhi) / tools.cosd(solar_zenith)
2356+
bad_values = (solar_zenith > max_zenith) | (ghi < 0) | (dni < 0)
2357+
dni = np.where(bad_values, 0, dni)
2358+
# ensure that closure relationship remains valid
2359+
dhi = np.where(bad_values, ghi, dhi)
2360+
2361+
data = OrderedDict()
2362+
data['dni'] = dni
2363+
data['dhi'] = dhi
2364+
data['kt'] = kt
2365+
2366+
if isinstance(datetime_or_doy, pd.DatetimeIndex):
2367+
data = pd.DataFrame(data, index=datetime_or_doy)
2368+
2369+
return data
2370+
2371+
22682372
def campbell_norman(zenith, transmittance, pressure=101325.0,
22692373
dni_extra=1367.0):
22702374
'''

‎pvlib/tests/test_irradiance.py‎

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -804,6 +804,22 @@ def test_erbs():
804804
assert_frame_equal(np.round(out, 0), np.round(expected, 0))
805805

806806

807+
def test_boland():
808+
index = pd.DatetimeIndex(['20190101']*3 + ['20190620'])
809+
ghi = pd.Series([0, 50, 1000, 1000], index=index)
810+
zenith = pd.Series([120, 85, 10, 10], index=index)
811+
expected = pd.DataFrame(np.array(
812+
[[0.0, 0.0, 0.0],
813+
[81.9448546, 42.8580353, 0.405723511],
814+
[723.764990, 287.230626, 0.718132729],
815+
[805.020419, 207.209650, 0.768214312]]),
816+
columns=['dni', 'dhi', 'kt'], index=index)
817+
818+
out = irradiance.boland(ghi, zenith, index)
819+
820+
assert np.allclose(out, expected)
821+
822+
807823
def test_erbs_min_cos_zenith_max_zenith():
808824
# map out behavior under difficult conditions with various
809825
# limiting kwargs settings

0 commit comments

Comments
 (0)