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

Calculating horizon profiles and associated shading losses #758

Open
wants to merge 29 commits into
base: main
Choose a base branch
from

Conversation

JPalakapillyKWH
Copy link
Contributor

@JPalakapillyKWH JPalakapillyKWH commented Jul 29, 2019

  • I am familiar with the contributing guidelines.
  • Fully tested. Added and/or modified tests to ensure correct behavior for all reasonable inputs. Tests (usually) must pass on the TravisCI and Appveyor testing services.
  • Updates entries to docs/sphinx/source/api.rst for API changes.
  • Adds description and name entries in the appropriate docs/sphinx/source/whatsnew file for all changes.
  • Code quality and style is sufficient. Passes LGTM and SticklerCI checks.
  • New code is fully documented. Includes sphinx/numpydoc compliant docstrings and comments in the code where necessary.
  • Pull request is nearly complete and ready for detailed review.

Code to get generate horizon profiles using the googlemaps API. Also has functionality to create shading losses due to the horizon profile. Currently only done with an isotropic sky diffuse model and the direct normal.
Follows much of the methodology in https://doi.org/10.1016/j.solener.2014.09.037

I've been developing this code in recent weeks. Is this something PVlib wants to incorporate?

@adriesse
Copy link
Member

I think this would make a great addition!

@cwhanse
Copy link
Member

cwhanse commented Jul 30, 2019

Is this something PVlib wants to incorporate?

Yes, I think so. horizon.py is OK, I would like to think a bit more about how this module relates to existing modules (irradiance.py, bifacial.py) and other process models not yet in pvlib, e.g., effects of partial shading on an array.

@cwhanse
Copy link
Member

cwhanse commented Jul 31, 2019

@JPalakapillyKWH I'm not familiar with Google maps elevation data. What is the spatial resolution? Do the data represent the near-field or only the far-field horizon? Near-field accounts for trees, buildings, etc., far-field is terrain such as mountains and hills.

@JPalakapillyKWH
Copy link
Contributor Author

