Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update irradiance.aoi to use more reliable formula #1191

Merged
merged 10 commits into from
May 14, 2021

Conversation

kandersolar
Copy link
Member

I decided the proposed option 2 in #1185 made the most sense, but happy to go another way if a better fix exists.

@kandersolar kandersolar added this to the 0.9.0 milestone Mar 10, 2021
@cwhanse
Copy link
Member

cwhanse commented Mar 10, 2021

I don't think option 4 is feasible. There will always be some combination of inputs that produce a round off error.

My only thought about using np.clip is, what if somehow in some really weird way tools.cosd (an untested pvlib function) changes and fails to produce values in [-1, 1] within roundoff? We wouldn't know.

An alternative is to truncate projection to 15 decimal places, since the round off surely won't be greater than the minimum floating point number (about 2.2e-16). In the example you found, that is exactly the amount that projection exceeds 1. I think this is Option 1.

In any case, I think the risk is small with the simple np.clip in this PR, and will not object if this goes forward.

Or, use np.where in a sort of complicated way, to replace (projection > 1) & (projection <= 1.0 + eps) with 1.0.

@kandersolar kandersolar changed the title Clamp AOI projection values to [-1, +1] Update irradiance.aoi to use more reliable formula Mar 13, 2021
def aoi(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth):
"""
Calculates the angle of incidence of the solar vector on a surface.
This is the angle between the solar vector and the surface normal.

Input all angles in degrees.

.. versionchanged:: 0.9.0
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I didn't know about this directive. Maybe we could be using it more often.

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

for arg in [surface_tilt, surface_azimuth, solar_zenith, solar_azimuth]:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it's OK to assume that these Series have a common index, and stop the loop when one index if found.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

does the new implementation preserve Series in --> Series out without this code?

I'm fine with dropping the name assignment and letting users deal with that.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

does the new implementation preserve Series in --> Series out without this code?

Unfortunately it does not; the _spherical_to_cartesian function turns Series into 2-D numpy arrays. I don't immediately see a way around it that's cleaner than the "fix types at the end" idiom.

@@ -160,6 +160,12 @@ def aoi_projection(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth):

Input all angles in degrees.

.. warning::
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would not be opposed to retiring the aoi_projection function. aoi could return cos(aoi) also, if that's of interest.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

+1 on retiring aoi_projection, but let the few users that need cos(aoi) call it themselves.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's not a big deal, but several of the sky diffuse transposition models use aoi_projection(), so that deprecation will need some other minor changes too.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

-1 on retiring aoi_projection because I use it.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

thanks @adriesse. I suggest we keep it and clip the output. Perhaps we should document the formula used by each function.

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

for arg in [surface_tilt, surface_azimuth, solar_zenith, solar_azimuth]:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

does the new implementation preserve Series in --> Series out without this code?

I'm fine with dropping the name assignment and letting users deal with that.

@@ -160,6 +160,12 @@ def aoi_projection(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth):

Input all angles in degrees.

.. warning::
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

+1 on retiring aoi_projection, but let the few users that need cos(aoi) call it themselves.

@adriesse
Copy link
Member

Further to my request to keep the function, I searched my code base and I make relatively few uses of it, mainly in code where I did not want a dependency on my private toolbox.

As to my usage patterns:

  • I don't always need aoi when I need cosinc (I like short names)
  • But I always need to clip cosinc to positive numbers
  • As a result I usually first calculate cosinc, then aoi just in case I need it, then clip cosinc to positive

@wholmgren
Copy link
Member

Should we merge this as is? Should we also add clip to the aoi projection function? Should we forget about the clever math and just use clip? I'm ok with any of these options.

.. versionchanged:: 0.9.0
Updated to a slower but more reliable formula that correctly handles
cases where the surface normal and solar position vectors are parallel.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This method treats aforementioned vectors as sides of a triangle, a and b, each length one. The value c in the code below is the length of the third side.

aoi_value.name = 'aoi'
except AttributeError:
pass
aoi = 2 * np.arctan(c / np.sqrt((1+(1+c))*((1-c)+1)))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think since a and b both equal 1, the sequencing of operations is not necessary, and I would suggest:
np.arctan2(c, np.sqrt((2.0 + c) * (2.0 - c)))

@adriesse
Copy link
Member

After finally following the links and looking at the reference, I must say I found these thoroughly interesting!

However, I am not convinced that a triangle formed with two unit vectors runs as many risks of numerical issues as and arbitrary triangle with sides a, b and c, and forming the vectors in the first place requires substantially more imprecise trig function calls. Furthermore, this method produces the aoi directly, so for most purposes a cos() call must follow.

I vote for a few well-placed clips to keep values sane.

@cwhanse cwhanse mentioned this pull request May 5, 2021
24 tasks
@kandersolar
Copy link
Member Author

I vote for a few well-placed clips to keep values sane.

This PR started out with a single clip near the end of the aoi_projection function -- anyone object to going back to 8069470? I mildly prefer it over maintaining an arcane equation I only sort of understand.

@wholmgren
Copy link
Member

No objection.

@adriesse
Copy link
Member

No objection. It was nice to learn about this other formula.

zenith = 89.26778228223463
azimuth = 60.932028605997004
projection = irradiance.aoi_projection(zenith, azimuth, zenith, azimuth)
assert projection == 1
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't know enough to intelligently comment on this, but I noticed that we're clipping to the int values -1 and 1. Is projection here also an int? I believe that if tested with realistic array input we'd always see floats returned. Might be worth testing that the clip behavior remains as expected in that case. Perhaps the important thing is that projection is always <= 1 and that np.isclose(projection, 1) returns True.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good suggestions. I also added a dtype check to make sure that clipping with int limits doesn't do something funky.

@cwhanse
Copy link
Member

cwhanse commented May 14, 2021

Lets merge

@wholmgren wholmgren merged commit 5bdad64 into pvlib:master May 14, 2021
@kandersolar kandersolar deleted the aoi_projection_clamp branch May 14, 2021 17:14
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging this pull request may close these issues.

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