diff --git a/docs/examples/shading/plot_shaded_fraction1d_ns_hsat_example.py b/docs/examples/shading/plot_shaded_fraction1d_ns_hsat_example.py new file mode 100644 index 0000000000..7eae2ed7aa --- /dev/null +++ b/docs/examples/shading/plot_shaded_fraction1d_ns_hsat_example.py @@ -0,0 +1,181 @@ +""" +Shaded fraction of a horizontal single-axis tracker +==================================================== + +This example illustrates how to calculate the 1D shaded fraction of three rows +in a North-South horizontal single axis tracker (HSAT) configuration. +""" + +# %% +# :py:func:`pvlib.shading.shaded_fraction1d` exposes a useful method for the +# calculation of the shaded fraction of the width of a solar collector. Here, +# the width is defined as the dimension perpendicular to the axis of rotation. +# This method for calculating the shaded fraction only requires minor +# modifications to be applicable for fixed-tilt systems. +# +# It is highly recommended to read :py:func:`pvlib.shading.shaded_fraction1d` +# documentation to understand the input parameters and the method's +# capabilities. For more in-depth information, please see the journal paper +# `10.1063/5.0202220 `_ describing +# the methodology. +# +# Let's start by obtaining the true-tracking angles for each of the rows and +# limiting the angles to the range of -50 to 50 degrees. This decision is +# done to create an example scenario with significant shading. +# +# Key functions used in this example are: +# +# 1. :py:func:`pvlib.tracking.singleaxis` to calculate the tracking angles. +# 2. :py:func:`pvlib.shading.projected_solar_zenith_angle` to calculate the +# projected solar zenith angle. +# 3. :py:func:`pvlib.shading.shaded_fraction1d` to calculate the shaded +# fractions. +# +# .. sectionauthor:: Echedey Luis + +import pvlib + +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from matplotlib.dates import DateFormatter + +# Define the solar system parameters +latitude, longitude = 28.51, -13.89 +altitude = pvlib.location.lookup_altitude(latitude, longitude) + +axis_tilt = 3 # degrees, positive is upwards in the axis_azimuth direction +axis_azimuth = 180 # degrees, N-S tracking axis +collector_width = 3.2 # m +pitch = 4.15 # m +gcr = collector_width / pitch +cross_axis_slope = -5 # degrees +surface_to_axis_offset = 0.07 # m + +# Generate a time range for the simulation +times = pd.date_range( + start="2024-01-01T05", + end="2024-01-01T21", + freq="5min", + tz="Atlantic/Canary", +) + +# Calculate the solar position +solar_position = pvlib.solarposition.get_solarposition( + times, latitude, longitude, altitude +) +solar_azimuth = solar_position["azimuth"] +solar_zenith = solar_position["apparent_zenith"] + +# Calculate the tracking angles +rotation_angle = pvlib.tracking.singleaxis( + solar_zenith, + solar_azimuth, + axis_tilt, + axis_azimuth, + max_angle=(-50, 50), # (min, max) degrees + backtrack=False, + gcr=gcr, + cross_axis_tilt=cross_axis_slope, +)["tracker_theta"] + +# %% +# Once the tracker angles have been obtained, the next step is to calculate +# the shaded fraction. Special care must be taken +# to ensure that the shaded or shading tracker roles are correctly assigned +# depending on the solar position. +# This means we will have a result for each row, ``eastmost_shaded_fraction``, +# ``middle_shaded_fraction``, and ``westmost_shaded_fraction``. +# Switching the parameters will be based on the +# sign of :py:func:`pvlib.shading.projected_solar_zenith_angle`. +# +# The following code is verbose to make it easier to understand the process, +# but with some effort you may be able to simplify it. This verbosity also +# allows to change the premises easily per case, e.g., in case of a tracker +# failure or with a different system configuration. + +psza = pvlib.shading.projected_solar_zenith_angle( + solar_zenith, solar_azimuth, axis_tilt, axis_azimuth +) + +# Calculate the shaded fraction for the eastmost row +eastmost_shaded_fraction = np.where( + psza < 0, + 0, # no shaded fraction in the morning + # shaded fraction in the evening + pvlib.shading.shaded_fraction1d( + solar_zenith, + solar_azimuth, + axis_azimuth, + shaded_row_rotation=rotation_angle, + axis_tilt=axis_tilt, + collector_width=collector_width, + pitch=pitch, + surface_to_axis_offset=surface_to_axis_offset, + cross_axis_slope=cross_axis_slope, + shading_row_rotation=rotation_angle, + ), +) + +# Calculate the shaded fraction for the middle row +middle_shaded_fraction = np.where( + psza < 0, + # shaded fraction in the morning + pvlib.shading.shaded_fraction1d( + solar_zenith, + solar_azimuth, + axis_azimuth, + shaded_row_rotation=rotation_angle, + axis_tilt=axis_tilt, + collector_width=collector_width, + pitch=pitch, + surface_to_axis_offset=surface_to_axis_offset, + cross_axis_slope=cross_axis_slope, + shading_row_rotation=rotation_angle, + ), + # shaded fraction in the evening + pvlib.shading.shaded_fraction1d( + solar_zenith, + solar_azimuth, + axis_azimuth, + shaded_row_rotation=rotation_angle, + axis_tilt=axis_tilt, + collector_width=collector_width, + pitch=pitch, + surface_to_axis_offset=surface_to_axis_offset, + cross_axis_slope=cross_axis_slope, + shading_row_rotation=rotation_angle, + ), +) + +# Calculate the shaded fraction for the westmost row +westmost_shaded_fraction = np.where( + psza < 0, + # shaded fraction in the morning + pvlib.shading.shaded_fraction1d( + solar_zenith, + solar_azimuth, + axis_azimuth, + shaded_row_rotation=rotation_angle, + axis_tilt=axis_tilt, + collector_width=collector_width, + pitch=pitch, + surface_to_axis_offset=surface_to_axis_offset, + cross_axis_slope=cross_axis_slope, + shading_row_rotation=rotation_angle, + ), + 0, # no shaded fraction in the evening +) + +# %% +# Plot the shaded fraction result for each row: +plt.plot(times, eastmost_shaded_fraction, label="East-most", color="blue") +plt.plot(times, middle_shaded_fraction, label="Middle", color="green", + linewidth=3, linestyle="--") # fmt: skip +plt.plot(times, westmost_shaded_fraction, label="West-most", color="red") +plt.title(r"$shaded\_fraction1d$ of each row vs time") +plt.xlabel("Time") +plt.gca().xaxis.set_major_formatter(DateFormatter("%H:%M")) +plt.ylabel("Shaded Fraction") +plt.legend() +plt.show() diff --git a/docs/sphinx/source/_images/Anderson_Jensen_2024_Fig3.png b/docs/sphinx/source/_images/Anderson_Jensen_2024_Fig3.png new file mode 100644 index 0000000000..d239d9d4f4 Binary files /dev/null and b/docs/sphinx/source/_images/Anderson_Jensen_2024_Fig3.png differ diff --git a/docs/sphinx/source/reference/effects_on_pv_system_output/shading.rst b/docs/sphinx/source/reference/effects_on_pv_system_output/shading.rst index a0fd74a795..57cfe0b806 100644 --- a/docs/sphinx/source/reference/effects_on_pv_system_output/shading.rst +++ b/docs/sphinx/source/reference/effects_on_pv_system_output/shading.rst @@ -11,3 +11,5 @@ Shading shading.masking_angle_passias shading.sky_diffuse_passias shading.projected_solar_zenith_angle + shading.shaded_fraction1d + diff --git a/docs/sphinx/source/whatsnew/v0.11.0.rst b/docs/sphinx/source/whatsnew/v0.11.0.rst index 8fe205a2a7..b5d543bff1 100644 --- a/docs/sphinx/source/whatsnew/v0.11.0.rst +++ b/docs/sphinx/source/whatsnew/v0.11.0.rst @@ -15,6 +15,10 @@ Deprecations Enhancements ~~~~~~~~~~~~ +* Add function :py:func:`pvlib.shading.shaded_fraction1d`, to calculate the + shade perpendicular to ``axis_azimuth``. The function is applicable to both + fixed-tilt and one-axis tracking systems. + (:issue:`1689`, :pull:`1725`, :pull:`1962`) Bug fixes @@ -36,4 +40,5 @@ Requirements Contributors ~~~~~~~~~~~~ * Cliff Hansen (:ghuser:`cwhanse`) -* Siddharth Kaul (:ghuser:`k10blogger`) \ No newline at end of file +* Mark Mikofski (:ghuser:`mikofski`) +* Siddharth Kaul (:ghuser:`k10blogger`) diff --git a/pvlib/shading.py b/pvlib/shading.py index 007e24e1b7..cbbf80a8f6 100644 --- a/pvlib/shading.py +++ b/pvlib/shading.py @@ -342,3 +342,214 @@ def projected_solar_zenith_angle(solar_zenith, solar_azimuth, # Eq. (5); angle between sun's beam and surface theta_T = np.degrees(np.arctan2(sx_prime, sz_prime)) return theta_T + + +def shaded_fraction1d( + solar_zenith, + solar_azimuth, + axis_azimuth, + shaded_row_rotation, + *, + collector_width, + pitch, + axis_tilt=0, + surface_to_axis_offset=0, + cross_axis_slope=0, + shading_row_rotation=None, +): + r""" + Shaded fraction in the vertical dimension of tilted rows, or perpendicular + to the axis of horizontal rows. + + If ``shading_row_rotation`` isn't provided, it is assumed that + both the shaded row and the shading row (the one blocking the + direct beam) have the same rotation and azimuth values. + + .. warning:: + The function assumes that the roles of the shaded and shading rows + remain the same during the day. In the case where the shading and + shaded rows change throughout the day, e.g. a N-S single-axis tracker, + the inputs must be switched depending on the sign of the projected + solar zenith angle. See the Examples section below. + + .. versionadded:: 0.11.0 + + Parameters + ---------- + solar_zenith : numeric + Solar zenith angle, in degrees. + solar_azimuth : numeric + Solar azimuth angle, in degrees. + axis_azimuth : numeric + Axis azimuth of the rotation axis of a tracker, in degrees. + Fixed-tilt arrays can be considered as a particular case of a tracker. + North=0º, South=180º, East=90º, West=270º. + shaded_row_rotation : numeric + Right-handed rotation of the row receiving the shade, with respect + to ``axis_azimuth``. In degrees :math:`^{\circ}`. + collector_width : numeric + Vertical length of a tilted row. The returned ``shaded_fraction`` + is the ratio of the shadow over this value. + pitch : numeric + Axis-to-axis horizontal spacing of the row. + axis_tilt : numeric, default 0 + Tilt of the rows axis from horizontal. In degrees :math:`^{\circ}`. + surface_to_axis_offset : numeric, default 0 + Distance between the rotating axis and the collector surface. + May be used to account for a torque tube offset. + cross_axis_slope : numeric, default 0 + Angle of the plane containing the rows' axes from + horizontal. Right-handed rotation with respect to the rows axes. + In degrees :math:`^{\circ}`. + shading_row_rotation : numeric, optional + Right-handed rotation of the row casting the shadow, with respect + to the row axis. In degrees :math:`^{\circ}`. + + Returns + ------- + shaded_fraction : numeric + The fraction of the collector width shaded by an adjacent row. A + value of 1 is completely shaded and 0 is no shade. + + Notes + ----- + All length parameters must have the same units. + + Parameters are defined as follow: + + .. figure:: ../../_images/Anderson_Jensen_2024_Fig3.png + :alt: Diagram showing the two rows and the parameters of the model. + + Figure 3 of [1]_. See correspondence between this nomenclature and the + function parameters in the table below. + + +------------------+----------------------------+---------------------+ + | Symbol | Parameter | Units | + +==================+============================+=====================+ + | :math:`\theta_1` | ``shading_row_rotation`` | | + +------------------+----------------------------+ | + | :math:`\theta_2` | ``shaded_row_rotation`` | Degrees | + +------------------+----------------------------+ :math:`^{\circ}` | + | :math:`\beta_c` | ``cross_axis_slope`` | | + +------------------+----------------------------+---------------------+ + | :math:`p` | ``pitch`` | Any consistent | + +------------------+----------------------------+ length unit across | + | :math:`\ell` | ``collector_width`` | all these | + +------------------+----------------------------+ parameters, e.g. | + | :math:`z_0` | ``surface_to_axis_offset`` | :math:`m`. | + +------------------+----------------------------+---------------------+ + | :math:`f_s` | Return value | Dimensionless | + +------------------+----------------------------+---------------------+ + + Examples + -------- + + **Fixed-tilt south-facing array on flat terrain** + + Tilted row with a pitch of 3 m, a collector width of + 2 m, and row rotations of 30°. In the morning. + + >>> shaded_fraction1d(solar_zenith=80, solar_azimuth=104.5, + ... axis_azimuth=90, shaded_row_rotation=30, shading_row_rotation=30, + ... collector_width=2, pitch=3, axis_tilt=0, + ... surface_to_axis_offset=0.05, cross_axis_slope=0) + 0.6827437712114521 + + **Fixed-tilt north-facing array on sloped terrain** + + Tilted row with a pitch of 4 m, a collector width of + 2.5 m, and row rotations of 50° for the shaded + row and 30° for the shading row. The rows are on a + 10° slope, where their axis is on the most inclined + direction (zero cross-axis slope). Shaded in the morning. + + >>> shaded_fraction1d(solar_zenith=65, solar_azimuth=75.5, + ... axis_azimuth=270, shaded_row_rotation=50, shading_row_rotation=30, + ... collector_width=2.5, pitch=4, axis_tilt=10, + ... surface_to_axis_offset=0.05, cross_axis_slope=0) + 0.6975923460352351 + + **N-S single-axis tracker on sloped terrain** + + Horizontal trackers with a pitch of 3 m, a collector width of + 1.4 m, and tracker rotations of 30° pointing east, + in the morning. Terrain slope is 7° west-east (east-most + tracker is higher than the west-most tracker). + + >>> shaded_fraction1d(solar_zenith=50, solar_azimuth=90, axis_azimuth=180, + ... shaded_row_rotation=-30, collector_width=1.4, pitch=3, axis_tilt=0, + ... surface_to_axis_offset=0.10, cross_axis_slope=7) + 0.5828961460616938 + + Note the previous example only is valid for the shaded fraction of the + west-most tracker in the morning, and assuming it is the + shaded tracker during all the day is incorrect. + During the afternoon, it is the one casting the shadow onto the + east-most tracker. + + To calculate the shaded fraction for the east-most + tracker, you must input the corresponding ``shaded_row_rotation`` + in the afternoon. + + >>> shaded_fraction1d(solar_zenith=50, solar_azimuth=270, axis_azimuth=180, + ... shaded_row_rotation=30, collector_width=1.4, pitch=3, axis_tilt=0, + ... surface_to_axis_offset=0.10, cross_axis_slope=7) + 0.4399034444363955 + + You must switch the input/output depending on the + sign of the projected solar zenith angle. See + :py:func:`~pvlib.shading.projected_solar_zenith_angle` and the example + :ref:`sphx_glr_gallery_shading_plot_shaded_fraction1d_ns_hsat_example.py` + + See also + -------- + pvlib.shading.projected_solar_zenith_angle + + References + ---------- + .. [1] Kevin S. Anderson, Adam R. Jensen; Shaded fraction and backtracking + in single-axis trackers on rolling terrain. J. Renewable Sustainable + Energy 1 March 2024; 16 (2): 023504. :doi:`10.1063/5.0202220` + """ + # For nomenclature you may refer to [1]. + + # rotation of row casting the shadow defaults to shaded row's one + if shading_row_rotation is None: + shading_row_rotation = shaded_row_rotation + + # projected solar zenith angle + projected_solar_zenith = projected_solar_zenith_angle( + solar_zenith, + solar_azimuth, + axis_tilt, + axis_azimuth, + ) + + # calculate repeated elements + thetas_1_S_diff = shading_row_rotation - projected_solar_zenith + thetas_2_S_diff = shaded_row_rotation - projected_solar_zenith + thetaS_rotation_diff = projected_solar_zenith - cross_axis_slope + + cos_theta_2_S_diff_abs = np.abs(cosd(thetas_2_S_diff)) + + # Eq. (12) of [1] + t_asterisk = ( + 0.5 + + np.abs(cosd(thetas_1_S_diff)) / cos_theta_2_S_diff_abs / 2 + + ( + np.sign(projected_solar_zenith) + * surface_to_axis_offset + / collector_width + / cos_theta_2_S_diff_abs + * (sind(thetas_2_S_diff) - sind(thetas_1_S_diff)) + ) + - ( + pitch + / collector_width + * cosd(thetaS_rotation_diff) + / cos_theta_2_S_diff_abs + / cosd(cross_axis_slope) + ) + ) + + return np.clip(t_asterisk, 0, 1) diff --git a/pvlib/tests/test_shading.py b/pvlib/tests/test_shading.py index 2f1e5b3410..b8bf8929ae 100644 --- a/pvlib/tests/test_shading.py +++ b/pvlib/tests/test_shading.py @@ -2,11 +2,12 @@ import pandas as pd from pandas.testing import assert_series_equal -from numpy.testing import assert_allclose +from numpy.testing import assert_allclose, assert_approx_equal import pytest from datetime import timezone, timedelta from pvlib import shading +from pvlib.tools import atand @pytest.fixture @@ -223,3 +224,106 @@ def test_projected_solar_zenith_angle_datatypes( ) psz = psz_func(sun_apparent_zenith, axis_azimuth, axis_tilt, axis_azimuth) assert isinstance(psz, cast_type) + + +@pytest.fixture +def sf1d_premises_and_expected(): + """Data comprised of solar position, rows parameters and terrain slope + with respective shade fractions (sf). Returns a 2-tuple with the premises + to be used directly in shaded_fraction1d(...) in the first element and + the expected shaded fractions in the second element. + See [1] in shaded_fraction1d() + Test data sourced from http://doi.org/10.5281/zenodo.10513987 + """ + test_data = pd.DataFrame( + columns=["x_L", "z_L", "theta_L", "x_R", "z_R", "theta_R", "z_0", "l", + "theta_s", "f_s"], + data=( + (1, 0.2, 50, 0, 0, 25, 0, 0.5, 80, 1), + (1, 0.1, 50, 0, 0, 25, 0.05, 0.5, 80, 0.937191), + (1, 0, 50, 0, 0.1, 25, 0, 0.5, 80, 0.30605), + (1, 0, 50, 0, 0.2, 25, 0, 0.5, 80, 0), + (1, 0.2, -25, 0, 0, -50, 0, 0.5, -80, 0), + (1, 0.1, -25, 0, 0, -50, 0, 0.5, -80, 0.30605), + (1, 0, -25, 0, 0.1, -50, 0.1, 0.5, -80, 0.881549), + (1, 0, -25, 0, 0.2, -50, 0, 0.5, -80, 1), + (1, 0.2, 5, 0, 0, 25, 0.05, 0.5, 80, 0.832499), + (1, 0.2, -25, 0, 0, 25, 0.05, 0.5, 80, 0.832499), + (1, 0.2, 5, 0, 0, -45, 0.05, 0.5, 80, 0.832499), + (1, 0.2, -25, 0, 0, -45, 0.05, 0.5, 80, 0.832499), + (1, 0, -25, 0, 0.2, 25, 0.05, 0.5, -80, 0.832499), + (1, 0, -25, 0, 0.2, -5, 0.05, 0.5, -80, 0.832499), + (1, 0, 45, 0, 0.2, 25, 0.05, 0.5, -80, 0.832499), + (1, 0, 45, 0, 0.2, -5, 0.05, 0.5, -80, 0.832499), + ), + ) # fmt: skip + + test_data["cross_axis_slope"] = atand( + (test_data["z_R"] - test_data["z_L"]) + / (test_data["x_L"] - test_data["x_R"]) + ) + test_data["pitch"] = test_data["x_L"] - test_data["x_R"] + # switch Left/Right rows if needed to make the right one the shaded + where_switch = test_data["theta_s"] >= 0 + test_data["theta_L"], test_data["theta_R"] = np.where( + where_switch, + (test_data["theta_L"], test_data["theta_R"]), + (test_data["theta_R"], test_data["theta_L"]), + ) + test_data.rename( + columns={ + "theta_L": "shading_row_rotation", + "theta_R": "shaded_row_rotation", + "z_0": "surface_to_axis_offset", + "l": "collector_width", + "theta_s": "solar_zenith", # for the projected solar zenith angle + "f_s": "shaded_fraction", + }, + inplace=True, + ) + test_data.drop(columns=["x_L", "z_L", "x_R", "z_R"], inplace=True) + # for the projected solar zenith angle + # this returns the same psz angle as test_data["solar_zenith"] + test_data["solar_azimuth"], test_data["axis_azimuth"] = 180, 90 + + # return 1st: premises dataframe first and 2nd: shaded fraction series + return ( + test_data.drop(columns=["shaded_fraction"]), + test_data["shaded_fraction"], + ) + + +def test_shaded_fraction1d(sf1d_premises_and_expected): + """Tests shaded_fraction1d""" + # unwrap sf_premises_and_expected values premises and expected results + premises, expected_sf_array = sf1d_premises_and_expected + # test scalar input + expected_result = expected_sf_array.iloc[0] + sf = shading.shaded_fraction1d(**premises.iloc[0]) + assert_approx_equal(sf, expected_result) + assert isinstance(sf, float) + + # test Series inputs + sf_vec = shading.shaded_fraction1d(**premises) + assert_allclose(sf_vec, expected_sf_array, atol=1e-6) + assert isinstance(sf_vec, pd.Series) + + +def test_shaded_fraction1d_unprovided_shading_row_rotation(): + """Tests shaded_fraction1d without providing shading_row_rotation""" + test_data = pd.DataFrame( + columns=[ + "shaded_row_rotation", "surface_to_axis_offset", "collector_width", + "solar_zenith", "cross_axis_slope", "pitch", "solar_azimuth", + "axis_azimuth", "expected_sf", + ], + data=[ + (30, 0, 5.7735, 60, 0, 5, 90, 180, 0), + (30, 0, 5.7735, 79, 0, 5, 90, 180, 0.5), + (30, 0, 5.7735, 90, 0, 5, 90, 180, 1), + ], + ) # fmt: skip + expected_sf = test_data["expected_sf"] + premises = test_data.drop(columns=["expected_sf"]) + sf = shading.shaded_fraction1d(**premises) + assert_allclose(sf, expected_sf, atol=1e-2)