calc_surface_orientation broadcasting error for 2-D inputs in 0.15.1 (regression vs 0.15.0)

[ArthurOnnoTerabase]

Describe the bug

pvlib.tracking.calc_surface_orientation regressed in 0.15.1: passing 2-D arrays for tracker_theta (and broadcasting-compatible 2-D axis_tilt / axis_azimuth, e.g. shape (timestamps, sites)) raises a broadcasting ValueError. The same call works in 0.15.0 and earlier.

The cause is the new _unit_normal helper introduced in #2702, which builds its return value with np.column_stack((x, y, z)). np.column_stack only stacks 1-D inputs into columns; for 2-D inputs of shape (T, N) it concatenates along axis=1, producing a (T, 3*N) array instead of the intended (T, N, 3). The subsequent unit_normal[:, 0] / unit_normal[:, 1] then collapse to shape (T,), and the downstream

python surface_azimuth = np.where(surface_tilt == 0., axis_azimuth - 90., surface_azimuth)

fails because the condition and true branch are still (T, N) while the false branch is now (T,).

To Reproduce

import numpy as np
import pvlib

T, N = 4655, 105
tracker_theta = np.zeros((T, N))
axis_tilt = np.zeros((T, N))
axis_azimuth = np.full((T, N), 180.0)

pvlib.tracking.calc_surface_orientation(tracker_theta, axis_tilt, axis_azimuth)
# pvlib 0.15.0: returns dict of (T, N) arrays — works
# pvlib 0.15.1: ValueError: operands could not be broadcast together with shapes (4655,105) (4655,105) (4655,)

Bisected against pip-installed releases on Python 3.11 / numpy 2.4.4 / pandas 3.0.2:

pvlib	result
0.11.2 → 0.15.0	works
0.15.1	broken

Expected behavior

Return surface_tilt and surface_azimuth arrays with the same shape as the broadcast inputs (the 0.15.0 behavior).

Versions

• pvlib.__version__: 0.15.1 (broken), 0.15.0 (works)
• pandas.__version__: 3.0.2
• numpy.__version__: 2.4.4
• python: 3.11.9
• platform: Windows

Additional context

Regression introduced in #2702 (milestone v0.15.1). Suggested one-line fix in
_unit_normal:

# before
return np.column_stack((x, y, z))
# after
return np.stack((x, y, z), axis=-1)

np.stack(..., axis=-1) preserves the leading dimensions for inputs of any rank and is equivalent to the old behavior for 1-D inputs, so it should be a drop-in replacement.

Hit in production while batch-processing tracker geometry across (timestamps × bays); pinned to pvlib<0.15.1 as a workaround.

[ArthurOnnoTerabase]

@cwhanse / @kandersolar — this looks like fallout from #2702 and would be a nice candidate for the v0.15.2 milestone. I can send a PR with the np.stack(..., axis=-1) fix and a 2-D regression test if it would help.

[cwhanse]

@ArthurOnnoTerabase absolutely! Reduction to 1d (time) was not intended. Although the previous behavior (T, N) was not tested nor intentional, I'm happy to hear that it worked :)

We'd welcome a PR and additions to the test suite.

[ArthurOnnoTerabase]

Submitted #2749 for this.

The fix turned out to be slightly larger than the one-line column_stack → stack(..., axis=-1) I suggested above. Two additional changes were needed:

• unit_normal[:, 0]/[:, 1] → unit_normal[..., 0]/[..., 1] in calc_surface_orientation (ellipsis indexing so it works for arbitrary input rank).
• np.atleast_1d(np.asarray(theta)) in _unit_normal to preserve the historical (1, 3) return shape for scalar input — without it the existing test__unit_normal regresses on the scalar case.

All 22 existing test_tracking.py tests still pass. Added a new parametrised regression test (test_calc_surface_orientation_2d) covering shapes (3, 5), (1, 7), (4, 1), and (2, 3, 4), which fails on main with the exact ValueError from the original report.