irradiance.aoi can return NaN when module orientation is perfectly aligned with solar position

[kandersolar]

Describe the bug
I was playing with a dual-axis tracking mount with #1176 and found that when the modules are perfectly aligned with the sun (i.e. AOI should be exactly zero), floating point round-off can result in aoi projection values slightly greater than one, resulting in NaN aoi. This only happens for some perfectly-aligned inputs (for example tilt=zenith=20, azimuth=180 returns aoi=0 as expected).

To Reproduce

import pvlib
zenith = 89.26778228223463
azimuth = 60.932028605997004
print(pvlib.irradiance.aoi_projection(zenith, azimuth, zenith, azimuth))
print(pvlib.irradiance.aoi(zenith, azimuth, zenith, azimuth))

# output:
1.0000000000000002
RuntimeWarning: invalid value encountered in arccos:  aoi_value = np.rad2deg(np.arccos(projection))
nan

Expected behavior
I expect aoi=0 whenever module orientation and solar position angles are identical.

Versions:

• pvlib.__version__: 0.9.0-alpha.4+14.g61650e9
• pandas.__version__: 0.25.1
• numpy.__version__: 1.17.0
• python: 3.7.7 (default, May 6 2020, 11:45:54) [MSC v.1916 64 bit (AMD64)]

Additional context
Some ideas for fixes:

1. In irradiance.aoi_projection, return a hard-coded 1.0 for inputs within some small tolerance
2. In irradiance.aoi_projection, clamp return value to [-1, +1]
3. In irradiance.aoi, clamp aoi_projection values to [-1, +1] before calling arccos
4. Rework the irradiance.aoi_projection trig equations to not generate impossible values?

[wholmgren]

notable that pyephem essentially clips based on aoi projection.

I found this stackexchange answer after some googling. I must have made a mistake because I couldn't get that formula to work. But the other answer does work (see page 15 of the linked reference for original source).

from pvlib.tools import sind, cosd


def spherical_to_cartesian(zenith, azimuth):
    sin_zen = sind(zenith)
    return np.array([sin_zen*cosd(azimuth), sin_zen*sind(azimuth), cosd(zenith)])


def aoi(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth):
    solar_vec = spherical_to_cartesian(solar_zenith, solar_azimuth)
    surface_vec = spherical_to_cartesian(surface_tilt, surface_azimuth)

    c_vec = solar_vec - surface_vec

    c_sqrd = np.dot(c_vec, c_vec)
    c = np.sqrt(c_sqrd)

    aoi = 2 * np.arctan2(c_sqrd, (1+(1+c))*((1-c)+1))
    aoi_deg = np.rad2deg(aoi)
    
    #print(solar_vec, surface_vec, c_vec, c_sqrd, c, aoi, aoi_deg, sep='\n')
    
    return aoi_deg


def aoi_sum_diff(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth):
    solar_vec = spherical_to_cartesian(solar_zenith, solar_azimuth)
    surface_vec = spherical_to_cartesian(surface_tilt, surface_azimuth)

    c_diff = solar_vec - surface_vec
    c_sum = solar_vec + surface_vec

    aoi = 2 * np.arctan2(np.linalg.norm(c_diff), np.linalg.norm(c_sum))
    aoi_deg = np.rad2deg(aoi)
    
    #print(solar_vec, surface_vec, c_diff, c_sum, aoi, aoi_deg, sep='\n')
    
    return aoi_deg

zenith = 89.26778228223463
azimuth = 60.932028605997004
print(aoi(zenith, azimuth, zenith, azimuth))
print(aoi_sum_diff(zenith, azimuth, zenith, azimuth))
print(pvlib.irradiance.aoi(zenith, azimuth, zenith, azimuth))

0.0
0.0
nan

surface_tilt, surface_azimuth, solar_zenith, solar_azimuth = 0, 0, 90, 90
print(aoi(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth))
print(aoi_sum_diff(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth))
print(pvlib.irradiance.aoi(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth))

89.99999999999999
89.99999999999999
90.0

surface_tilt, surface_azimuth, solar_zenith, solar_azimuth = 45, 180, 0, 180
print(aoi(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth))
print(aoi_sum_diff(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth))
print(pvlib.irradiance.aoi(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth))

19.471220634490688  # what did I do wrong?
45.0
45.0

I don't know what the best approach is, but maybe you can dig at it more.

[kandersolar]

Cool, thanks! The problem was a missing sqrt. Here's an updated version with some additional tweaks so that it works for arrays in addition to scalars:

Click to expand!

import pvlib
from pvlib.tools import sind, cosd
import numpy as np


def spherical_to_cartesian(zenith, azimuth):
    sin_zen = sind(zenith)
    return np.array([sin_zen*cosd(azimuth), sin_zen*sind(azimuth), cosd(zenith)])


def aoi(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth):
    solar_vec = spherical_to_cartesian(solar_zenith, solar_azimuth)
    surface_vec = spherical_to_cartesian(surface_tilt, surface_azimuth)

    c_vec = solar_vec - surface_vec
    c = np.linalg.norm(c_vec, axis=0)
    c_sqrd = np.sum(c_vec*c_vec, axis=0)  # couldn't get np.dot working for vectors
    c = np.sqrt(c_sqrd)

    aoi = 2 * np.arctan(c / np.sqrt((1+(1+c))*((1-c)+1)))
    aoi_deg = np.rad2deg(aoi)
    return aoi_deg


def aoi_sum_diff(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth):
    solar_vec = spherical_to_cartesian(solar_zenith, solar_azimuth)
    surface_vec = spherical_to_cartesian(surface_tilt, surface_azimuth)

    c_diff = solar_vec - surface_vec
    c_sum = solar_vec + surface_vec

    aoi = 2 * np.arctan2(np.linalg.norm(c_diff, axis=0), np.linalg.norm(c_sum, axis=0))
    aoi_deg = np.rad2deg(aoi)
    return aoi_deg

Testing the three functions with many random inputs (where solar position and module orientation are parallel) shows that the two methods from that SO post reliably return aoi=zero while the current pvlib function returns nan for ~3.7% of inputs. The downside is that the two SO functions are half as fast as the current pvlib function. Avoiding redundant conversions to radians (by not using cosd and sind) squeak a bit of extra speed but not much.

import time

def count_nan(x):
    return np.sum(np.isnan(x))

np.random.seed(0)
N = 1_000_000
surface_tilt = solar_zenith = np.random.uniform(0, 180, size=N)
surface_azimuth = solar_azimuth = np.random.uniform(0, 360, size=N)

for function in [aoi, aoi_sum_diff, pvlib.irradiance.aoi]:
    st = time.time()
    result = function(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth)
    ed = time.time()
    print(all(result == 0), count_nan(result), ed - st)

Output:

True 0 0.26395654678344727
True 0 0.2891671657562256
False 37653 0.12442183494567871

I don't love cutting execution speed in half, but it's still pretty fast. Ok with me to use one of these functions instead of np.clip. @cwhanse?

Edit to add: using this wouldn't fix the underlying issue of irradiance.aoi_projection returning values > 1. Is aoi_projection useful for anything other than calculating aoi?

[cwhanse]

I'll take slow, correct and reliable over fast, mostly correct and not reliable.

[mikofski]

do a benchmark before you merge this so we can document when the performance hit occured

[cwhanse]

I have a mild preference for the aoi version, it is somewhat more readable and looks more likely to be further optimized.

[adriesse]

I've been using clip for a long time. All these trig functions are numerical approximations and converting to and from radians doesn't help. I also use vector representations frequently, but I think it is entirely possible that these variations have other corner cases where they fail. In clip I trust. :)