Move surface orientation calculation from tracking.singleaxis to new function (#1480)

* pvlib.tools.acosd

* add calc_surface_orientation

* prune singleaxis(), call calc_surface_orientation

* tests

* whatsnew

* add calc_surface_orientation to sphinx api list

* pr number

* a bit of docstring cleanup

* reference doi

* include 0.9.2 on whatsnew page

* Apply suggestions from code review

Co-authored-by: Cliff Hansen <[email protected]>

* mention right-handed rotation requirement

* fix sphinx rendering issue

* Apply suggestions from code review

Co-authored-by: Cliff Hansen <[email protected]>

Co-authored-by: Cliff Hansen <[email protected]>

Author: kandersolar

docs/sphinx/source/reference/tracking.rst

@@ -25,3 +25,4 @@ Functions
    tracking.singleaxis
    tracking.calc_axis_tilt
    tracking.calc_cross_axis_tilt
+   tracking.calc_surface_orientation

docs/sphinx/source/whatsnew.rst

@@ -6,6 +6,7 @@ What's New
 
 These are new features and improvements of note in each release.
 
+.. include:: whatsnew/v0.9.2.rst
 .. include:: whatsnew/v0.9.1.rst
 .. include:: whatsnew/v0.9.0.rst
 .. include:: whatsnew/v0.8.1.rst

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

@@ -8,6 +8,9 @@ Deprecations
 
 Enhancements
 ~~~~~~~~~~~~
+* Add :py:func:`pvlib.tracking.calc_surface_orientation` for calculating
+  single-axis tracker ``surface_tilt`` and ``surface_azimuth`` from
+  rotation angles. (:issue:`1471`, :pull:`1480`)
 
 Bug fixes
 ~~~~~~~~~

pvlib/tests/test_tracking.py

@@ -517,3 +517,72 @@ def test_singleaxis_aoi_gh1221():
     fixed = pvlib.irradiance.aoi(90, 180, sp['apparent_zenith'], sp['azimuth'])
     fixed[np.isnan(tr['aoi'])] = np.nan
     assert np.allclose(tr['aoi'], fixed, equal_nan=True)
+
+
+def test_calc_surface_orientation_types():
+    # numpy arrays
+    rotations = np.array([-10, 0, 10])
+    expected_tilts = np.array([10, 0, 10], dtype=float)
+    expected_azimuths = np.array([270, 90, 90], dtype=float)
+    out = tracking.calc_surface_orientation(tracker_theta=rotations)
+    np.testing.assert_allclose(expected_tilts, out['surface_tilt'])
+    np.testing.assert_allclose(expected_azimuths, out['surface_azimuth'])
+
+    # pandas Series
+    rotations = pd.Series(rotations)
+    expected_tilts = pd.Series(expected_tilts).rename('surface_tilt')
+    expected_azimuths = pd.Series(expected_azimuths).rename('surface_azimuth')
+    out = tracking.calc_surface_orientation(tracker_theta=rotations)
+    assert_series_equal(expected_tilts, out['surface_tilt'])
+    assert_series_equal(expected_azimuths, out['surface_azimuth'])
+
+    # float
+    for rotation, expected_tilt, expected_azimuth in zip(
+            rotations, expected_tilts, expected_azimuths):
+        out = tracking.calc_surface_orientation(rotation)
+        assert out['surface_tilt'] == pytest.approx(expected_tilt)
+        assert out['surface_azimuth'] == pytest.approx(expected_azimuth)
+
+
+def test_calc_surface_orientation_kwargs():
+    # non-default axis tilt & azimuth
+    rotations = np.array([-10, 0, 10])
+    expected_tilts = np.array([22.2687445, 20.0, 22.2687445])
+    expected_azimuths = np.array([152.72683041, 180.0, 207.27316959])
+    out = tracking.calc_surface_orientation(rotations,
+                                            axis_tilt=20,
+                                            axis_azimuth=180)
+    np.testing.assert_allclose(out['surface_tilt'], expected_tilts)
+    np.testing.assert_allclose(out['surface_azimuth'], expected_azimuths)
+
+
+def test_calc_surface_orientation_special():
+    # special cases for rotations
+    rotations = np.array([-180, -90, -0, 0, 90, 180])
+    expected_tilts = np.array([180, 90, 0, 0, 90, 180], dtype=float)
+    expected_azimuths = [270, 270, 90, 90, 90, 90]
+    out = tracking.calc_surface_orientation(rotations)
+    np.testing.assert_allclose(out['surface_tilt'], expected_tilts)
+    np.testing.assert_allclose(out['surface_azimuth'], expected_azimuths)
+
+    # special case for axis_tilt
+    rotations = np.array([-10, 0, 10])
+    expected_tilts = np.array([90, 90, 90], dtype=float)
+    expected_azimuths = np.array([350, 0, 10], dtype=float)
+    out = tracking.calc_surface_orientation(rotations, axis_tilt=90)
+    np.testing.assert_allclose(out['surface_tilt'], expected_tilts)
+    np.testing.assert_allclose(out['surface_azimuth'], expected_azimuths)
+
+    # special cases for axis_azimuth
+    rotations = np.array([-10, 0, 10])
+    expected_tilts = np.array([10, 0, 10], dtype=float)
+    expected_azimuth_offsets = np.array([-90, 90, 90], dtype=float)
+    for axis_azimuth in [0, 90, 180, 270, 360]:
+        expected_azimuths = (axis_azimuth + expected_azimuth_offsets) % 360
+        out = tracking.calc_surface_orientation(rotations,
+                                                axis_azimuth=axis_azimuth)
+        np.testing.assert_allclose(out['surface_tilt'], expected_tilts)
+        # the rounding is a bit ugly, but necessary to test approximately equal
+        # in a modulo-360 sense.
+        np.testing.assert_allclose(np.round(out['surface_azimuth'], 4) % 360,
+                                   expected_azimuths, rtol=1e-5, atol=1e-5)

pvlib/tools.py

@@ -85,6 +85,25 @@ def asind(number):
     return res
 
 
+def acosd(number):
+    """
+    Inverse Cosine returning an angle in degrees
+
+    Parameters
+    ----------
+    number : float
+        Input number
+
+    Returns
+    -------
+    result : float
+        arccos result
+    """
+
+    res = np.degrees(np.arccos(number))
+    return res
+
+
 def localize_to_utc(time, location):
     """
     Converts or localizes a time series to UTC.

pvlib/tracking.py

@@ -1,7 +1,7 @@
 import numpy as np
 import pandas as pd
 
-from pvlib.tools import cosd, sind, tand
+from pvlib.tools import cosd, sind, tand, acosd, asind
 from pvlib.pvsystem import (
     PVSystem, Array, SingleAxisTrackerMount, _unwrap_single_value
 )
@@ -334,9 +334,9 @@ def singleaxis(apparent_zenith, apparent_azimuth,
     Returns
     -------
     dict or DataFrame with the following columns:
-        * `tracker_theta`: The rotation angle of the tracker.
-          tracker_theta = 0 is horizontal, and positive rotation angles are
-          clockwise. [degrees]
+        * `tracker_theta`: The rotation angle of the tracker is a right-handed
+          rotation defined by `axis_azimuth`.
+          tracker_theta = 0 is horizontal. [degrees]
         * `aoi`: The angle-of-incidence of direct irradiance onto the
           rotated panel surface. [degrees]
         * `surface_tilt`: The angle between the panel surface and the earth
@@ -349,6 +349,7 @@ def singleaxis(apparent_zenith, apparent_azimuth,
     --------
     pvlib.tracking.calc_axis_tilt
     pvlib.tracking.calc_cross_axis_tilt
+    pvlib.tracking.calc_surface_orientation
 
     References
     ----------
@@ -396,9 +397,10 @@ def singleaxis(apparent_zenith, apparent_azimuth,
     cos_axis_tilt = cosd(axis_tilt)
     sin_axis_tilt = sind(axis_tilt)
     xp = x*cos_axis_azimuth - y*sin_axis_azimuth
-    yp = (x*cos_axis_tilt*sin_axis_azimuth
-          + y*cos_axis_tilt*cos_axis_azimuth
-          - z*sin_axis_tilt)
+    # 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)
@@ -446,88 +448,79 @@ def singleaxis(apparent_zenith, apparent_azimuth,
     # system-plane normal
     tracker_theta = np.clip(tracker_theta, -max_angle, max_angle)
 
-    # Calculate panel normal vector in panel-oriented x, y, z coordinates.
-    # y-axis is axis of tracker rotation. tracker_theta is a compass angle
-    # (clockwise is positive) rather than a trigonometric angle.
-    # NOTE: the *0 is a trick to preserve NaN values.
-    panel_norm = np.array([sind(tracker_theta),
-                           tracker_theta*0,
-                           cosd(tracker_theta)])
-
-    # sun position in vector format in panel-oriented x, y, z coordinates
-    sun_vec = np.array([xp, yp, zp])
-
-    # calculate angle-of-incidence on panel
-    # TODO: use irradiance.aoi
-    projection = np.clip(np.sum(sun_vec*panel_norm, axis=0), -1, 1)
-    aoi = np.degrees(np.arccos(projection))
-
-    # Calculate panel tilt and azimuth in a coordinate system where the panel
-    # tilt is the angle from horizontal, and the panel azimuth is the compass
-    # angle (clockwise from north) to the projection of the panel's normal to
-    # the earth's surface. These outputs are provided for convenience and
-    # comparison with other PV software which use these angle conventions.
-
-    # Project normal vector to earth surface. First rotate about x-axis by
-    # angle -axis_tilt so that y-axis is also parallel to earth surface, then
-    # project.
-
-    # Calculate standard rotation matrix
-    rot_x = np.array([[1, 0, 0],
-                      [0, cosd(-axis_tilt), -sind(-axis_tilt)],
-                      [0, sind(-axis_tilt), cosd(-axis_tilt)]])
-
-    # panel_norm_earth contains the normal vector expressed in earth-surface
-    # coordinates (z normal to surface, y aligned with tracker axis parallel to
-    # earth)
-    panel_norm_earth = np.dot(rot_x, panel_norm).T
-
-    # projection to plane tangent to earth surface, in earth surface
-    # coordinates
-    projected_normal = np.array([panel_norm_earth[:, 0],
-                                 panel_norm_earth[:, 1],
-                                 panel_norm_earth[:, 2]*0]).T
-
-    # calculate vector magnitudes
-    projected_normal_mag = np.sqrt(np.nansum(projected_normal**2, axis=1))
-
-    # renormalize the projected vector, avoid creating nan values.
-    non_zeros = projected_normal_mag != 0
-    projected_normal[non_zeros] = (projected_normal[non_zeros].T /
-                                   projected_normal_mag[non_zeros]).T
-
-    # calculation of surface_azimuth
-    surface_azimuth = \
-        np.degrees(np.arctan2(projected_normal[:, 1], projected_normal[:, 0]))
-
-    # Rotate 0 reference from panel's x-axis to its y-axis and then back to
-    # north.
-    surface_azimuth = 90 - surface_azimuth + axis_azimuth
-
-    # Map azimuth into [0,360) domain.
-    with np.errstate(invalid='ignore'):
-        surface_azimuth = surface_azimuth % 360
-
-    # Calculate surface_tilt
-    dotproduct = (panel_norm_earth * projected_normal).sum(axis=1)
-    # for edge cases like axis_tilt=90, numpy's SIMD can produce values like
-    # dotproduct = (1 + 2e-16). Clip off the excess so that arccos works:
-    dotproduct = np.clip(dotproduct, -1, 1)
-    surface_tilt = 90 - np.degrees(np.arccos(dotproduct))
+    # Calculate auxiliary angles
+    surface = calc_surface_orientation(tracker_theta, axis_tilt, axis_azimuth)
+    surface_tilt = surface['surface_tilt']
+    surface_azimuth = surface['surface_azimuth']
+    aoi = irradiance.aoi(surface_tilt, surface_azimuth,
+                         apparent_zenith, apparent_azimuth)
 
     # Bundle DataFrame for return values and filter for sun below horizon.
     out = {'tracker_theta': tracker_theta, 'aoi': aoi,
-           'surface_tilt': surface_tilt, 'surface_azimuth': surface_azimuth}
+           'surface_azimuth': surface_azimuth, 'surface_tilt': surface_tilt}
     if index is not None:
         out = pd.DataFrame(out, index=index)
-        out = out[['tracker_theta', 'aoi', 'surface_azimuth', 'surface_tilt']]
         out[zen_gt_90] = np.nan
     else:
         out = {k: np.where(zen_gt_90, np.nan, v) for k, v in out.items()}
 
     return out
 
 
+def calc_surface_orientation(tracker_theta, axis_tilt=0, axis_azimuth=0):
+    """
+    Calculate the surface tilt and azimuth angles for a given tracker rotation.
+
+    Parameters
+    ----------
+    tracker_theta : numeric
+        Tracker rotation angle as a right-handed rotation around
+        the axis defined by ``axis_tilt`` and ``axis_azimuth``.  For example,
+        with ``axis_tilt=0`` and ``axis_azimuth=180``, ``tracker_theta > 0``
+        results in ``surface_azimuth`` to the West while ``tracker_theta < 0``
+        results in ``surface_azimuth`` to the East. [degree]
+    axis_tilt : float, default 0
+        The tilt of the axis of rotation with respect to horizontal. [degree]
+    axis_azimuth : float, default 0
+        A value denoting the compass direction along which the axis of
+        rotation lies. Measured east of north. [degree]
+
+    Returns
+    -------
+    dict or DataFrame
+        Contains keys ``'surface_tilt'`` and ``'surface_azimuth'`` representing
+        the module orientation accounting for tracker rotation and axis
+        orientation. [degree]
+
+    References
+    ----------
+    .. [1] William F. Marion and Aron P. Dobos, "Rotation Angle for the Optimum
+       Tracking of One-Axis Trackers", Technical Report NREL/TP-6A20-58891,
+       July 2013. :doi:`10.2172/1089596`
+    """
+    with np.errstate(invalid='ignore', divide='ignore'):
+        surface_tilt = acosd(cosd(tracker_theta) * cosd(axis_tilt))
+
+        # clip(..., -1, +1) to prevent arcsin(1 + epsilon) issues:
+        azimuth_delta = asind(np.clip(sind(tracker_theta) / sind(surface_tilt),
+                                      a_min=-1, a_max=1))
+        # Combine Eqs 2, 3, and 4:
+        azimuth_delta = np.where(abs(tracker_theta) < 90,
+                                 azimuth_delta,
+                                 -azimuth_delta + np.sign(tracker_theta) * 180)
+        # handle surface_tilt=0 case:
+        azimuth_delta = np.where(sind(surface_tilt) != 0, azimuth_delta, 90)
+        surface_azimuth = (axis_azimuth + azimuth_delta) % 360
+
+    out = {
+        'surface_tilt': surface_tilt,
+        'surface_azimuth': surface_azimuth,
+    }
+    if hasattr(tracker_theta, 'index'):
+        out = pd.DataFrame(out)
+    return out
+
+
 def calc_axis_tilt(slope_azimuth, slope_tilt, axis_azimuth):
     """
     Calculate tracker axis tilt in the global reference frame when on a sloped