Consider extracting the surface orientation calculation in pvlib.tracking.singleaxis() to its own function #1471

kandersolar · 2022-06-14T01:53:45Z

Is your feature request related to a problem? Please describe.
The usual workflow for modeling single-axis tracking in pvlib is to treat tracker rotation (tracker_theta) as an unknown to be calculated from solar position and array geometry. However, sometimes a user might have their own tracker rotations but not have the corresponding surface_tilt and surface_azimuth values. Here are a few motivating examples:

Using measured rotation angles
Post-processing the output of tracking.singleaxis() to include wind stow events or tracker stalls
Other tracking algorithms that determine rotation differently from the astronomical method

Assuming I have my tracker rotations already in hand, getting the corresponding surface_tilt and surface_azimuth angles is not as easy as it should be. For the specific case of horizontal N-S axis the math isn't so bad, but either way it's annoying to have to DIY when pvlib already has code to calculate those angles from tracker rotation.

Describe the solution you'd like
A function pvlib.tracking.rotation_to_orientation that implements the same math in pvlib.tracking.singleaxis to go from tracker_theta to surface_tilt and surface_azimuth. Basically extract out the second half of tracking.singleaxis into a new function. Suggestions for the function name are welcome. To be explicit, this is more or less what I'm imagining:

def rotation_to_orientation(tracker_theta, axis_tilt=0, axis_azimuth=0, max_angle=90):
    # insert math from second half of tracking.singleaxis() here
    out = {'tracker_theta': tracker_theta, 'aoi': aoi,
           'surface_tilt': surface_tilt, 'surface_azimuth': surface_azimuth}
    return pandas_if_needed(out)

Describe alternatives you've considered
Continue suffering

Additional context
This is one step towards a broader goal I have for pvlib.tracking to house other methods to determine tracker rotation in addition to the current astronomical method, the same way we have multiple temperature and transposition models. These functions would be responsible for determining tracker rotations, and they'd all use this rotation_to_orientation function to convert rotation to module orientation.

Separately, I wonder if the code could be simplified using the tilt and azimuth equations in Bill's technical report (https://www.nrel.gov/docs/fy13osti/58891.pdf) -- seems like what we're doing is overly complicated, although maybe I've just not studied it closely enough.

cc @williamhobbs @spaneja

The text was updated successfully, but these errors were encountered:

williamhobbs · 2022-06-14T21:08:31Z

I like this.

This is related to an issue submitted for NREL SAM, NREL/SAM#850, and I think @mjprilliman is looking at something related.

kurt-rhee · 2022-06-16T23:05:18Z

@kanderso-nrel Nice meeting you at PVSC the other day. I've been working on this a bit using the tilt and azimuth equations in the technical report by Bill Marion mentioned. I would like to use this for the custom backtracking schedules Nevados generates for the terrain following trackers.

I have surface tilt working. It requires a dataframe df of axis tilt angles. I use the dataframe because Nevados has a different axis tilt for each bay in our tracker. stdf (surface tilt dataframe) starts as a dataframe of rotation angles (theta) indexed by timestep and is transformed into surface tilts. col.name is used to match the tracker and bay's rotation angle to their corresponding axis tilt.

        def calc_surface_tilt(col):
            axis_tilt = df.at[col.name, 'axis_tilt']
            surface_tilt = np.rad2deg(
                np.arccos(
                    np.cos(np.deg2rad(col)) * np.cos(np.deg2rad(axis_tilt))
                )
            )
            return surface_tilt

        stdf = stdf.apply(calc_surface_tilt, axis=0)

Unfortunately I can't seem to get surface azimuth working correctly. sadf (surface angle data frame) is almost equal to surface azimuth as calculated by pvlib, but not quite in the middle of the day. ts2 is the output of pvlib.tracking.singleaxis

    sadf = np.rad2deg(
        np.deg2rad(180) +
        np.arcsin(
            (
                np.sin(np.deg2rad(ts2['tracker_theta'])) /
                np.sin(np.deg2rad(ts2['surface_tilt']))
            ).clip(upper=1, lower=-1)
        )
    )

kandersolar · 2022-06-17T12:34:40Z

Thanks @kurt-rhee for this investigation. Trying some simple examples on my end, things seem to line up. Here's a complete copy/pasteable example where I get negligible difference between the current pvlib approach and your code. Note that I did replace the hard-coded axis_azimuth of 180 in the surface_azimuth calculation.

Click to expand!

