forked from pvlib/pvlib-python
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathplot_diffuse_fraction.py
213 lines (190 loc) · 7.88 KB
/
plot_diffuse_fraction.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
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