diff --git a/docs/sphinx/source/reference/tracking.rst b/docs/sphinx/source/reference/tracking.rst index 7f2599b6a0..f4bbcc9aa9 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 4219f9c3cb..d1c5290738 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 893cc65639..d162c68719 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 ~~~~~~~~~ diff --git a/pvlib/tests/test_tracking.py b/pvlib/tests/test_tracking.py index c88c92b248..0fbace0f17 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) diff --git a/pvlib/tools.py b/pvlib/tools.py index 5c6e1dd293..61e6d170c4 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 951f2e886e..d9cd2b7853 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