import pvlib
import pandas as pd
import numpy as np

axis_tilt = 20
axis_azimuth = 230
loc = pvlib.location.Location(40, -80)
times = pd.date_range('2019-06-01', '2019-06-02', freq='5T', tz='Etc/GMT+5')
sp = loc.get_solarposition(times)
tr = pvlib.tracking.singleaxis(sp.apparent_zenith, sp.azimuth,
                               axis_tilt=axis_tilt, axis_azimuth=axis_azimuth)


def rotation_to_orientation(tracker_theta, axis_tilt=0, axis_azimuth=0, max_angle=90):
    surface_tilt = np.rad2deg(
        np.arccos(
            np.cos(np.deg2rad(tracker_theta)) * np.cos(np.deg2rad(axis_tilt))
        )
    )
    surface_azimuth = np.rad2deg(
        np.deg2rad(axis_azimuth) +
        np.arcsin(
            (
                np.sin(np.deg2rad(tracker_theta)) /
                np.sin(np.deg2rad(surface_tilt))
            ).clip(upper=1, lower=-1)
        )
    )
    return pd.DataFrame({
        'tracker_theta': tracker_theta,
        'surface_tilt': surface_tilt,
        'surface_azimuth': surface_azimuth,
    })

tr2 = rotation_to_orientation(tr.tracker_theta, axis_tilt=axis_tilt, axis_azimuth=axis_azimuth)

In [53]: (tr[['surface_tilt', 'surface_azimuth']] - tr2[['surface_tilt', 'surface_azimuth']]).describe()
Out[53]: 
       surface_tilt  surface_azimuth
count  1.780000e+02     1.780000e+02
mean  -6.586492e-16     3.193450e-15
std    8.916369e-15     2.864187e-14
min   -2.842171e-14    -5.684342e-14
25%   -7.105427e-15     0.000000e+00
50%    0.000000e+00     0.000000e+00
75%    3.552714e-15     2.842171e-14
max    2.131628e-14     5.684342e-14

kandersolar · 2022-06-17T17:08:19Z

Not that there was much doubt, but I've convinced myself that Bill's surface orientation equations are mathematically equivalent to the approach pvlib takes. Here are some notes if anyone is interested: https://gist.github.com/kanderso-nrel/ac3051de41261df317180c794144d6a9

If we do switch to Bill's equations we should be sure to preserve the handling of NaN and edge cases of the current implementation.

kurt-rhee · 2022-06-18T01:09:47Z

I realize that my small error metric is due to a small timeshift that I had in my data that changed my answer when resampling / averaging.

Cheers

kurt-rhee · 2022-07-05T21:39:23Z

@kanderso-nrel

I tried this code with a negative axis tilt and had some issues:

import pvlib
import pandas as pd
import numpy as np

axis_tilt = -5
axis_azimuth = 180
loc = pvlib.location.Location(40, -80)
times = pd.date_range('2019-06-01', '2019-06-02', freq='5T', tz='Etc/GMT+5')
sp = loc.get_solarposition(times)
tr = pvlib.tracking.singleaxis(sp.apparent_zenith, sp.azimuth,
                               axis_tilt=axis_tilt, axis_azimuth=axis_azimuth)


def rotation_to_orientation(tracker_theta, axis_tilt=0, axis_azimuth=0, max_angle=90):
    surface_tilt = np.rad2deg(
        np.arccos(
            np.cos(np.deg2rad(tracker_theta)) * np.cos(np.deg2rad(axis_tilt))
        )
    )
    surface_azimuth = np.rad2deg(
        np.deg2rad(axis_azimuth) +
        np.arcsin(
            (
                np.sin(np.deg2rad(tracker_theta)) /
                np.sin(np.deg2rad(surface_tilt))
            ).clip(upper=1, lower=-1)
        )
    )
    return pd.DataFrame({
        'tracker_theta': tracker_theta,
        'surface_tilt': surface_tilt,
        'surface_azimuth': surface_azimuth,
    })

tr2 = rotation_to_orientation(tr.tracker_theta, axis_tilt=axis_tilt, axis_azimuth=axis_azimuth)

print(
(tr[['surface_tilt', 'surface_azimuth']] - tr2[['surface_tilt', 'surface_azimuth']]).describe()
)

       surface_tilt  surface_azimuth
