Pvlib-python: solarposition.solar_azimuth_analytical returns nan

[adriesse]

The value passed to arccos() can end up being greater than 1 due to numerical imprecision. When this happens, nan is returned.
return np.sign(hour_angle) * np.abs(np.arccos((np.cos(zenith) * np.sin( latitude) - np.sin(declination)) / (np.sin(zenith) * np.cos( latitude)))) + np.pi
My solution would be to put in a clip(-1,1) (and also do some reformatting).

A search for arccos and arcsin turns up a few other code locations that could be susceptible to this.

[adriesse]

This could be a good practice problem for my first PR...

[cwhanse]

I would test for a tolerance rather than clip, to alert the user if there's an error in the solar angles passed to this function.

[adriesse]

Actually it's not a problem with the inputs.

[cwhanse]

Yes, I understood that. I'm saying that the error can be caused by other conditions, which would be masked by a clip. I'm in favor of clip when we've tested a tolerance to determine that the error is due to numerical error, rather than other causes.

[cwhanse]

And +1 for the practice problem, I needed one of those myself :)

[wholmgren]

I'd like to understand this a little better. Is the issue that np.cos(zenith) * np.sin(latitude) - np.sin(declination) can be slightly outside of [0, 1] even for correct inputs with arbitrary numerical precision and you want to make this approximate function work for such inputs? If so, maybe something like

arg = np.cos(zenith) * np.sin(latitude) - np.sin(declination)
eps = 1e-6  # or whatever the maximum error is for reasonable inputs
arg = np.where(arg > 1, arg - eps, angle)

[adriesse]

Yes, that's correct. Your suggestion would be a less drastic fix indeed. np.round(arg, 6) would also do it.

Now I notice now that a divide by zero can also occur...

(By the way, that last angle in your code should be arg, right?)

[adriesse]

Here is some code to demonstrate the problem:

import numpy as np
from pvlib.solarposition import solar_azimuth_analytical, solar_zenith_analytical

lat  = 0
decl = np.linspace(-20, 20, 9)
ha   = np.linspace(-180, 180, 9)

ha, decl = np.meshgrid(ha, decl)

ha = ha.flatten()
decl = decl.flatten()

D2R = np.pi/180.

zen = solar_zenith_analytical(lat*D2R, ha*D2R, decl*D2R)
azi = solar_azimuth_analytical(lat*D2R, ha*D2R, decl*D2R, zen)

bad_zen = np.isnan(zen)
print('Bad zenith values:', bad_zen.sum() )

bad_azi = np.isnan(azi)
print('Bad azimuth values:', bad_azi.sum())

out = np.array([ha, decl, zen/D2R, azi/D2R]).T.round(5)
print(' [ha,  decl,  zen,  azi] ')
print (out[bad_azi])

[adriesse]

While putting together the list of test cases, it became apparent that solar_zenith_analytical() contributes to the problem. A case might have a calculated zenith angle of 20.00000000? degrees and produce a NaN value for azimuth. If I then put in the literal '20' for zenith, the NaN disappears. In case anyone ever wonders...

[wholmgren]

Is it the difference between float(20) and int(20) or that solar_zenith_analytical returns a float that is not quite equal to float(20)?

[adriesse]

solar_zenith_analytical returns a float that is not quite equal to float(20).

[adriesse]

closed by #431