From 229a187265d8cdd59a39dc631b1236890b110119 Mon Sep 17 00:00:00 2001 From: Echedey Luis <80125792+echedey-ls@users.noreply.github.com> Date: Thu, 20 Jun 2024 17:53:09 +0200 Subject: [PATCH] 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 :knife: * 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 * 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 --------- 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 Co-authored-by: Cliff Hansen <5393711+cwhanse@users.noreply.github.com> --- docs/examples/agrivoltaics/README.rst | 2 + .../plot_diffuse_PAR_Spitters_relationship.py | 161 ++++++++++++++++++ .../reference/irradiance/components.rst | 1 + docs/sphinx/source/whatsnew/v0.11.0.rst | 4 + pvlib/irradiance.py | 73 ++++++++ pvlib/tests/test_irradiance.py | 20 +++ 6 files changed, 261 insertions(+) create mode 100644 docs/examples/agrivoltaics/README.rst create mode 100644 docs/examples/agrivoltaics/plot_diffuse_PAR_Spitters_relationship.py diff --git a/docs/examples/agrivoltaics/README.rst b/docs/examples/agrivoltaics/README.rst new file mode 100644 index 0000000000..9c3c2ad045 --- /dev/null +++ b/docs/examples/agrivoltaics/README.rst @@ -0,0 +1,2 @@ +Agrivoltaic Systems Modelling +----------------------------- diff --git a/docs/examples/agrivoltaics/plot_diffuse_PAR_Spitters_relationship.py b/docs/examples/agrivoltaics/plot_diffuse_PAR_Spitters_relationship.py new file mode 100644 index 0000000000..97ebe860a6 --- /dev/null +++ b/docs/examples/agrivoltaics/plot_diffuse_PAR_Spitters_relationship.py @@ -0,0 +1,161 @@ +""" +Calculating daily diffuse PAR using Spitter's relationship +========================================================== + +This example demonstrates how to calculate the diffuse photosynthetically +active radiation (PAR) from diffuse fraction of broadband insolation. +""" + +# %% +# The photosynthetically active radiation (PAR) is a key metric in quantifying +# the photosynthesis process of plants. As with broadband irradiance, PAR can +# be divided into direct and diffuse components. The diffuse fraction of PAR +# with respect to the total PAR is important in agrivoltaic systems, where +# crops are grown under solar panels. The diffuse fraction of PAR can be +# calculated using the Spitter's relationship [1]_ implemented in +# :py:func:`~pvlib.irradiance.diffuse_par_spitters`. +# This model requires the average daily solar zenith angle and the +# daily fraction of the broadband insolation that is diffuse as inputs. +# +# .. note:: +# Understanding the distinction between the broadband insolation and the PAR +# is a key concept. Broadband insolation is the total amount of solar +# energy that gets to a surface, often used in PV applications, while PAR +# is a measurement of a narrower spectrum of wavelengths that are involved +# in photosynthesis. See section on *Photosynthetically Active insolation* +# in pp. 222-223 of [1]_. +# +# References +# ---------- +# .. [1] C. J. T. Spitters, H. A. J. M. Toussaint, and J. Goudriaan, +# 'Separating the diffuse and direct component of global radiation and its +# implications for modeling canopy photosynthesis Part I. Components of +# incoming radiation', Agricultural and Forest Meteorology, vol. 38, +# no. 1, pp. 217-229, Oct. 1986, :doi:`10.1016/0168-1923(86)90060-2`. +# +# Read some example data +# ^^^^^^^^^^^^^^^^^^^^^^ +# Let's read some weather data from a TMY3 file and calculate the solar +# position. + +import pvlib +import pandas as pd +import matplotlib.pyplot as plt +from matplotlib.dates import AutoDateLocator, ConciseDateFormatter +from pathlib import Path + +# Datafile found in the pvlib distribution +DATA_FILE = Path(pvlib.__path__[0]).joinpath("data", "723170TYA.CSV") + +tmy, metadata = pvlib.iotools.read_tmy3( + DATA_FILE, coerce_year=2002, map_variables=True +) +tmy = tmy.filter( + ["ghi", "dhi", "dni", "pressure", "temp_air"] +) # remaining columns are not needed +tmy = tmy["2002-09-06":"2002-09-21"] # select some days + +solar_position = pvlib.solarposition.get_solarposition( + # TMY timestamp is at end of hour, so shift to center of interval + tmy.index.shift(freq="-30T"), + latitude=metadata["latitude"], + longitude=metadata["longitude"], + altitude=metadata["altitude"], + pressure=tmy["pressure"] * 100, # convert from millibar to Pa + temperature=tmy["temp_air"], +) +solar_position.index = tmy.index # reset index to end of the hour + +# %% +# Calculate daily values +# ^^^^^^^^^^^^^^^^^^^^^^ +# The daily average solar zenith angle and the daily diffuse fraction of +# broadband insolation are calculated as follows: + +daily_solar_zenith = solar_position["zenith"].resample("D").mean() +# integration over the day with a time step of 1 hour +daily_tmy = tmy[["ghi", "dhi"]].resample("D").sum() * 1 +daily_tmy["diffuse_fraction"] = daily_tmy["dhi"] / daily_tmy["ghi"] + +# %% +# Calculate Photosynthetically Active Radiation +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# The total PAR can be approximated as 0.50 times the broadband horizontal +# insolation (integral of GHI) for an average solar elevation higher that 10°. +# See section on *Photosynthetically Active Radiation* in pp. 222-223 of [1]_. + +par = pd.DataFrame({"total": 0.50 * daily_tmy["ghi"]}, index=daily_tmy.index) +if daily_solar_zenith.min() < 10: + raise ValueError( + "The total PAR can't be assumed to be half the broadband insolation " + + "for average zenith angles lower than 10°." + ) + +# Calculate broadband insolation diffuse fraction, input of the Spitter's model +daily_tmy["diffuse_fraction"] = daily_tmy["dhi"] / daily_tmy["ghi"] + +# Calculate diffuse PAR fraction using Spitter's relationship +par["diffuse_fraction"] = pvlib.irradiance.diffuse_par_spitters( + solar_position["zenith"], daily_tmy["diffuse_fraction"] +) + +# Finally, calculate the diffuse PAR +par["diffuse"] = par["total"] * par["diffuse_fraction"] + +# %% +# Plot the results +# ^^^^^^^^^^^^^^^^ +# Insolation on left axis, diffuse fraction on right axis + +fig, ax_l = plt.subplots(figsize=(12, 6)) +ax_l.set( + xlabel="Time", + ylabel="Daily insolation $[Wh/m^2/day]$", + title="Diffuse PAR using Spitter's relationship", +) +ax_l.xaxis.set_major_formatter( + ConciseDateFormatter(AutoDateLocator(), tz=daily_tmy.index.tz) +) +ax_l.plot( + daily_tmy.index, + daily_tmy["ghi"], + label="Broadband: total", + color="deepskyblue", +) +ax_l.plot( + daily_tmy.index, + daily_tmy["dhi"], + label="Broadband: diffuse", + color="skyblue", + linestyle="-.", +) +ax_l.plot(daily_tmy.index, par["total"], label="PAR: total", color="orangered") +ax_l.plot( + daily_tmy.index, + par["diffuse"], + label="PAR: diffuse", + color="coral", + linestyle="-.", +) +ax_l.grid() +ax_l.legend(loc="upper left") + +ax_r = ax_l.twinx() +ax_r.set(ylabel="Diffuse fraction") +ax_r.plot( + daily_tmy.index, + daily_tmy["diffuse_fraction"], + label="Broadband diffuse fraction", + color="plum", + linestyle=":", +) +ax_r.plot( + daily_tmy.index, + par["diffuse_fraction"], + label="PAR diffuse fraction", + color="chocolate", + linestyle=":", +) +ax_r.legend(loc="upper right") + +plt.show() diff --git a/docs/sphinx/source/reference/irradiance/components.rst b/docs/sphinx/source/reference/irradiance/components.rst index 70e94ab38d..ce75d9d083 100644 --- a/docs/sphinx/source/reference/irradiance/components.rst +++ b/docs/sphinx/source/reference/irradiance/components.rst @@ -14,3 +14,4 @@ Decomposing and combining irradiance irradiance.get_ground_diffuse irradiance.dni irradiance.complete_irradiance + irradiance.diffuse_par_spitters diff --git a/docs/sphinx/source/whatsnew/v0.11.0.rst b/docs/sphinx/source/whatsnew/v0.11.0.rst index 5806df2cf2..a404fe9908 100644 --- a/docs/sphinx/source/whatsnew/v0.11.0.rst +++ b/docs/sphinx/source/whatsnew/v0.11.0.rst @@ -52,6 +52,10 @@ Enhancements * Added extraterrestrial and direct spectra of the ASTM G173-03 standard with the new function :py:func:`pvlib.spectrum.get_reference_spectra`. (:issue:`1963`, :pull:`2039`) +* Add function :py:func:`pvlib.irradiance.diffuse_par_spitters` to calculate the + diffuse fraction of Photosynthetically Active Radiation (PAR) from the + global diffuse fraction and the solar zenith. + (:issue:`2047`, :pull:`2048`) Bug fixes ~~~~~~~~~ diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py index 73448b6d6d..1f16bdf854 100644 --- a/pvlib/irradiance.py +++ b/pvlib/irradiance.py @@ -3761,3 +3761,76 @@ def louche(ghi, solar_zenith, datetime_or_doy, max_zenith=90): data = pd.DataFrame(data, index=datetime_or_doy) return data + + +def diffuse_par_spitters(daily_solar_zenith, global_diffuse_fraction): + r""" + Derive daily diffuse fraction of Photosynthetically Active Radiation (PAR) + from daily average solar zenith and diffuse fraction of daily insolation. + + The relationship is based on the work of Spitters et al. (1986) [1]_. The + resulting value is the fraction of daily PAR that is diffuse. + + .. note:: + The diffuse fraction is defined as the ratio of + diffuse to global daily insolation, in J m⁻² day⁻¹ or equivalent. + + Parameters + ---------- + daily_solar_zenith : numeric + Average daily solar zenith angle. In degrees [°]. + + global_diffuse_fraction : numeric + Fraction of daily global broadband insolation that is diffuse. + Unitless [0, 1]. + + Returns + ------- + par_diffuse_fraction : numeric + Fraction of daily photosynthetically active radiation (PAR) that is + diffuse. Unitless [0, 1]. + + Notes + ----- + The relationship is given by equations (9) & (10) in [1]_ and (1) in [2]_: + + .. math:: + + k_{diffuse\_PAR}^{model} = \frac{PAR_{diffuse}}{PAR_{total}} = + \frac{\left[1 + 0.3 \left(1 - \left(k_d\right) ^2\right)\right] + k_d} + {1 + \left(1 - \left(k_d\right)^2\right) \cos ^2 (90 - \beta) + \cos ^3 \beta} + + where :math:`k_d` is the diffuse fraction of the global insolation, and + :math:`\beta` is the daily average of the solar elevation angle in degrees. + + A comparison using different models for the diffuse fraction of + the global insolation can be found in [2]_ in the context of Sweden. + + References + ---------- + .. [1] C. J. T. Spitters, H. A. J. M. Toussaint, and J. Goudriaan, + 'Separating the diffuse and direct component of global radiation and its + implications for modeling canopy photosynthesis Part I. Components of + incoming radiation', Agricultural and Forest Meteorology, vol. 38, + no. 1, pp. 217-229, Oct. 1986, :doi:`10.1016/0168-1923(86)90060-2`. + .. [2] S. Ma Lu et al., 'Photosynthetically active radiation decomposition + models for agrivoltaic systems applications', Solar Energy, vol. 244, + pp. 536-549, Sep. 2022, :doi:`10.1016/j.solener.2022.05.046`. + """ + # notation change: + # cosd(90-x) = sind(x) and 90-solar_elevation = solar_zenith + cosd_solar_zenith = tools.cosd(daily_solar_zenith) + cosd_solar_elevation = tools.sind(daily_solar_zenith) + par_diffuse_fraction = ( + (1 + 0.3 * (1 - global_diffuse_fraction**2)) + * global_diffuse_fraction + / ( + 1 + + (1 - global_diffuse_fraction**2) + * cosd_solar_zenith**2 + * cosd_solar_elevation**3 + ) + ) + return par_diffuse_fraction diff --git a/pvlib/tests/test_irradiance.py b/pvlib/tests/test_irradiance.py index 19ce0790ee..8e9f99f11f 100644 --- a/pvlib/tests/test_irradiance.py +++ b/pvlib/tests/test_irradiance.py @@ -1419,3 +1419,23 @@ def test_SURFACE_ALBEDOS_deprecated(): @pytest.mark.filterwarnings("ignore:SURFACE_ALBEDOS") def test_SURFACE_ALBEDO_equals(): assert irradiance.SURFACE_ALBEDOS == albedo.SURFACE_ALBEDOS + + +def test_diffuse_par_spitters(): + solar_zenith, global_diffuse_fraction = np.meshgrid( + [90, 85, 75, 60, 40, 30, 10, 0], [0.01, 0.1, 0.3, 0.6, 0.8, 0.99] + ) + solar_zenith = solar_zenith.ravel() + global_diffuse_fraction = global_diffuse_fraction.ravel() + result = irradiance.diffuse_par_spitters( + solar_zenith, global_diffuse_fraction + ) + expected = np.array([ + 0.01300, 0.01290, 0.01226, 0.01118, 0.01125, 0.01189, 0.01293, 0.01300, + 0.12970, 0.12874, 0.12239, 0.11174, 0.11236, 0.11868, 0.12905, 0.12970, + 0.38190, 0.37931, 0.36201, 0.33273, 0.33446, 0.35188, 0.38014, 0.38190, + 0.71520, 0.71178, 0.68859, 0.64787, 0.65033, 0.67472, 0.71288, 0.71520, + 0.88640, 0.88401, 0.86755, 0.83745, 0.83931, 0.85746, 0.88478, 0.88640, + 0.99591, 0.99576, 0.99472, 0.99270, 0.99283, 0.99406, 0.99581, 0.99591, + ]) # fmt: skip + assert_allclose(result, expected, atol=1e-5)