count  1.780000e+02       178.000000
mean  -1.741428e-15       -21.709086
std    1.299374e-14        64.389051
min   -5.773160e-14      -276.967526
25%   -7.105427e-15       -20.887249
50%    0.000000e+00       -10.025403
75%    4.218847e-15         7.843488
max    6.217249e-14        72.438867

kurt-rhee · 2022-07-05T22:13:23Z

@kanderso-nrel

Here is the difference between the two

Surface azimuth as calculated by pvlib in blue
Surface azimuth as calculated by the above script in orange

kandersolar · 2022-07-06T00:21:20Z

Oh interesting, thanks for pointing that out. Clearly for a tracker inclined slightly facing north as I interpret axis_tilt = -5 and axis_azimuth = 180 to mean, surface_azimuth should be 0/360 at solar noon. But negative axis tilt isn't something I've considered before. I've always thought of axis_tilt as being in the range [0, 90], and "negative" axis tilt is achieved by keeping axis_tilt positive and adding 180 to axis_azimuth. Bill's technical report defines axis tilt as being in [0, 90] as well.

As a point of reference, SAM/SSC claims MIN=0/MAX=90 for tilt as well, although "tilt" is overloaded there to refer to both fixed-tilt racking and SAT, so not an exact comparison. Interestingly though SAM does accept negative tilts for SAT simulations, and produces the same (incorrect) surface_azimuth curve:

@kurt-rhee do you think people will want to use pvlib with negative axis_tilt values? I'm wondering if this is a valid use case worth supporting or if we should just document that axis_tilt must be >= 0 and call it a day. I slightly object to negative axis_tilt on principle, but maybe that view isn't universal.

kurt-rhee · 2022-07-06T14:09:37Z

.

kurt-rhee · 2022-07-06T14:39:22Z

Hey @kanderso-nrel,

Thanks for that, very illuminating and clear as always. I think just a note in the documentation to keep axis_tilt >= 0 is enough. Better to keep 1 right way of doing things to avoid confusion.

Another thing I noticed about the calculation is that keeping axis_tilt positive and rotating axis azimuth by 180 degrees gives different answers using the Marion equations. This should be apparent given the below equation where azimuth is added to, but I was curious about the difference between pvlib's implementation and Marion's at axis azimuth >270 degrees. and <0 degrees.

