Repository: pvlib/pvlib-python
Pull request #1470 (2022-06-09 20:54:43)

Gallery example for modeling fixed-tilt arrays with pvfactors
 - ~[ ] Closes #xxxx~
 - [x] I am familiar with the [contributing guidelines](https://pvlib-python.readthedocs.io/en/latest/contributing.html)
 - ~[ ] Tests added~
 - ~[ ] Updates entries in [`docs/sphinx/source/reference`](https://github.com/pvlib/pvlib-python/blob/master/docs/sphinx/source/reference) for API changes.~
 - [x] Adds description and name entries in the appropriate "what's new" file in [`docs/sphinx/source/whatsnew`](https://github.com/pvlib/pvlib-python/tree/master/docs/sphinx/source/whatsnew) for all changes. Includes link to the GitHub Issue with `` :issue:`num` `` or this Pull Request with `` :pull:`num` ``. Includes contributor name and/or GitHub username (link with `` :ghuser:`user` ``).
 - [x] New code is fully documented. Includes [numpydoc](https://numpydoc.readthedocs.io/en/latest/format.html) compliant docstrings, examples, and comments where necessary.
 - [x] Pull request is nearly complete and ready for detailed review.
 - [x] Maintainer: Appropriate GitHub Labels and Milestone are assigned to the Pull Request and linked Issue.

Motivated by two interactions on the google group where the role of the `axis_azimuth` input for the pvfactors model was not clear when modeling fixed-tilt arrays:
- https://groups.google.com/g/pvlib-python/c/rPFn9G94Qf4/m/a1PQX5QPAgAJ
- https://groups.google.com/g/pvlib-python/c/0MFG0J59Xnc/m/ZqZyWmeJCQAJ

---

Rendered example here: https://pvlib-python--1470.org.readthedocs.build/en/1470/gallery/bifacial/plot_pvfactors_fixed_tilt.html


Patch:
diff --git a/.landscape.yml b/.landscape.yml
deleted file mode 100644
index 2a56757..0000000
--- a/.landscape.yml
+++ /dev/null
@@ -1,11 +0,0 @@
-doc-warnings: yes
-test-warnings: yes
-strictness: medium
-max-line-length: 79
-autodetect: yes
-ignore-paths:
-    - ci
-    - docs
-pep257:
-    disable:
-        - D202
diff --git a/README.md b/README.md
index 211eb7f..88c1326 100644
--- a/README.md
+++ b/README.md
@@ -28,8 +28,11 @@
     <a href="http://pvlib-python.readthedocs.org/en/stable/">
     <img src="https://readthedocs.org/projects/pvlib-python/badge/?version=stable" alt="documentation build status" />
     </a>
-    <a href="https://dev.azure.com/solararbiter/pvlib%20python/_build/latest?definitionId=4&branchName=master">
-      <img src="https://dev.azure.com/solararbiter/pvlib%20python/_apis/build/status/pvlib.pvlib-python?branchName=master" alt="Azure Pipelines build status" />
+    <a href="https://github.com/pvlib/pvlib-python/actions/workflows/pytest.yml?query=branch%3Amaster">
+      <img src="https://github.com/pvlib/pvlib-python/actions/workflows/pytest.yml/badge.svg?branch=master" alt="GitHub Actions Testing Status" />
+    </a>
+    <a href="https://codecov.io/gh/pvlib/pvlib-python">
+    <img src="https://codecov.io/gh/pvlib/pvlib-python/branch/master/graph/badge.svg" alt="codecov coverage" />
     </a>
   </td>
 </tr>
@@ -44,17 +47,6 @@
     </a>
   </td>
 </tr>
-<tr>
-  <td>Coverage</td>
-  <td>
-    <a href="https://coveralls.io/r/pvlib/pvlib-python">
-    <img src="https://img.shields.io/coveralls/pvlib/pvlib-python.svg" alt="coveralls coverage" />
-    </a>
-    <a href="https://codecov.io/gh/pvlib/pvlib-python">
-    <img src="https://codecov.io/gh/pvlib/pvlib-python/branch/master/graph/badge.svg" alt="codecov coverage" />
-    </a>
-  </td>
-</tr>
 <tr>
   <td>Benchmarks</td>
   <td>
diff --git a/docs/examples/bifacial/plot_pvfactors_fixed_tilt.py b/docs/examples/bifacial/plot_pvfactors_fixed_tilt.py
new file mode 100644
index 0000000..5525c49
--- /dev/null
+++ b/docs/examples/bifacial/plot_pvfactors_fixed_tilt.py
@@ -0,0 +1,68 @@
+"""
+Fixed-Tilt Simulation with pvfactors
+====================================
+
+Modeling the irradiance on the rear side of a fixed-tilt array.
+"""
+
+# %%
+# Because pvfactors was originally designed for modeling single-axis
+# tracking systems, it's not necessarily obvious how to use it to model
+# fixed-tilt systems correctly.
+# This example shows how to model rear-side irradiance on a fixed-tilt
+# array using :py:func:`pvlib.bifacial.pvfactors.pvfactors_timeseries`.
+
+import pandas as pd
+from pvlib import location
+from pvlib.bifacial.pvfactors import pvfactors_timeseries
+import matplotlib.pyplot as plt
+import warnings
+
+# supressing shapely warnings that occur on import of pvfactors
+warnings.filterwarnings(action='ignore', module='pvfactors')
+
+# %%
+# First, generate the usual modeling inputs:
+
+times = pd.date_range('2021-06-21', '2021-06-22', freq='1T', tz='Etc/GMT+5')
+loc = location.Location(latitude=40, longitude=-80, tz=times.tz)
+sp = loc.get_solarposition(times)
+cs = loc.get_clearsky(times)
+
+# example array geometry
+pvrow_height = 1
+pvrow_width = 4
+pitch = 10
+gcr = pvrow_width / pitch
+axis_azimuth = 180
+albedo = 0.2
+
+# %%
+# Now the trick: since pvfactors only wants to model single-axis tracking
+# arrays, we have to pretend our fixed tilt array is a single-axis tracking
+# array that never rotates.  In that case, the "axis of rotation" is
+# along the length of the row, with ``axis_azimuth`` 90 degrees offset from the
+# fixed ``surface_azimuth``.
+
+irrad = pvfactors_timeseries(
+    solar_azimuth=sp['azimuth'],
+    solar_zenith=sp['apparent_zenith'],
+    surface_azimuth=180,  # south-facing array
+    surface_tilt=20,
+    axis_azimuth=90,  # 90 degrees off from surface_azimuth.  270 is ok too
+    timestamps=times,
+    dni=cs['dni'],
+    dhi=cs['dhi'],
+    gcr=gcr,
+    pvrow_height=pvrow_height,
+    pvrow_width=pvrow_width,
+    albedo=albedo,
+    n_pvrows=3,
+    index_observed_pvrow=1
+)
+
+# turn into pandas DataFrame
+irrad = pd.concat(irrad, axis=1)
+
+irrad[['total_inc_back', 'total_abs_back']].plot()
+plt.ylabel('Irradiance [W m$^{-2}$]')
diff --git a/docs/sphinx/source/reference/tracking.rst b/docs/sphinx/source/reference/tracking.rst
index 7f2599b..f4bbcc9 100644
--- a/docs/sphinx/source/reference/tracking.rst
+++ b/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
diff --git a/docs/sphinx/source/whatsnew.rst b/docs/sphinx/source/whatsnew.rst
index 4219f9c..d1c5290 100644
--- a/docs/sphinx/source/whatsnew.rst
+++ b/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
diff --git a/docs/sphinx/source/whatsnew/v0.9.2.rst b/docs/sphinx/source/whatsnew/v0.9.2.rst
index 893cc65..0d1e14b 100644
--- a/docs/sphinx/source/whatsnew/v0.9.2.rst
+++ b/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
 ~~~~~~~~~
@@ -27,6 +30,9 @@ Testing
 
 Documentation
 ~~~~~~~~~~~~~
+* Add gallery example of simulating rearside irradiance for a fixed-tilt
+  array with pvfactors (:pull:`1470`)
+
 
 Benchmarking
 ~~~~~~~~~~~~~
diff --git a/pvlib/bifacial/pvfactors.py b/pvlib/bifacial/pvfactors.py
index 3c70b77..ed09879 100644
--- a/pvlib/bifacial/pvfactors.py
+++ b/pvlib/bifacial/pvfactors.py
@@ -35,6 +35,8 @@ def pvfactors_timeseries(
     axis_azimuth: float
         Azimuth angle of the rotation axis of the PV modules, using pvlib's
         convention (deg). This is supposed to be fixed for all timestamps.
+        When modeling fixed-tilt arrays, set this value to be 90 degrees
+        clockwise from ``surface_azimuth``.
     timestamps: datetime or DatetimeIndex
         List of simulation timestamps
     dni: numeric
diff --git a/pvlib/tools.py b/pvlib/tools.py
index 5c6e1dd..61e6d17 100644
--- a/pvlib/tools.py
+++ b/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.
diff --git a/pvlib/tracking.py b/pvlib/tracking.py
index 951f2e8..d9cd2b7 100644
--- a/pvlib/tracking.py
+++ b/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,81 +448,18 @@ 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()}
@@ -528,6 +467,60 @@ def singleaxis(apparent_zenith, apparent_azimuth,
     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


Test patch:
diff --git a/pvlib/tests/test_tracking.py b/pvlib/tests/test_tracking.py
index c88c92b..0fbace0 100644
--- a/pvlib/tests/test_tracking.py
+++ b/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)