Skip to content

Commit 229a187

Browse files
echedey-lsAdamRJensenkandersolarcwhanse
authored
Agrivoltaics - PAR diffuse fraction model (#2048)
* New PAR module with spitters_relationship * Update v0.11.0.rst * linter * API update * Example rendering * Update test_par.py * Apply suggestions from code review (Adam) Co-authored-by: Adam R. Jensen <39184289+AdamRJensen@users.noreply.github.com> * Update par.py * Move function to spectrum (mismatch.py) * Improve units formatting * Split legends * Remove api page, move to spectrum index * Update v0.11.0.rst * Move to ``irradiance.py`` * Flake8 🔪 * Fix trigonometry - double testing with a spreadsheet Co-Authored-By: Kevin Anderson <57452607+kandersolar@users.noreply.github.com> * Move section of PAR Co-Authored-By: Kevin Anderson <57452607+kandersolar@users.noreply.github.com> * I should read more carefully Co-Authored-By: Kevin Anderson <57452607+kandersolar@users.noreply.github.com> * Update decomposition.rst * Apply suggestions from code review (Cliff) Co-authored-by: Cliff Hansen <cwhanse@sandia.gov> * More docs refurbishment Co-Authored-By: Cliff Hansen <5393711+cwhanse@users.noreply.github.com> * Rename to `diffuse_par_spitters` Instead of `spitters_relationship` Co-Authored-By: Cliff Hansen <5393711+cwhanse@users.noreply.github.com> * `global` -> `broadband` Co-Authored-By: Cliff Hansen <5393711+cwhanse@users.noreply.github.com> * Code review from Adam, first batch Co-Authored-By: Adam R. Jensen <39184289+AdamRJensen@users.noreply.github.com> * Apply trigonometric property * Fix merge - linter * Forgot to apply this comment Co-Authored-By: Adam R. Jensen <39184289+AdamRJensen@users.noreply.github.com> * Remove model from eq Co-Authored-By: Adam R. Jensen <39184289+AdamRJensen@users.noreply.github.com> * Dailies, insolation instead of instant, irradiance values Co-Authored-By: Adam R. Jensen <39184289+AdamRJensen@users.noreply.github.com> * More docs refurbishment Co-Authored-By: Adam R. Jensen <39184289+AdamRJensen@users.noreply.github.com> * Review from Cliff Co-authored-by: Cliff Hansen <cwhanse@sandia.gov> --------- Co-authored-by: Adam R. Jensen <39184289+AdamRJensen@users.noreply.github.com> Co-authored-by: Kevin Anderson <57452607+kandersolar@users.noreply.github.com> Co-authored-by: Cliff Hansen <cwhanse@sandia.gov> Co-authored-by: Cliff Hansen <5393711+cwhanse@users.noreply.github.com>
1 parent 418d6d0 commit 229a187

File tree

6 files changed

+261
-0
lines changed

6 files changed

+261
-0
lines changed

‎docs/examples/agrivoltaics/README.rst‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,2 @@
1+
Agrivoltaic Systems Modelling
2+
-----------------------------

‎docs/examples/agrivoltaics/plot_diffuse_PAR_Spitters_relationship.py‎

Lines changed: 161 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,161 @@
1+
"""
2+
Calculating daily diffuse PAR using Spitter's relationship
3+
==========================================================
4+
5+
This example demonstrates how to calculate the diffuse photosynthetically
6+
active radiation (PAR) from diffuse fraction of broadband insolation.
7+
"""
8+
9+
# %%
10+
# The photosynthetically active radiation (PAR) is a key metric in quantifying
11+
# the photosynthesis process of plants. As with broadband irradiance, PAR can
12+
# be divided into direct and diffuse components. The diffuse fraction of PAR
13+
# with respect to the total PAR is important in agrivoltaic systems, where
14+
# crops are grown under solar panels. The diffuse fraction of PAR can be
15+
# calculated using the Spitter's relationship [1]_ implemented in
16+
# :py:func:`~pvlib.irradiance.diffuse_par_spitters`.
17+
# This model requires the average daily solar zenith angle and the
18+
# daily fraction of the broadband insolation that is diffuse as inputs.
19+
#
20+
# .. note::
21+
# Understanding the distinction between the broadband insolation and the PAR
22+
# is a key concept. Broadband insolation is the total amount of solar
23+
# energy that gets to a surface, often used in PV applications, while PAR
24+
# is a measurement of a narrower spectrum of wavelengths that are involved
25+
# in photosynthesis. See section on *Photosynthetically Active insolation*
26+
# in pp. 222-223 of [1]_.
27+
#
28+
# References
29+
# ----------
30+
# .. [1] C. J. T. Spitters, H. A. J. M. Toussaint, and J. Goudriaan,
31+
# 'Separating the diffuse and direct component of global radiation and its
32+
# implications for modeling canopy photosynthesis Part I. Components of
33+
# incoming radiation', Agricultural and Forest Meteorology, vol. 38,
34+
# no. 1, pp. 217-229, Oct. 1986, :doi:`10.1016/0168-1923(86)90060-2`.
35+
#
36+
# Read some example data
37+
# ^^^^^^^^^^^^^^^^^^^^^^
38+
# Let's read some weather data from a TMY3 file and calculate the solar
39+
# position.
40+
41+
import pvlib
42+
import pandas as pd
43+
import matplotlib.pyplot as plt
44+
from matplotlib.dates import AutoDateLocator, ConciseDateFormatter
45+
from pathlib import Path
46+
47+
# Datafile found in the pvlib distribution
48+
DATA_FILE = Path(pvlib.__path__[0]).joinpath("data", "723170TYA.CSV")
49+
50+
tmy, metadata = pvlib.iotools.read_tmy3(
51+
DATA_FILE, coerce_year=2002, map_variables=True
52+
)
53+
tmy = tmy.filter(
54+
["ghi", "dhi", "dni", "pressure", "temp_air"]
55+
) # remaining columns are not needed
56+
tmy = tmy["2002-09-06":"2002-09-21"] # select some days
57+
58+
solar_position = pvlib.solarposition.get_solarposition(
59+
# TMY timestamp is at end of hour, so shift to center of interval
60+
tmy.index.shift(freq="-30T"),
61+
latitude=metadata["latitude"],
62+
longitude=metadata["longitude"],
63+
altitude=metadata["altitude"],
64+
pressure=tmy["pressure"] * 100, # convert from millibar to Pa
65+
temperature=tmy["temp_air"],
66+
)
67+
solar_position.index = tmy.index # reset index to end of the hour
68+
69+
# %%
70+
# Calculate daily values
71+
# ^^^^^^^^^^^^^^^^^^^^^^
72+
# The daily average solar zenith angle and the daily diffuse fraction of
73+
# broadband insolation are calculated as follows:
74+
75+
daily_solar_zenith = solar_position["zenith"].resample("D").mean()
76+
# integration over the day with a time step of 1 hour
77+
daily_tmy = tmy[["ghi", "dhi"]].resample("D").sum() * 1
78+
daily_tmy["diffuse_fraction"] = daily_tmy["dhi"] / daily_tmy["ghi"]
79+
80+
# %%
81+
# Calculate Photosynthetically Active Radiation
82+
# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
83+
# The total PAR can be approximated as 0.50 times the broadband horizontal
84+
# insolation (integral of GHI) for an average solar elevation higher that 10°.
85+
# See section on *Photosynthetically Active Radiation* in pp. 222-223 of [1]_.
86+
87+
par = pd.DataFrame({"total": 0.50 * daily_tmy["ghi"]}, index=daily_tmy.index)
88+
if daily_solar_zenith.min() < 10:
89+
raise ValueError(
90+
"The total PAR can't be assumed to be half the broadband insolation "
91+
+ "for average zenith angles lower than 10°."
92+
)
93+
94+
# Calculate broadband insolation diffuse fraction, input of the Spitter's model
95+
daily_tmy["diffuse_fraction"] = daily_tmy["dhi"] / daily_tmy["ghi"]
96+
97+
# Calculate diffuse PAR fraction using Spitter's relationship
98+
par["diffuse_fraction"] = pvlib.irradiance.diffuse_par_spitters(
99+
solar_position["zenith"], daily_tmy["diffuse_fraction"]
100+
)
101+
102+
# Finally, calculate the diffuse PAR
103+
par["diffuse"] = par["total"] * par["diffuse_fraction"]
104+
105+
# %%
106+
# Plot the results
107+
# ^^^^^^^^^^^^^^^^
108+
# Insolation on left axis, diffuse fraction on right axis
109+
110+
fig, ax_l = plt.subplots(figsize=(12, 6))
111+
ax_l.set(
112+
xlabel="Time",
113+
ylabel="Daily insolation $[Wh/m^2/day]$",
114+
title="Diffuse PAR using Spitter's relationship",
115+
)
116+
ax_l.xaxis.set_major_formatter(
117+
ConciseDateFormatter(AutoDateLocator(), tz=daily_tmy.index.tz)
118+
)
119+
ax_l.plot(
120+
daily_tmy.index,
121+
daily_tmy["ghi"],
122+
label="Broadband: total",
123+
color="deepskyblue",
124+
)
125+
ax_l.plot(
126+
daily_tmy.index,
127+
daily_tmy["dhi"],
128+
label="Broadband: diffuse",
129+
color="skyblue",
130+
linestyle="-.",
131+
)
132+
ax_l.plot(daily_tmy.index, par["total"], label="PAR: total", color="orangered")
133+
ax_l.plot(
134+
daily_tmy.index,
135+
par["diffuse"],
136+
label="PAR: diffuse",
137+
color="coral",
138+
linestyle="-.",
139+
)
140+
ax_l.grid()
141+
ax_l.legend(loc="upper left")
142+
143+
ax_r = ax_l.twinx()
144+
ax_r.set(ylabel="Diffuse fraction")
145+
ax_r.plot(
146+
daily_tmy.index,
147+
daily_tmy["diffuse_fraction"],
148+
label="Broadband diffuse fraction",
149+
color="plum",
150+
linestyle=":",
151+
)
152+
ax_r.plot(
153+
daily_tmy.index,
154+
par["diffuse_fraction"],
155+
label="PAR diffuse fraction",
156+
color="chocolate",
157+
linestyle=":",
158+
)
159+
ax_r.legend(loc="upper right")
160+
161+
plt.show()

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

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -14,3 +14,4 @@ Decomposing and combining irradiance
1414
irradiance.get_ground_diffuse
1515
irradiance.dni
1616
irradiance.complete_irradiance
17+
irradiance.diffuse_par_spitters

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

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -52,6 +52,10 @@ Enhancements
5252
* Added extraterrestrial and direct spectra of the ASTM G173-03 standard with
5353
the new function :py:func:`pvlib.spectrum.get_reference_spectra`.
5454
(:issue:`1963`, :pull:`2039`)
55+
* Add function :py:func:`pvlib.irradiance.diffuse_par_spitters` to calculate the
56+
diffuse fraction of Photosynthetically Active Radiation (PAR) from the
57+
global diffuse fraction and the solar zenith.
58+
(:issue:`2047`, :pull:`2048`)
5559

5660
Bug fixes
5761
~~~~~~~~~

‎pvlib/irradiance.py‎

Lines changed: 73 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3761,3 +3761,76 @@ def louche(ghi, solar_zenith, datetime_or_doy, max_zenith=90):
37613761
data = pd.DataFrame(data, index=datetime_or_doy)
37623762

37633763
return data
3764+
3765+
3766+
def diffuse_par_spitters(daily_solar_zenith, global_diffuse_fraction):
3767+
r"""
3768+
Derive daily diffuse fraction of Photosynthetically Active Radiation (PAR)
3769+
from daily average solar zenith and diffuse fraction of daily insolation.
3770+
3771+
The relationship is based on the work of Spitters et al. (1986) [1]_. The
3772+
resulting value is the fraction of daily PAR that is diffuse.
3773+
3774+
.. note::
3775+
The diffuse fraction is defined as the ratio of
3776+
diffuse to global daily insolation, in J m⁻² day⁻¹ or equivalent.
3777+
3778+
Parameters
3779+
----------
3780+
daily_solar_zenith : numeric
3781+
Average daily solar zenith angle. In degrees [°].
3782+
3783+
global_diffuse_fraction : numeric
3784+
Fraction of daily global broadband insolation that is diffuse.
3785+
Unitless [0, 1].
3786+
3787+
Returns
3788+
-------
3789+
par_diffuse_fraction : numeric
3790+
Fraction of daily photosynthetically active radiation (PAR) that is
3791+
diffuse. Unitless [0, 1].
3792+
3793+
Notes
3794+
-----
3795+
The relationship is given by equations (9) & (10) in [1]_ and (1) in [2]_:
3796+
3797+
.. math::
3798+
3799+
k_{diffuse\_PAR}^{model} = \frac{PAR_{diffuse}}{PAR_{total}} =
3800+
\frac{\left[1 + 0.3 \left(1 - \left(k_d\right) ^2\right)\right]
3801+
k_d}
3802+
{1 + \left(1 - \left(k_d\right)^2\right) \cos ^2 (90 - \beta)
3803+
\cos ^3 \beta}
3804+
3805+
where :math:`k_d` is the diffuse fraction of the global insolation, and
3806+
:math:`\beta` is the daily average of the solar elevation angle in degrees.
3807+
3808+
A comparison using different models for the diffuse fraction of
3809+
the global insolation can be found in [2]_ in the context of Sweden.
3810+
3811+
References
3812+
----------
3813+
.. [1] C. J. T. Spitters, H. A. J. M. Toussaint, and J. Goudriaan,
3814+
'Separating the diffuse and direct component of global radiation and its
3815+
implications for modeling canopy photosynthesis Part I. Components of
3816+
incoming radiation', Agricultural and Forest Meteorology, vol. 38,
3817+
no. 1, pp. 217-229, Oct. 1986, :doi:`10.1016/0168-1923(86)90060-2`.
3818+
.. [2] S. Ma Lu et al., 'Photosynthetically active radiation decomposition
3819+
models for agrivoltaic systems applications', Solar Energy, vol. 244,
3820+
pp. 536-549, Sep. 2022, :doi:`10.1016/j.solener.2022.05.046`.
3821+
"""
3822+
# notation change:
3823+
# cosd(90-x) = sind(x) and 90-solar_elevation = solar_zenith
3824+
cosd_solar_zenith = tools.cosd(daily_solar_zenith)
3825+
cosd_solar_elevation = tools.sind(daily_solar_zenith)
3826+
par_diffuse_fraction = (
3827+
(1 + 0.3 * (1 - global_diffuse_fraction**2))
3828+
* global_diffuse_fraction
3829+
/ (
3830+
1
3831+
+ (1 - global_diffuse_fraction**2)
3832+
* cosd_solar_zenith**2
3833+
* cosd_solar_elevation**3
3834+
)
3835+
)
3836+
return par_diffuse_fraction

‎pvlib/tests/test_irradiance.py‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1419,3 +1419,23 @@ def test_SURFACE_ALBEDOS_deprecated():
14191419
@pytest.mark.filterwarnings("ignore:SURFACE_ALBEDOS")
14201420
def test_SURFACE_ALBEDO_equals():
14211421
assert irradiance.SURFACE_ALBEDOS == albedo.SURFACE_ALBEDOS
1422+
1423+
1424+
def test_diffuse_par_spitters():
1425+
solar_zenith, global_diffuse_fraction = np.meshgrid(
1426+
[90, 85, 75, 60, 40, 30, 10, 0], [0.01, 0.1, 0.3, 0.6, 0.8, 0.99]
1427+
)
1428+
solar_zenith = solar_zenith.ravel()
1429+
global_diffuse_fraction = global_diffuse_fraction.ravel()
1430+
result = irradiance.diffuse_par_spitters(
1431+
solar_zenith, global_diffuse_fraction
1432+
)
1433+
expected = np.array([
1434+
0.01300, 0.01290, 0.01226, 0.01118, 0.01125, 0.01189, 0.01293, 0.01300,
1435+
0.12970, 0.12874, 0.12239, 0.11174, 0.11236, 0.11868, 0.12905, 0.12970,
1436+
0.38190, 0.37931, 0.36201, 0.33273, 0.33446, 0.35188, 0.38014, 0.38190,
1437+
0.71520, 0.71178, 0.68859, 0.64787, 0.65033, 0.67472, 0.71288, 0.71520,
1438+
0.88640, 0.88401, 0.86755, 0.83745, 0.83931, 0.85746, 0.88478, 0.88640,
1439+
0.99591, 0.99576, 0.99472, 0.99270, 0.99283, 0.99406, 0.99581, 0.99591,
1440+
]) # fmt: skip
1441+
assert_allclose(result, expected, atol=1e-5)

0 commit comments

Comments
 (0)