surface_azimuth = np.rad2deg(
    np.deg2rad(azimuth) + (
        np.arcsin(
            (
                np.sin(np.deg2rad(-tracker_theta)) /
                np.sin(np.deg2rad(surface_tilt))
            ).clip(upper=1, lower=-1)
        )
    )

The left column of these plots are surface angles calculated by pvlib, the right column are calculated by equations. The top row is at 0 degrees, the bottom row is at 360 degrees.

williamhobbs · 2022-07-06T14:42:16Z

I just submitted this issue for NREL SAM: NREL/SAM#1087.

cwhanse · 2022-07-06T14:43:13Z

@kurt-rhee are all four of those surface_azimuth curves equivalent (once the angle is constrained to be in [0, 360))? Maybe I'm missing something.

kandersolar · 2022-07-06T14:55:02Z

@kurt-rhee instead of the quick mockup code from earlier in this thread, I suggest testing the implementation in #1480: https://github.com/pvlib/pvlib-python/pull/1480/files#diff-85ef163d90257651a76daf9d07052201a78139739b429ea96947451caad3dfb8R470-R521

One relevant difference is that surface_azimuth is constrained to [0, 360) as @cwhanse mentions.

kurt-rhee · 2022-07-06T14:57:07Z

@cwhanse They are equivalent! Thank you for pointing that out. Calculating POA using both methods gives slightly different answers.

@kanderso-nrel Nice, I will try that out right now.

(poa - poa2).describe()
Out[16]: 
         poa_global    poa_direct  ...  poa_sky_diffuse  poa_ground_diffuse
count  1.600000e+01  1.600000e+01  ...      8760.000000        1.600000e+01
mean   1.222035e+00  9.835606e-01  ...         0.000436       -5.551115e-16
std    4.888140e+00  3.934242e+00  ...         0.040767        2.220446e-15
min    0.000000e+00  0.000000e+00  ...         0.000000       -8.881784e-15
25%    0.000000e+00  0.000000e+00  ...         0.000000        0.000000e+00
50%    0.000000e+00  0.000000e+00  ...         0.000000        0.000000e+00
75%    1.847411e-13  1.278977e-13  ...         0.000000        0.000000e+00
max    1.955256e+01  1.573697e+01  ...         3.815592        0.000000e+00

kurt-rhee · 2022-07-06T15:48:19Z

@kanderso-nrel Works perfectly!

(poa - poa2).describe()
Out[4]: 
         poa_global    poa_direct  ...  poa_sky_diffuse  poa_ground_diffuse
count  1.600000e+01  1.600000e+01  ...     8.760000e+03        1.600000e+01
mean   6.039613e-14  4.263256e-14  ...     3.568936e-17       -5.551115e-16
std    9.617295e-14  7.903757e-14  ...     1.391200e-15        2.220446e-15
min    0.000000e+00  0.000000e+00  ...     0.000000e+00       -8.881784e-15
25%    0.000000e+00  0.000000e+00  ...     0.000000e+00        0.000000e+00
50%    0.000000e+00  0.000000e+00  ...     0.000000e+00        0.000000e+00
75%    1.278977e-13  2.842171e-14  ...     0.000000e+00        0.000000e+00
max    2.273737e-13  2.273737e-13  ...     5.684342e-14        0.000000e+00

kandersolar · 2022-07-06T15:56:48Z

Great, thanks for checking! Any interest in submitting a PR for clarifying the >=0 requirement in the pvlib.tracking docstrings?

kurt-rhee · 2022-07-06T16:28:12Z

I would love to!

kandersolar added the enhancement label Jun 15, 2022

kandersolar added this to the 0.9.2 milestone Jun 22, 2022

kandersolar mentioned this issue Jun 22, 2022

Move surface orientation calculation from tracking.singleaxis to new function; switch to Marion & Dobos 2013 equations #1480

Merged

8 tasks

kandersolar closed this as completed in #1480 Jul 5, 2022

williamhobbs mentioned this issue Jul 6, 2022

Sign (negative) ignored for 1 axis tilt NREL/SAM#1087

Closed

kurt-rhee mentioned this issue Jul 6, 2022

docstring updated to indicate limits of axis_tilt #1489

Merged

1 task

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Consider extracting the surface orientation calculation in pvlib.tracking.singleaxis() to its own function #1471

Consider extracting the surface orientation calculation in pvlib.tracking.singleaxis() to its own function #1471

kandersolar commented Jun 14, 2022

williamhobbs commented Jun 14, 2022

kurt-rhee commented Jun 16, 2022 •

edited

Loading

kandersolar commented Jun 17, 2022 •

edited

Loading

kandersolar commented Jun 17, 2022

kurt-rhee commented Jun 18, 2022

kurt-rhee commented Jul 5, 2022 •

edited

Loading

kurt-rhee commented Jul 5, 2022

kandersolar commented Jul 6, 2022

kurt-rhee commented Jul 6, 2022 via email •

edited

Loading

kurt-rhee commented Jul 6, 2022

williamhobbs commented Jul 6, 2022

cwhanse commented Jul 6, 2022

kandersolar commented Jul 6, 2022

kurt-rhee commented Jul 6, 2022

kurt-rhee commented Jul 6, 2022

kandersolar commented Jul 6, 2022

kurt-rhee commented Jul 6, 2022

Consider extracting the surface orientation calculation in pvlib.tracking.singleaxis() to its own function #1471

Consider extracting the surface orientation calculation in pvlib.tracking.singleaxis() to its own function #1471

Comments

kandersolar commented Jun 14, 2022

williamhobbs commented Jun 14, 2022

kurt-rhee commented Jun 16, 2022 • edited Loading

kandersolar commented Jun 17, 2022 • edited Loading

kandersolar commented Jun 17, 2022

kurt-rhee commented Jun 18, 2022

kurt-rhee commented Jul 5, 2022 • edited Loading

kurt-rhee commented Jul 5, 2022

kandersolar commented Jul 6, 2022

kurt-rhee commented Jul 6, 2022 via email • edited Loading

kurt-rhee commented Jul 6, 2022

williamhobbs commented Jul 6, 2022

cwhanse commented Jul 6, 2022

kandersolar commented Jul 6, 2022

kurt-rhee commented Jul 6, 2022

kurt-rhee commented Jul 6, 2022

kandersolar commented Jul 6, 2022

kurt-rhee commented Jul 6, 2022

kurt-rhee commented Jun 16, 2022 •

edited

Loading

kandersolar commented Jun 17, 2022 •

edited

Loading

kurt-rhee commented Jul 5, 2022 •

edited

Loading

kurt-rhee commented Jul 6, 2022 via email •

edited

Loading