Agrivoltaics - PAR diffuse fraction model

docs/examples/agrivoltaics/README.rst

@@ -0,0 +1,2 @@
+Agrivoltaic Systems Modelling
+-----------------------------

docs/examples/agrivoltaics/plot_diffuse_PAR_Spitters_relationship.py

@@ -0,0 +1,145 @@
+"""
+Calculating the diffuse PAR using Spitter's relationship
+=========================================================
+
+This example demonstrates how to calculate the diffuse photosynthetically
+active radiation (PAR) from total PAR using Spitter's relationship.
+"""
+
+# %%
+# The photosynthetically active radiation (PAR) is a key component in the
+# photosynthesis process of plants. As in photovoltaic systems, 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.spitters_relationship`.
+# This model requires the solar zenith angle and the fraction of the global
+# radiation that is diffuse as inputs.
+#
+# .. note::
+#    Understanding the distinction between the global radiation and the PAR is
+#    a key concept. The global radiation is the total amount of solar radiation
+#    that is usually accounted for in PV applications, while the PAR is a
+#    measurement of a narrower range of wavelengths that are used in
+#    photosynthesis. See section on *Photosynthetically Active Radiation* in
+#    pp. 222-223 of [1]_.
+#
+# The key function used in this example is
+# :py:func:`pvlib.irradiance.spitters_relationship` to calculate the diffuse
+# PAR fraction, as a function of global diffuse fraction and solar zenith.
+#
+# 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
+# ^^^^^^^^^^^^^^^^^^^^^^
+
+import pvlib
+import pandas as pd
+import matplotlib.pyplot as plt
+from matplotlib.dates import DateFormatter
+from pathlib import Path
+
+# Read some sample data
+DATA_FILE = Path(pvlib.__path__[0]).joinpath("data", "723170TYA.CSV")
+
+tmy, metadata = pvlib.iotools.read_tmy3(
+    DATA_FILE, coerce_year=1990, map_variables=True
+)
+tmy = tmy.filter(
+    ["ghi", "dhi", "dni", "pressure", "temp_air"]
+)  # remaining data is not needed
+tmy = tmy[
+    "1990-04-11T06":"1990-04-11T22"
+]  # select a single day for this example
+
+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 Photosynthetically Active Radiation
+# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+# The total PAR can be approximated as 0.50 times the global horizontal
+# irradiance (GHI) for solar elevation higher that 10°.
+# See section on *Photosynthetically Active Radiation* in pp. 222-223 of [1]_.
+
+par = pd.DataFrame({"total": 0.50 * tmy["ghi"]}, index=tmy.index)
+
+# Calculate global irradiance diffuse fraction, input of the Spitter's model
+tmy["diffuse_fraction"] = tmy["dhi"] / tmy["ghi"]
+
+# Calculate diffuse PAR fraction using Spitter's relationship
+par["diffuse_fraction"] = pvlib.irradiance.spitters_relationship(
+    solar_position["zenith"], tmy["diffuse_fraction"]
+)
+
+# Finally, calculate the diffuse PAR
+par["diffuse"] = par["total"] * par["diffuse_fraction"]
+par[solar_position["zenith"] > 80] = (
+    0  # set to zero for elevation < 10 degrees
+)
+
+# %%
+# Plot the results
+# ^^^^^^^^^^^^^^^^
+# Irradiances on left axis, diffuse fractions on right axis
+
+fig, ax_l = plt.subplots(figsize=(12, 6))
+ax_l.set(
+    xlabel="Time",
+    ylabel="Irradiance $[W/m^2]$",
+    title="Diffuse PAR using Spitter's relationship",
+)
+ax_l.xaxis.set_major_formatter(DateFormatter("%H:%M", tz=tmy.index.tz))
+ax_l.plot(tmy.index, tmy["ghi"], label="Global: total", color="deepskyblue")
+ax_l.plot(
+    tmy.index,
+    tmy["dhi"],
+    label="Global: diffuse",
+    color="skyblue",
+    linestyle="-.",
+)
+ax_l.plot(tmy.index, par["total"], label="PAR: total", color="orangered")
+ax_l.plot(
+    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(
+    tmy.index,
+    tmy["diffuse_fraction"],
+    label="Global diffuse fraction",
+    color="plum",
+    linestyle=":",
+)
+ax_r.plot(
+    tmy.index,
+    par["diffuse_fraction"],
+    label="PAR diffuse fraction",
+    color="chocolate",
+    linestyle=":",
+)
+ax_r.legend(loc="upper right")
+
+plt.show()

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

@@ -14,3 +14,12 @@ Decomposing and combining irradiance
    irradiance.get_ground_diffuse
    irradiance.dni
    irradiance.complete_irradiance
+
+
+Photosynthetically Active Radiation
+-----------------------------------
+
+.. autosummary::
+   :toctree: ../generated/
+
+   irradiance.spitters_relationship

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

@@ -19,6 +19,10 @@ Enhancements
   shade perpendicular to ``axis_azimuth``. The function is applicable to both
   fixed-tilt and one-axis tracking systems.
   (:issue:`1689`, :pull:`1725`, :pull:`1962`)
+- Add function :py:func:`pvlib.irradiance.spitters_relationship` 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

pvlib/irradiance.py

@@ -3766,3 +3766,69 @@ def louche(ghi, solar_zenith, datetime_or_doy, max_zenith=90):
         data = pd.DataFrame(data, index=datetime_or_doy)
 
     return data
+
+
+def spitters_relationship(solar_zenith, global_diffuse_fraction):
+    r"""
+    Derive the diffuse fraction of photosynthetically active radiation (PAR)
+    respect to the global radiation diffuse fraction.
+
+    The relationship is based on the work of Spitters et al. (1986) [1]_.
+
+    Parameters
+    ----------
+    solar_zenith : numeric
+        Solar zenith angle. Degrees.
+
+    global_diffuse_fraction : numeric
+        Fraction of the global radiation that is diffuse. Unitless [0, 1].
+
+    Returns
+    -------
+    par_diffuse_fraction : numeric
+        Photosynthetically Active Radiation diffuse fraction. 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^{model}\right) ^2\right)\right]
+        k_d^{model}}
+        {1 + \left(1 - \left(k_d^{model}\right)^2\right) \cos ^2 (90 - \beta)
+        \cos ^3 \beta}
+
+    where :math:`k_d^{model}` is the diffuse fraction of the global radiation,
+    provided by some model.
+
+    A comparison of different models performance for the diffuse fraction of
+    the global irradiance 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
+    sind_solar_zenith = tools.sind(solar_zenith)
+    cosd_solar_elevation = tools.cosd(90 - solar_zenith)
+    par_diffuse_fraction = (
+        (1 + 0.3 * (1 - global_diffuse_fraction**2))
+        * global_diffuse_fraction
+        / (
+            1
+            + (1 - global_diffuse_fraction**2)
+            * sind_solar_zenith**2
+            * cosd_solar_elevation**3
+        )
+    )
+    return par_diffuse_fraction

pvlib/tests/test_irradiance.py

@@ -1406,3 +1406,25 @@ def test_louche():
     out = irradiance.louche(ghi, zenith, index)
 
     assert_frame_equal(out, expected)
+
+
+def test_spitters_relationship():
+    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.spitters_relationship(
+        solar_zenith, global_diffuse_fraction
+    )
+    expected = np.array([
+        0.00650018, 0.00656213, 0.00706211, 0.00874170, 0.01171437, 0.01260581,
+        0.01299765, 0.01299970, 0.06517588, 0.06579393, 0.07077986, 0.08750105,
+        0.11699064, 0.12580782, 0.12967973, 0.12970000, 0.19994764, 0.20176275,
+        0.21635259, 0.26460255, 0.34722693, 0.37134002, 0.38184514, 0.38190000,
+        0.43609756, 0.43933488, 0.46497584, 0.54521789, 0.66826809, 0.70117647,
+        0.71512774, 0.71520000, 0.65176471, 0.65503875, 0.68042968, 0.75414541,
+        0.85271445, 0.87653894, 0.88634962, 0.88640000, 0.97647838, 0.97683827,
+        0.97952006, 0.98634857, 0.99374028, 0.99529135, 0.99590717, 0.9959103
+    ])  # fmt: skip
+    assert_allclose(result, expected, atol=1e-8)