Projected zenith convenience function (#1904)

* Function prototype

* Update shading.rst

* Update shading.py

* Minimal test

* Implementation

From NREL paper

* Fix, fix, fix, fix & format

* Format issues

* Extend tests (compare with singleaxis) & format with ruff

* Format fixes

* Upgrade tests

* Array -> Axis

* type

* Whatsnew

* xd

* bruh

* Minor Python optimization a la tracking.singleaxis

* Comment and minor optimizations

* Typo found by Mikofski

Reported at: #1725 (comment)

Confirmed via "Slope-Aware Backtracking for Single-Axis Trackers", paragraph after Eq. 1

Co-Authored-By: Mark Mikofski <[email protected]>

* Surface -> Axis

Co-Authored-By: Kevin Anderson <[email protected]>

* Elevation -> Zenith

Co-Authored-By: Kevin Anderson <[email protected]>

* Elev -> Zenith

Co-Authored-By: Kevin Anderson <[email protected]>

* Update shading.py

* Update docstring

Co-Authored-By: Anton Driesse <[email protected]>

* Add comments from `tracking.singleaxis`

Co-Authored-By: Will Holmgren <[email protected]>
Co-Authored-By: Mark Mikofski <[email protected]>

* Singleaxis implementation port & test addition, based on old pvlib.tracking.singleaxis

* Update v0.10.4.rst

* Linter

* Code review

Co-Authored-By: Cliff Hansen <[email protected]>

* Add Fig 5 [1] (still gotta check the built output)

* Add caption, change size and describe in alternate text

* rST fixes ?

* Figures have captions, images do not

https://pandemic-overview.readthedocs.io/en/latest/myGuides/reStructuredText-Images-and-Figures-Examples.html#id18

* Flip arguments order

* I forgot 💀

* Linter are you happy now?

* Remove port test and add edge cases test

Co-Authored-By: Kevin Anderson <[email protected]>

* Update test_shading.py

Co-Authored-By: Kevin Anderson <[email protected]>

* Indentation xd

* Update test_shading.py

* I forgot how to code

* Align data

* Docstring suggestion from Kevin

Co-Authored-By: Kevin Anderson <[email protected]>

* Update link to example?

* Link, please work

* Update shading.py

* Update shading.py

* Update shading.py

* Update shading.py

* Update shading.py

* Update shading.py

* Update shading.py

* Update shading.py

* Lintaaaaaaarrrgh

Fixed the link finally

* Update pvlib/shading.py

Co-authored-by: Kevin Anderson <[email protected]>

---------

Co-authored-by: Mark Mikofski <[email protected]>
Co-authored-by: Kevin Anderson <[email protected]>
Co-authored-by: Anton Driesse <[email protected]>
Co-authored-by: Will Holmgren <[email protected]>
Co-authored-by: Cliff Hansen <[email protected]>
Co-authored-by: Kevin Anderson <[email protected]>

Author: 6 people

docs/sphinx/source/_images/Anderson_Mikofski_2020_Fig5.jpg

@@ -9,4 +9,5 @@ Shading
    shading.ground_angle
    shading.masking_angle
    shading.masking_angle_passias
-   shading.sky_diffuse_passias
+   shading.sky_diffuse_passias
+   shading.projected_solar_zenith_angle

docs/sphinx/source/reference/effects_on_pv_system_output/shading.rst

@@ -8,11 +8,14 @@ v0.10.4 (Anticipated March, 2024)
 Enhancements
 ~~~~~~~~~~~~
 * Added the Huld PV model used by PVGIS (:pull:`1940`)
+* Added function :py:func:`pvlib.shading.projected_solar_zenith_angle`,
+  a common calculation in shading and tracking. (:issue:`1734`, :pull:`1904`)
 * Added :py:func:`~pvlib.iotools.get_solrad` for fetching irradiance data from
   the SOLRAD ground station network. (:pull:`1967`)
 * Added metadata parsing to :py:func:`~pvlib.iotools.read_solrad` to follow the standard iotools
   convention of returning a tuple of (data, meta). Previously the function only returned a dataframe. (:pull:`1968`)
 
+
 Bug fixes
 ~~~~~~~~~
 * Fixed an error in solar position calculations when using

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

@@ -232,3 +232,113 @@ def sky_diffuse_passias(masking_angle):
        Available at https://www.nrel.gov/docs/fy18osti/67399.pdf
     """
     return 1 - cosd(masking_angle/2)**2
+
+
+def projected_solar_zenith_angle(solar_zenith, solar_azimuth,
+                                 axis_tilt, axis_azimuth):
+    r"""
+    Calculate projected solar zenith angle in degrees.
+
+    This solar zenith angle is projected onto the plane whose normal vector is
+    defined by ``axis_tilt`` and ``axis_azimuth``. The normal vector is in the
+    direction of ``axis_azimuth`` (clockwise from north) and tilted from
+    horizontal by ``axis_tilt``. See Figure 5 in [1]_:
+
+    .. figure:: ../../_images/Anderson_Mikofski_2020_Fig5.jpg
+       :alt: Wire diagram of coordinates systems to obtain the projected angle.
+       :align: center
+       :scale: 50 %
+
+       Fig. 5, [1]_: Solar coordinates projection onto tracker rotation plane.
+
+    Parameters
+    ----------
+    solar_zenith : numeric
+        Sun's apparent zenith in degrees.
+    solar_azimuth : numeric
+        Sun's azimuth in degrees.
+    axis_tilt : numeric
+        Axis tilt angle in degrees. From horizontal plane to array plane.
+    axis_azimuth : numeric
+        Axis azimuth angle in degrees.
+        North = 0°; East = 90°; South = 180°; West = 270°
+
+    Returns
+    -------
+    Projected_solar_zenith : numeric
+        In degrees.
+
+    Notes
+    -----
+    This projection has a variety of applications in PV. For example:
+
+    - Projecting the sun's position onto the plane perpendicular to
+      the axis of a single-axis tracker (i.e. the plane
+      whose normal vector coincides with the tracker torque tube)
+      yields the tracker rotation angle that maximizes direct irradiance
+      capture. This tracking strategy is called *true-tracking*. Learn more
+      about tracking in
+      :ref:`sphx_glr_gallery_solar-tracking_plot_single_axis_tracking.py`.
+
+    - Self-shading in large PV arrays is often modeled by assuming
+      a simplified 2-D array geometry where the sun's position is
+      projected onto the plane perpendicular to the PV rows.
+      The projected zenith angle is then used for calculations
+      regarding row-to-row shading.
+
+    Examples
+    --------
+    Calculate the ideal true-tracking angle for a horizontal north-south
+    single-axis tracker:
+
+    >>> rotation = projected_solar_zenith_angle(solar_zenith, solar_azimuth,
+    >>>                                         axis_tilt=0, axis_azimuth=180)
+
+    Calculate the projected zenith angle in a south-facing fixed tilt array
+    (note: the ``axis_azimuth`` of a fixed-tilt row points along the length
+    of the row):
+
+    >>> psza = projected_solar_zenith_angle(solar_zenith, solar_azimuth,
+    >>>                                     axis_tilt=0, axis_azimuth=90)
+
+    References
+    ----------
+    .. [1] K. Anderson and M. Mikofski, 'Slope-Aware Backtracking for
+       Single-Axis Trackers', National Renewable Energy Lab. (NREL), Golden,
+       CO (United States);
+       NREL/TP-5K00-76626, Jul. 2020. :doi:`10.2172/1660126`.
+
+    See Also
+    --------
+    pvlib.solarposition.get_solarposition
+    """
+    # Assume the tracker reference frame is right-handed. Positive y-axis is
+    # oriented along tracking axis; from north, the y-axis is rotated clockwise
+    # by the axis azimuth and tilted from horizontal by the axis tilt. The
+    # positive x-axis is 90 deg clockwise from the y-axis and parallel to
+    # horizontal (e.g., if the y-axis is south, the x-axis is west); the
+    # positive z-axis is normal to the x and y axes, pointed upward.
+
+    # Since elevation = 90 - zenith, sin(90-x) = cos(x) & cos(90-x) = sin(x):
+    # Notation from [1], modified to use zenith instead of elevation
+    # cos(elevation) = sin(zenith) and sin(elevation) = cos(zenith)
+    # Avoid recalculating these values
+    sind_solar_zenith = sind(solar_zenith)
+    cosd_axis_azimuth = cosd(axis_azimuth)
+    sind_axis_azimuth = sind(axis_azimuth)
+    sind_axis_tilt = sind(axis_tilt)
+
+    # Sun's x, y, z coords
+    sx = sind_solar_zenith * sind(solar_azimuth)
+    sy = sind_solar_zenith * cosd(solar_azimuth)
+    sz = cosd(solar_zenith)
+    # Eq. (4); sx', sz' values from sun coordinates projected onto surface
+    sx_prime = sx * cosd_axis_azimuth - sy * sind_axis_azimuth
+    sz_prime = (
+        sx * sind_axis_azimuth * sind_axis_tilt
+        + sy * sind_axis_tilt * cosd_axis_azimuth
+        + sz * cosd(axis_tilt)
+    )
+    # Eq. (5); angle between sun's beam and surface
+    theta_T = np.degrees(np.arctan2(sx_prime, sz_prime))
+    return theta_T

pvlib/shading.py

@@ -2,28 +2,31 @@
 import pandas as pd
 
 from pandas.testing import assert_series_equal
+from numpy.testing import assert_allclose
 import pytest
+from datetime import timezone, timedelta
 
 from pvlib import shading
 
 
 @pytest.fixture
 def test_system():
-    syst = {'height': 1.0,
-            'pitch': 2.,
-            'surface_tilt': 30.,
-            'surface_azimuth': 180.,
-            'rotation': -30.}  # rotation of right edge relative to horizontal
-    syst['gcr'] = 1.0 / syst['pitch']
+    syst = {
+        "height": 1.0,
+        "pitch": 2.0,
+        "surface_tilt": 30.0,
+        "surface_azimuth": 180.0,
+        "rotation": -30.0,
+    }  # rotation of right edge relative to horizontal
+    syst["gcr"] = 1.0 / syst["pitch"]
     return syst
 
 
 def test__ground_angle(test_system):
     ts = test_system
-    x = np.array([0., 0.5, 1.0])
-    angles = shading.ground_angle(
-        ts['surface_tilt'], ts['gcr'], x)
-    expected_angles = np.array([0., 5.866738789543952, 9.896090638982903])
+    x = np.array([0.0, 0.5, 1.0])
+    angles = shading.ground_angle(ts["surface_tilt"], ts["gcr"], x)
+    expected_angles = np.array([0.0, 5.866738789543952, 9.896090638982903])
     assert np.allclose(angles, expected_angles)
 
 
@@ -37,7 +40,7 @@ def test__ground_angle_zero_gcr():
 
 @pytest.fixture
 def surface_tilt():
-    idx = pd.date_range('2019-01-01', freq='h', periods=3)
+    idx = pd.date_range("2019-01-01", freq="h", periods=3)
     return pd.Series([0, 20, 90], index=idx)
 
 
@@ -104,3 +107,119 @@ def test_sky_diffuse_passias_scalar(average_masking_angle, shading_loss):
     for angle, loss in zip(average_masking_angle, shading_loss):
         actual_loss = shading.sky_diffuse_passias(angle)
         assert np.isclose(loss, actual_loss)
+
+
+@pytest.fixture
+def true_tracking_angle_and_inputs_NREL():
+    # data from NREL 'Slope-Aware Backtracking for Single-Axis Trackers'
+    # doi.org/10.2172/1660126 ; Accessed on 2023-11-06.
+    tzinfo = timezone(timedelta(hours=-5))
+    axis_tilt_angle = 9.666  # deg
+    axis_azimuth_angle = 195.0  # deg
+    timedata = pd.DataFrame(
+        columns=("Apparent Elevation", "Solar Azimuth", "True-Tracking"),
+        data=(
+            (2.404287, 122.791770, -84.440),
+            (11.263058, 133.288729, -72.604),
+            (18.733558, 145.285552, -59.861),
+            (24.109076, 158.939435, -45.578),
+            (26.810735, 173.931802, -28.764),
+            (26.482495, 189.371536, -8.475),
+            (23.170447, 204.136810, 15.120),
+            (17.296785, 217.446538, 39.562),
+            (9.461862, 229.102218, 61.587),
+            (0.524817, 239.330401, 79.530),
+        ),
+    )
+    timedata.index = pd.date_range(
+        "2019-01-01T08", "2019-01-01T17", freq="1H", tz=tzinfo
+    )
+    timedata["Apparent Zenith"] = 90.0 - timedata["Apparent Elevation"]
+    return (axis_tilt_angle, axis_azimuth_angle, timedata)
+
+
+@pytest.fixture
+def projected_solar_zenith_angle_edge_cases():
+    premises_and_result_matrix = pd.DataFrame(
+        data=[
+            # s_zen | s_azm | ax_tilt | ax_azm | psza
+            [   0,       0,      0,        0,      0],
+            [   0,     180,      0,        0,      0],
+            [   0,       0,      0,      180,      0],
+            [   0,     180,      0,      180,      0],
+            [  45,       0,      0,      180,      0],
+            [  45,      90,      0,      180,    -45],
+            [  45,     270,      0,      180,     45],
+            [  45,      90,     90,      180,    -90],
+            [  45,     270,     90,      180,     90],
+            [  45,      90,     90,        0,     90],
+            [  45,     270,     90,        0,    -90],
+            [  45,      45,     90,      180,   -135],
+            [  45,     315,     90,      180,    135],
+        ],
+        columns=["solar_zenith", "solar_azimuth", "axis_tilt", "axis_azimuth",
+                 "psza"],
+    )
+    return premises_and_result_matrix
+
+
+def test_projected_solar_zenith_angle_numeric(
+    true_tracking_angle_and_inputs_NREL,
+    projected_solar_zenith_angle_edge_cases
+):
+    psza_func = shading.projected_solar_zenith_angle
+    axis_tilt, axis_azimuth, timedata = true_tracking_angle_and_inputs_NREL
+    # test against data provided by NREL
+    psz = psza_func(
+        timedata["Apparent Zenith"],
+        timedata["Solar Azimuth"],
+        axis_tilt,
+        axis_azimuth,
+    )
+    assert_allclose(psz, timedata["True-Tracking"], atol=1e-3)
+    # test by changing axis azimuth and tilt
+    psza = psza_func(
+        timedata["Apparent Zenith"],
+        timedata["Solar Azimuth"],
+        -axis_tilt,
+        axis_azimuth - 180,
+    )
+    assert_allclose(psza, -timedata["True-Tracking"], atol=1e-3)
+
+    # test edge cases
+    solar_zenith, solar_azimuth, axis_tilt, axis_azimuth, psza_expected = (
+        v for _, v in projected_solar_zenith_angle_edge_cases.items()
+    )
+    psza = psza_func(
+        solar_zenith,
+        solar_azimuth,
+        axis_tilt,
+        axis_azimuth,
+    )
+    assert_allclose(psza, psza_expected, atol=1e-9)
+
+
+@pytest.mark.parametrize(
+    "cast_type, cast_func",
+    [
+        (float, lambda x: float(x)),
+        (np.ndarray, lambda x: np.array([x])),
+        (pd.Series, lambda x: pd.Series(data=[x])),
+    ],
+)
+def test_projected_solar_zenith_angle_datatypes(
+    cast_type, cast_func, true_tracking_angle_and_inputs_NREL
+):
+    psz_func = shading.projected_solar_zenith_angle
+    axis_tilt, axis_azimuth, timedata = true_tracking_angle_and_inputs_NREL
+    sun_apparent_zenith = timedata["Apparent Zenith"].iloc[0]
+    sun_azimuth = timedata["Solar Azimuth"].iloc[0]
+
+    axis_tilt, axis_azimuth, sun_apparent_zenith, sun_azimuth = (
+        cast_func(sun_apparent_zenith),
+        cast_func(sun_azimuth),
+        cast_func(axis_tilt),
+        cast_func(axis_azimuth),
+    )
+    psz = psz_func(sun_apparent_zenith, axis_azimuth, axis_tilt, axis_azimuth)
+    assert isinstance(psz, cast_type)

pvlib/tests/test_shading.py

@@ -3,6 +3,7 @@
 
 from pvlib.tools import cosd, sind, tand, acosd, asind
 from pvlib import irradiance
+from pvlib import shading
 
 
 def singleaxis(apparent_zenith, apparent_azimuth,
@@ -126,51 +127,20 @@ def singleaxis(apparent_zenith, apparent_azimuth,
     if apparent_azimuth.ndim > 1 or apparent_zenith.ndim > 1:
         raise ValueError('Input dimensions must not exceed 1')
 
-    # Calculate sun position x, y, z using coordinate system as in [1], Eq 1.
-
-    # NOTE: solar elevation = 90 - solar zenith, then use trig identities:
-    # sin(90-x) = cos(x) & cos(90-x) = sin(x)
-    sin_zenith = sind(apparent_zenith)
-    x = sin_zenith * sind(apparent_azimuth)
-    y = sin_zenith * cosd(apparent_azimuth)
-    z = cosd(apparent_zenith)
-
-    # Assume the tracker reference frame is right-handed. Positive y-axis is
-    # oriented along tracking axis; from north, the y-axis is rotated clockwise
-    # by the axis azimuth and tilted from horizontal by the axis tilt. The
-    # positive x-axis is 90 deg clockwise from the y-axis and parallel to
-    # horizontal (e.g., if the y-axis is south, the x-axis is west); the
-    # positive z-axis is normal to the x and y axes, pointed upward.
-
-    # Calculate sun position (xp, yp, zp) in tracker coordinate system using
-    # [1] Eq 4.
-
-    cos_axis_azimuth = cosd(axis_azimuth)
-    sin_axis_azimuth = sind(axis_azimuth)
-    cos_axis_tilt = cosd(axis_tilt)
-    sin_axis_tilt = sind(axis_tilt)
-    xp = x*cos_axis_azimuth - y*sin_axis_azimuth
-    # not necessary to calculate y'
-    # yp = (x*cos_axis_tilt*sin_axis_azimuth
-    #       + y*cos_axis_tilt*cos_axis_azimuth
-    #       - z*sin_axis_tilt)
-    zp = (x*sin_axis_tilt*sin_axis_azimuth
-          + y*sin_axis_tilt*cos_axis_azimuth
-          + z*cos_axis_tilt)
-
     # The ideal tracking angle wid is the rotation to place the sun position
-    # vector (xp, yp, zp) in the (y, z) plane, which is normal to the panel and
+    # vector (xp, yp, zp) in the (x, z) plane, which is normal to the panel and
     # contains the axis of rotation.  wid = 0 indicates that the panel is
     # horizontal. Here, our convention is that a clockwise rotation is
     # positive, to view rotation angles in the same frame of reference as
     # azimuth. For example, for a system with tracking axis oriented south, a
     # rotation toward the east is negative, and a rotation to the west is
     # positive. This is a right-handed rotation around the tracker y-axis.
-
-    # Calculate angle from x-y plane to projection of sun vector onto x-z plane
-    # using [1] Eq. 5.
-
-    wid = np.degrees(np.arctan2(xp, zp))
+    wid = shading.projected_solar_zenith_angle(
+        axis_tilt=axis_tilt,
+        axis_azimuth=axis_azimuth,
+        solar_zenith=apparent_zenith,
+        solar_azimuth=apparent_azimuth,
+    )
 
     # filter for sun above panel horizon
     zen_gt_90 = apparent_zenith > 90