@cwhanse The spatial resolution varies since googlemaps (as I've gathered) uses a patchwork of different signals but is generally at least 90m. The code will capture the near-field horizon to the extent that the googlemaps AP is affected by buildings/trees etc.

I included the googlemaps elevation calls because it was a simple way to get elevation data for sites distributed across the world, but the bulk of the horizon module is agnostic to the source of the elevation data.

@cwhanse
Copy link
Member

cwhanse commented Jul 31, 2019

Ok, was just thinking how to describe the output of this module, as near-field or far-field shading. Horizon is usually interpreted as far field.

@JPalakapillyKWH
Copy link
Contributor Author

It is designed for far-field.

@cwhanse
Copy link
Member

cwhanse commented Aug 1, 2019

@JPalakapillyKWH super content for pvlib. I've read the reference and look through the code. Some high-level comments about modules and connecting with the rest of pvlib. would appreciate @wholmgren @mikofski input here.

  • horizon.py main purpose is to provide functions that return the horizon profile (list of (azimuth, dip angle) tuples) and sky diffuse loss factor dtf. I'd put calculate_dtf and collection_plane_dip_angle in horizon.py. I'd add the Boolean "in-front" variable as an output if that is being calculated somewhere.
  • Some of the helper functions e.g. pol2cart, latitude_to_geocentric could go into pvlib.tools. I'd leave the grid, sample, filter and gmaps functions in horizon.py.
  • There's a mix of "to" and "from" function names, I'd prefer "to".
  • I am wondering if there's a way to use a modified calculation of dtf in the current isotropic sky diffuse model, rather than making a new model horizon_adjusted. Could we change the denominator in the dtf calculation to sum over only the cells in front of the array? Then sky diffuse = dhi * (1 + cos(tilt))/2 * dtf and dtf could be an optional kwarg for the current isotropic function.
  • it's too bad the reference isn't more explicit about the time series simulation.

@wholmgren
Copy link
Member

I skimmed the code and it looks like a great start - thanks @JPalakapillyKWH! I've not yet looked at the reference.

In addition to expanding some of the function documentation, I think we're going to want some documentation that synthesizes the proposed code.

There are a handful of places that generate tuples or loop over them in calculations. Can the code make more use of numpy arrays to improve performance and readability?

I would encourage marking some of the helper functions as private with an _ and/or putting some of these functions in pvlib.tools (which is not explicitly private but is also not documented in the api.rst file).

Great that the code is agnostic about elevation data source. I am not so sure about including the googlemaps api wrapper code in pvlib python, but I'm definitely ok using it in example documentation. I'd probably want to use the functions with a gridded data source such as SRTM.

@JPalakapillyKWH
Copy link
Contributor Author

@cwhanse and @wholmgren I made a lot of the code restructuring changes y'all suggested. Many functions were moved to tools.py.

I am currently working on making the code more numpy-friendly. I'm also working on creating documentation by creating a .rst file in the sphinx source folder and updating api.rst. Is this right way to create documentation?

I did include a few functions for visualizing horizon profiles but they require matplotlib to run. Do these seem worth including in the horizon module itself or is this something that would go in the documentation?

@cwhanse I don't think summing over the front of the panel will work since the horizon could be at a higher elevation angle than the back of the panel (e.g. panel tilt is 20 deg but the horizon is at a 30 deg elevation angle behind it). I just included the horizon adjustment code in the isotropic model with horizon profile as a kwarg. Let me know what you think.

pvlib/irradiance.py Outdated Show resolved Hide resolved
pvlib/irradiance.py Outdated Show resolved Hide resolved
pvlib/irradiance.py Outdated Show resolved Hide resolved
pvlib/irradiance.py Outdated Show resolved Hide resolved
pvlib/test/test_irradiance.py Outdated Show resolved Hide resolved
@mikofski
Copy link
Member

mikofski commented Aug 6, 2019

Sorry to jump in late. Great addition @JPalakapillyKWH.

One caveat that should probably be mentioned somewhere, is that I believe some irradiance sources, such as NSRDB may already include the far field horizon loss, since they are calibrated by ground stations, please correct me if I'm wrong.

I too would like to see examples from other geographical sources like NASA SRTM and ESRI, but I think those should be separate PRs .

I'll try to take a closer look later today.

Thanks!

@cwhanse
Copy link
Member

cwhanse commented Aug 6, 2019

@cwhanse I don't think summing over the front of the panel will work

I was not clear - I meant sum over the sky cells visible from the plane of array. The revised implementation accomplishes the same end for the isotropic sky and is arguably closer to the published reference.

Not in this PR, but it appears to me that we can extend the use of horizon data to other transposition models. That's what I had in mind when I made the comment.

@JPalakapillyKWH
Copy link
Contributor Author

@mikofski You are correct about ground station data already incorporating far field losses. This is designed for satellite irradiance data (PSM, TWC etc.)
@adriesse I took many of your suggestions and ran with them.

There's been a lot of restructuring in the code to move towards numpy. I'm still updating docstrings though.

@mikofski
Copy link
Member

mikofski commented Aug 6, 2019

ground station data already incorporating far field losses. This is designed for satellite irradiance data (PSM, TWC etc.)

Awesome! Would you consider adding a advisory note in the docstring with this warning on usage?

@JPalakapillyKWH
Copy link
Contributor Author

Awesome! Would you consider adding a advisory note in the docstring with this warning on usage?

Which docstring do you mean?

@mikofski
Copy link
Member

mikofski commented Aug 6, 2019

I would recommend adding a warning note in any method/function that uses the new horizon_* args/kwargs:

  • pvlib.irradiance.isotropic
  • pvlib.irradiance.get_sky_diffuse
  • pvlib.irradiance.get_total_irradiance

But this is not a blocker for me. It's just a suggestion, my opinion, take it or leave it :)

@JPalakapillyKWH
Copy link
Contributor Author

Working my way through the docstrings and code. I am a little confused about the purpose of the "sample" functions - if a grid is being supplied as input, why does the grid need to be sampled? Is the input grid not regularly spaced?

The input grid will be regularly spaced. The sampling functions are to produce a smooth horizon profile.
Without sampling, the "smoothness" of the horizon profile is highly dependent on the size of your grid. If you don't have enough points on your grid then you will have a very choppy horizon profile. Some azimuths might not even have elevation angles. So the sampling functions are a way of pseudo-generating more elevation data to make the horizon smoother.

@wholmgren
Copy link
Member

I stumbled across earthpy and this hillshade example. Thought people interested in this PR might want to check it out.

@cwhanse
Copy link
Member

cwhanse commented Aug 9, 2019

Working my way through the docstrings and code. I am a little confused about the purpose of the "sample" functions - if a grid is being supplied as input, why does the grid need to be sampled? Is the input grid not regularly spaced?

The input grid will be regularly spaced. The sampling functions are to produce a smooth horizon profile.
Without sampling, the "smoothness" of the horizon profile is highly dependent on the size of your grid. If you don't have enough points on your grid then you will have a very choppy horizon profile. Some azimuths might not even have elevation angles. So the sampling functions are a way of pseudo-generating more elevation data to make the horizon smoother.

Makes sense now, thanks.

pvlib/horizon.py Outdated Show resolved Hide resolved
pvlib/horizon.py Outdated Show resolved Hide resolved
pvlib/horizon.py Outdated Show resolved Hide resolved
pvlib/horizon.py Outdated Show resolved Hide resolved
pvlib/horizon.py Outdated Show resolved Hide resolved
pvlib/horizon.py Outdated Show resolved Hide resolved
pvlib/horizon.py Outdated Show resolved Hide resolved
pvlib/horizon.py Show resolved Hide resolved
pvlib/horizon.py Outdated Show resolved Hide resolved
pvlib/horizon.py Outdated Show resolved Hide resolved
pvlib/horizon.py Outdated
"""
assert(horizon_azimuths.shape[0] == horizon_angles.shape[0])

wedges = {}
Copy link
Member

Choose a reason for hiding this comment

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

I'll guess this loop may be a bit slow. What about instead:

rounded_azimuths = tools.round_to_nearest(horizon_azimuths, bin_size)
bins = np.unique(rounded_azimuths)

filtered = np.column_stack((bins, np.nan * bins))

for i in range(filtered.shape[0]):
    idx = np.nonzero(rounded_azimuths == filtered[i, 0])
    filtered[i, 1] = np.max(horizon_angles[idx])

return filtered[:, 0], filtered[:, 1]

Now that I've thought that through, there will always be a bin at 0 degrees, so the user ought to be informed of this point. I'd add to the note about bin_size.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

A couple of questions on this.
Is np.nonzero necessary? I think you can just index with the boolean array.
Why will there always be a bin at 0 degrees?

Copy link
Member

Choose a reason for hiding this comment

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

A couple of questions on this.
Is np.nonzero necessary? I think you can just index with the boolean array.

You're probably right. I started with argwhere, decided it wasn't the right approach, stopped at nonzero without thinking further.

Why will there always be a bin at 0 degrees?

You are right again, there won't always be a bin at 0. But the bins will be a subset of 0, +/-bin_size, +/-2*bin_size, etc. That's what I meant to describe and what I think we shoudl tell the user.

Copy link
Member

Choose a reason for hiding this comment

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

Needs tests for the functions added to tools.py then I'm OK with the review.

normal = np.atleast_2d(np.stack([2*x1/a**2, 2*y1/a**2, 2*z1/b**2], axis=1))

# Take the dot product of corresponding vectors
dot = np.sum(np.multiply(delta, normal), axis=1)
Copy link
Member

Choose a reason for hiding this comment

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

Would np.dot() do the trick?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Okay so I'm taking two Nx3 arrays and trying to do an "element-wise" dot product to get back an Nx1 array. I couldn't get np.dot() to do this but let me know if you know how.

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 np.dot should work if the inputs are properly aligned, but let's address the potentially unnecessary array operations above before figuring out this one.

pvlib/horizon.py Outdated
"""
The ``horizon`` module contains functions for horizon profile modeling.
There are various geometric utilities that are useful in horizon calculations
as well as a method that uses the googlemaps elevation API to create a
Copy link
Member

Choose a reason for hiding this comment

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

Was the gmaps function removed?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yep. I'll fix the docstring

pvlib/horizon.py Show resolved Hide resolved
@adriesse
Copy link
Member

I read the paper mentioned at the top of this PR and my impression is that there is not much overlap between the methods discussed there, and those programmed here. Is this correct?

I would very much like to see notes in the docstrings regarding the nature of the algorithms used, their origins, and their validation. If there are new original methods, please also identify and take credit for them! If there is a forthcoming publication, I would love to see it; but if you're still working on one, I would suggest pausing the PR until it is finished and you can share it.

@JPalakapillyKWH
Copy link
Contributor Author

@adriesse
You are right in that there is not a whole lot of overlap. Much of the code in this PR is taken from the geometry used in viewshed analysis. The sampling functions and binning of the horizon profile were ideas that I got from http://sherrytowers.com/2014/04/13/archeoastronomy-calculating-the-horizon-profile-using-online-us-geographic-survey-data/.
I don't have any plans for a publication, just working on this for the summer.

For my own validation of horizon profile modeling, I used the algorithms in this module with elevation data taken from the googlemaps elevation API to calculate horizon profiles. These were mostly compared by eye to horizon profiles from other sources and the online tool http://www.heywhatsthat.com/.

I also did validation of shading losses due to the horizon. I found shading loss by comparing the output of a modelchain run with the horizon adjusted irradiance vs without. I then compared this to far field losses outputted by PVSyst for a few modeled projects and found that they agreed.

@adriesse
Copy link
Member

Perhaps a good next step would be to add (or expand) Notes sections for the doc strings, where you explain the algorithms used, and give references to the sources you used, where applicable.

If validation is not provided in some formal publication, then perhaps it can be done in a jupyter notebook or something, but I don't know what the best place would be to store such a thing. Of course you could also reconsider whether maybe you could make a publication out of this work...

@cwhanse
Copy link
Member

cwhanse commented Aug 13, 2019

Perhaps a good next step would be to add (or expand) Notes sections for the doc strings, where you explain the algorithms used, and give references to the sources you used, where applicable.

calculate_dtf can reference the Goss et al. Solar Energy paper for the idea of dtf. We need a reference for the summarization formula (lines 598 - 601). I thought the code was implementing the sum over equal area sky patches? I may have missed a decision along the way.

I'd reference the ACM transactions paper instead of the class notes for uniformly_sample_triangle.

Add the reference to elevation_and_azimuth - it is called out as [1] but not listed. The rest of the functions seem straight-forward enough.

If validation is not provided in some formal publication, then perhaps it can be done in a jupyter notebook or something, but I don't know what the best place would be to store such a thing. Of course you could also reconsider whether maybe you could make a publication out of this work...

I vote for a page on readthedocs. Separate PR though.

Copy link
Member

@wholmgren wholmgren left a comment

Choose a reason for hiding this comment

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

@JPalakapillyKWH thanks for the significant refactoring. I think we're a lot closer. I like the scope of the PR now.

I agree that we need some kind of validation of the functions before we can merge, but this does not need to be a publication. It could be a jupyter notebook that you link to, an rst page that's built with our documentation, expansion of the module docstring, examples sections in the function docstrings, figures pasted into the comments in this thread, or some combination of the above.

lon : numeric
The longitude of the location that is to be the center of the grid.

grid_radius : numeric
Copy link
Member

Choose a reason for hiding this comment

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

"radius" is a surprising name for an argument in a function that returns a square grid. Not opposed to keeping it because I'm guessing that what we really would like is a circular grid.

Copy link
Member

Choose a reason for hiding this comment

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

Getting ahead of myself but perhaps the function could return nan for grid points outside grid radius (near the corners of the square). Probably should use a library like shapely or geopandas if we going to do something more than the trivial approach.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

A lot of elevation data (SRTM) is already gridded rectangularly. I figured a rectangular grid would be easiest to interface with different sources of elevation data.
I agree radius isn't a great term but I'm not sure what's better..."half-width"?

Copy link
Member

Choose a reason for hiding this comment

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

How about specifying width (=length) and resolution in km, and then making your rectangular lat/lon grid to suit?

elevation angles are to be calculated. Longitude should be given in
degrees East of the Prime Meridian and latitude in degrees North of the
Equator. Units are [deg, deg, meters]

Copy link
Member

Choose a reason for hiding this comment

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

delete line

The elevation angles that the points in pt2 make with the horizontal
as observed from the points in pt1. Given in degrees above the
horizontal.

Copy link
Member

Choose a reason for hiding this comment

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

Brief Examples section would be helpful to avoid confusion about array formats.

pvlib/horizon.py Outdated
# Polar Radius of the Earth (ellipsoid model) in meters
b = 6356752.0

lat1 = np.atleast_1d(pt1.T[0])
Copy link
Member

Choose a reason for hiding this comment

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

why the calls to atleast_1d and atleast_2d here and below? We've specified that the input is a Nx3 ndarray, so I don't understand why this is necessary.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I had this function accept [3,] shaped arrays earlier.

pvlib/horizon.py Outdated

delta = np.atleast_2d(np.subtract(v1, v2))

normal = np.atleast_2d(np.stack([2*x1/a**2, 2*y1/a**2, 2*z1/b**2], axis=1))
Copy link
Member

Choose a reason for hiding this comment

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

define a_sqrd and b_sqrd instead of a, b to avoid repeated values here. Multiply the whole array by 2 instead of each component.

pvlib/horizon.py Outdated
return samples


def uniformly_sample_triangle(p1, p2, p3, num_samples):
Copy link
Member

Choose a reason for hiding this comment

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

What's the purpose of this public method?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

will make private

pvlib/tools.py Outdated
# Equatorial Radius of the Earth (ellipsoid model) in meters
a = 6378137.0
# Polar Radius of the Earth (ellipsoid model) in meters
b = 6356752.0
Copy link
Member

Choose a reason for hiding this comment

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

Let's define these once at the top of tools.py as earth_radius_equatorial and earth_radius_polar. You can set a = and b = within the functions if you want to keep the expressions short. Or a_sqrd = , b_sqrd = in the functions as suggested elsewhere.

setup.py Outdated
'pandas >= 0.18.1',
'pytz',
'requests']
'requests',
'scipy',
Copy link
Member

