Fix calc_surface_orientation broadcasting for >1-D tracker_theta

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

@@ -18,6 +18,10 @@ Bug fixes
   data type integer, users can expect modeled cell temperature values to
   increase slightly.
   (:issue:`2608`, :pull:`2745`)
+* Fixes a regression in :py:func:`pvlib.tracking.calc_surface_orientation`
+  introduced in v0.15.1 (:pull:`2702`) that caused a broadcasting
+  ``ValueError`` when ``tracker_theta`` was a 2-D (or higher rank) array.
+  (:issue:`2747`, :pull:`2749`)
 
 Enhancements
 ~~~~~~~~~~~~
@@ -53,3 +57,4 @@ Contributors
 ~~~~~~~~~~~~
 * :ghuser:`Omesh37`
 * Cliff Hansen (:ghuser:`cwhanse`)
+* Arthur Onno (:ghuser:`ArthurOnnoTerabase`)

pvlib/tracking.py

@@ -233,10 +233,11 @@ def _unit_normal(axis_azimuth, axis_tilt, theta):
     Returns
     -------
     ndarray
-        Shape (N,3) where theta has length N
+        Shape ``theta.shape + (3,)``, with a minimum rank of 2.  For 1-D
+        ``theta`` of length N this is ``(N, 3)``.
     """
 
-    theta = np.asarray(theta)
+    theta = np.atleast_1d(np.asarray(theta))
 
     cA, sA = cosd(-axis_azimuth), sind(-axis_azimuth)
     cT, sT = cosd(-axis_tilt), sind(-axis_tilt)
@@ -248,7 +249,7 @@ def _unit_normal(axis_azimuth, axis_tilt, theta):
     y = sA * sTh - cA * sT * cTh
     z = cT * cTh
 
-    result = np.column_stack((x, y, z))
+    result = np.stack((x, y, z), axis=-1)
 
     return result
 
@@ -296,7 +297,7 @@ def calc_surface_orientation(tracker_theta, axis_tilt=0, axis_azimuth=0):
 
     # project unit_normal to x-y plane to calculate azimuth
     surface_azimuth = np.degrees(
-        np.arctan2(unit_normal[:, 0], unit_normal[:, 1]))
+        np.arctan2(unit_normal[..., 0], unit_normal[..., 1]))
 
     surface_azimuth = np.where(surface_tilt == 0., axis_azimuth - 90.,
                                surface_azimuth)

tests/test_tracking.py

@@ -557,3 +557,32 @@ def test_calc_surface_orientation_special():
         # in a modulo-360 sense.
         np.testing.assert_allclose(np.round(out['surface_azimuth'], 4) % 360,
                                    expected_azimuths, rtol=1e-5, atol=1e-5)
+
+
+@pytest.mark.parametrize('shape', [(3, 5), (1, 7), (4, 1), (2, 3, 4)])
+def test_calc_surface_orientation_2d(shape):
+    # Regression test for GH#2747: calc_surface_orientation must accept
+    # tracker_theta of arbitrary rank, not just 1-D.  Compare the >1-D result
+    # to the 1-D result computed on the flattened input.
+    rng = np.random.default_rng(0)
+    rotations_flat = rng.uniform(-90, 90, size=int(np.prod(shape)))
+    rotations_nd = rotations_flat.reshape(shape)
+
+    out_1d = tracking.calc_surface_orientation(rotations_flat,
+                                               axis_tilt=20,
+                                               axis_azimuth=180)
+    out_nd = tracking.calc_surface_orientation(rotations_nd,
+                                               axis_tilt=20,
+                                               axis_azimuth=180)
+
+    assert out_nd['surface_tilt'].shape == shape
+    assert out_nd['surface_azimuth'].shape == shape
+    np.testing.assert_allclose(out_nd['surface_tilt'].reshape(-1),
+                               out_1d['surface_tilt'])
+    np.testing.assert_allclose(out_nd['surface_azimuth'].reshape(-1),
+                               out_1d['surface_azimuth'])
+
+    # _unit_normal must preserve the input rank, appending a trailing axis
+    # of length 3.
+    unorm = tracking._unit_normal(180., 20., rotations_nd)
+    assert unorm.shape == shape + (3,)