Choose a reason for hiding this comment

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

revert this change

setup.py Outdated
@@ -37,10 +37,12 @@
MAINTAINER_EMAIL = 'holmgren@email.arizona.edu'
URL = 'https://github.com/pvlib/pvlib-python'

INSTALL_REQUIRES = ['numpy >= 1.10.4',
INSTALL_REQUIRES = ['numpy >= 1.13.3',
Copy link
Member

Choose a reason for hiding this comment

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

what are you using that requires numpy 1.13? Increasing the required version is allowed, but I'd like to understand why.

pvlib/horizon.py Outdated
The ``horizon`` module contains functions for horizon profile modeling.
There are various geometric utilities that are useful in horizon calculations.
"""
from __future__ import division
Copy link
Member

Choose a reason for hiding this comment

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

python 3 only now, so this is not necessary

over the visible section of the sky dome. The integrand of the surface
integral is the cosine of the angle between the incoming radiation and the
vector normal to the surface. The method calculates a sum of integrations
from the "peak" of the sky dome down to the elevation angle of the horizon.
Copy link
Member

Choose a reason for hiding this comment

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

Could you mention the relevant section or equation(s) in [1]?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Do you mean [2]?

Copy link
Member

Choose a reason for hiding this comment

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

I only see one reference...

Copy link
Contributor Author

Choose a reason for hiding this comment

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

my bad. latest changes have the new reference. Old one was wrong.

Copy link
Member

Choose a reason for hiding this comment

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

Ok, I read this paper too. Which equation(s) are you implementing?


[1] Goss et al. (2014) Solar Energy 110, 410-419

[2] Wright D. (2019) IEEE Journal of Photovoltaics 9(2), 391-396
Copy link
Member

Choose a reason for hiding this comment

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

Looks interesting. Could you send me a copy?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

sent an email

@cwhanse
Copy link
Member

cwhanse commented Jan 31, 2020

@JPalakapillyKWH could you contact me privately at cwhanse@sandia.gov? I have a few a questions I'd appreciate your input.

@mikofski
Copy link
Member

I would like to see this topic picked up again, I added an issue #1290 to spur more interest, that issue also adds the need to download the horizon data (or figure out how to generate it form GIS elevation data). I couldn't tell if this PR already encompasses downloading horizon data or just applying it to existing solar resource.

@cwhanse
Copy link
Member

cwhanse commented Aug 25, 2021

Most of this PR involved turning a [lat, long, elev] grid into [azimuth, elevation] for an observer location. I don't see any code to download the topographic data, but I could have missed it.

The original submitter @JPalakapillyKWH isn't available to complete this PR, so feel free (anyone) to pick it up from here.

@mikofski
Copy link
Member

mikofski commented Aug 25, 2021

Maybe probably needs 3 or 4 separate tasks then?

  1. either get horizon data or elevation data
  2. if elevation data, convert it to horizon data
  3. use horizon data to calculate effect on DHI and DNI (I think? I might be wrong)

@AdamRJensen
Copy link
Member

Maybe probably needs 3 or 4 separate tasks then?

  1. either get horizon data or elevation data
  2. if elevation data, convert it to horizon data
  3. use horizon data to calculate effect on DHI and DNI (I think? I might be wrong)

Tasks 1 and 3 listed by @mikofski have been addressed in #1395 and #1849. I suggest we close this issue as it is unlikely to get merged and already somewhat addressed by other PRs in the mean time. Thoughts @pvlib/pvlib-core ?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

7 participants