Skip to content

Safely calculate dni from dhi and ghi #250

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

Merged
merged 44 commits into from
May 29, 2017
Merged

Safely calculate dni from dhi and ghi #250

merged 44 commits into from
May 29, 2017

Conversation

uvchik
Copy link
Contributor

@uvchik uvchik commented Oct 18, 2016

This PR is a little bit confusing, because it works on top of PR #239. The changes of this PR start with commit 6907b6c. It will be easier to read once PR #239 is merged.

The aim of this PR is to safely calculate DNI from GHI and DHI. There are different approaches and we may collect them within this function. A method parameter (or something similar) can be used to switch between the different approaches.

As this PR is a spin-off from PR #239 there are already two comments regarding this PR.

@adriesse wrote:

For your specific example, you can keep DNI values somewhat sane by not allowing them to be greater than a value obtained from a clear-sky model. Calculating the sun angles for the midpoint of the averaging interval is common practice (I hope). For long time steps like one hour, you can split the intervals that includes sunrise/set into two parts (one before and one after the event).

@cwhanse wrote:

Ø Calculating the angles for the midpoint of the averaging interval is common practice (I hope). It ought to be common practice. SAM does it, as does PVsyst I think, but other applications? We ought to build this capability into pvlib I think, as part of model chain – a kwarg switch perhaps?

A rough approach is to use the horizontal value for low angles, as there might not be relevant irradiation at small angles (high zenith angles) or set it to zero.

dni = (ghi - dhi) / tools.cosd(zenith)
dni[zenith > 88] = ghi[zenith > 88] - dhi[zenith > 88]
# or
dni[zenith > 88] = 0
  1. What methods do we want to implement?
  2. How do we want to switch between methods?

@birgits will help to implement this feature.

@uvchik uvchik force-pushed the master branch 2 times, most recently from e3eca71 to b9f6365 Compare November 24, 2016 13:42
@wholmgren wholmgren modified the milestones: 0.4.3, 0.4.2 Dec 2, 2016
@wholmgren
Copy link
Member

sorry for overlooking this @uvchik. I'd suggest some kind of filtering for the simple calculation but you could convince me that it's the user's responsibility. The comments above suggested to me that you were thinking of doing something more for time averages and sunrise/sunset transitions. Is that the case?

@birgits
Copy link
Contributor

birgits commented Dec 17, 2016

Sorry for having neglected this PR for so long. I have been working on it and am still having some problems, but I will post what I've got so far at the beginning of next week.

The comments above suggested to me that you were thinking of doing something more for time averages and sunrise/sunset transitions.

This is right. When calculating the DNI from GHI and DHI it may happen that for zenith angles close to 90°, so sunrise/sunset, the calculated DNI is unreasonable high or negative. The method we want to implement will "correct" those wrong values using different approaches.

@birgits
Copy link
Contributor

birgits commented Dec 20, 2016

This is a first suggestion for how to correct calculated DNI values that are unreasonably high or even negative in sunrise/sunset transitions. For the dataset I used the two methods implemented right now give nearly the same results.
@uvchik and I decided that it would be good to let the user know for how many timesteps the DNI had to be corrected, maybe serving as an indication on how good the quality of the given weather dataset is. We also found out that if corrections were only conducted for sunrise or sunset hours, a likely reason is that the timezone is not correct, and decided to also give out a warning in that case.

@adriesse
Copy link
Member

I am afraid that you could see problematic DNI values well before 88 degrees with real measurements due to cosine and alignment errors. Maybe you could make a gradual constraint that goes from the clear sky value at 80 degrees and zero at 90 degrees? In real data you may also need to constrain the DNI to non-negative values.

@birgits
Copy link
Contributor

birgits commented Dec 21, 2016

Thank you @adriesse for your comments and suggestions, I will try to implement them!
I was also wondering, since you suggested to take the clearsky DNI as an upper bound to get "sane" values during sunrise/sunset transitions, if the clearsky DNI should always be higher than the calculated or measured DNI, as is the case for GHI? I couldn't really find anything on that.

@adriesse
Copy link
Member

Since the clear sky DNI you use here is from a model, it is going to represent some kind of average conditions, so certainly the true DNI could be a bit higher than the model values. I suppose you could even add 10% to the model value and still achieve your objective of clipping unrealistic values.

A separate comment: in some cases it might be preferable to replace suspicious values with nan's.

@birgits
Copy link
Contributor

birgits commented Dec 21, 2016

A separate comment: in some cases it might be preferable to replace suspicious values with nan's.

I could add a parameter that can be used to specify if values should be corrected or set to nan. What do you think?

About setting the clearsky DNI as upper bound - to me it seems a bit random to add 10% to the modelled clearsky value, but maybe this could also be an optional parameter?

@adriesse
Copy link
Member

I think an option for nan's would be nice. I have no opinion about the 10%.

@wholmgren
Copy link
Member

Thanks @birgits. I agree with @adriesse's comments. A couple of implementation things:

The function signature should use basic inputs such as clearsky_dni and zenith rather than a Location. This gives users more more flexibility, keeps the function as simple as possible, and makes it faster. I'm guessing that most users will already have that data available in their script or interpreter session, too.

I'd also prefer that we not add the timezone speculation code. Generally speaking, users have a responsibility to get their timezones right. However, we could consider a more robust, dedicated function for checking timezones in a separate issue/pr (there was also some talk of this in one of the IO package/module threads).

It might help if you (and others) provide some examples of some real data that you're trying to clean up.

@birgits
Copy link
Contributor

birgits commented Dec 22, 2016

The function signature should use basic inputs such as clearsky_dni and zenith rather than a Location.

I agree with you. I just thought it would be easier for the user to only have to create a location object and not have to worry about how to get the zenith and which clearsky model to use, but have the definition take care of that. It took me quite a while anyway to go through all the clearsky models and try to understand how I can use them and what inputs I needed and I was very happy when I found that the location object does all that for me already. But you are right, it will be more straight forward the other way.

@birgits
Copy link
Contributor

birgits commented Apr 7, 2017

I took out the different methods and instead added two new parameters for the zenith angles between which the unreasonable values are set to NaN or corrected, as suggested by @wholmgren and @uvchik. What do you think?

A few things I want to point out:

  • I chose the default values for lower_cutoff_zenith and upper_cutoff_zenith from an early comment by @adriesse. They are open for discussion.
  • Negative DNI values as well as nonzero values above the upper_cutoff_zenith (above 88 for example) are set to NaN. Before there was also the option to set them to zero.
  • In the modelchain the function irradiance.dni is now by default called with the calculated clearsky_dni.

@uvchik
Copy link
Contributor Author

uvchik commented Apr 10, 2017

  1. I like the dni-function.

  2. It might be controversial but I would remove the **kwargs attribute from both functions (dni and complete_irradiance) and change the following lines in complete_irradiance:

-                clearsky_dni=kwargs.get('clearsky_dni',
-                                        self.location.get_clearsky(times).dni),
-                clearsky_tolerance=kwargs.get('clearsky_tolerance', 1))
+                clearsky_dni=self.location.get_clearsky(times).dni),
+                clearsky_tolerance=1))

People should use this function or do their own stuff.

  1. A NaN to Zero is not important for me. Anyway we might want to add it to the ModelChain but not to the dni function.

  2. You have to rebase the actual pvlib-dev branch on top of yours to solve the conflict.

@birgits
Copy link
Contributor

birgits commented May 7, 2017

Thank you @uvchik for your comments. @wholmgren and @adriesse do you have any more comments or suggestions?

birgits added 3 commits May 10, 2017 17:26
The only conflict was that both me and the upstream master added a test
at the end of a test file so the solution was easy: just take both.
I had a mishap which removed these files from being version controlled.
This should fix it.
Again, the only conflict was that both me and the upstream master added
a test at the end of a test file
Solution: just take both.

I have no idea, what happened during the first merges. Everything should
be fine now.
@wholmgren
Copy link
Member

I am generally ok with these additions. I'll make minor comments later this week.

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.

here are the minor comments

Default: 88.

lower_cutoff_zenith : float
Zenith angle above which DNI values greater the clearsky DNI (plus
Copy link
Member

Choose a reason for hiding this comment

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

"plus" should be "times"

Copy link
Contributor

Choose a reason for hiding this comment

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

Done

"""

# calculate DNI
dni_tmp = (ghi - dhi) / tools.cosd(zenith)
Copy link
Member

Choose a reason for hiding this comment

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

why do you need dni_tmp and dni? maybe a holdover from a previous implementation?

# upper_cutoff_zenith that are greater than the clearsky tolerance
if clearsky_dni is not None:
dni[(zenith >= lower_cutoff_zenith) & (zenith <= upper_cutoff_zenith) &
(dni > (clearsky_dni * clearsky_tolerance))] = (clearsky_dni *
Copy link
Member

Choose a reason for hiding this comment

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

You're doing the same calculation twice, so it may be slightly faster and slightly more readable if you define and use max_dni = clearsky_dni * clearsky_tolerance.

Copy link
Contributor

Choose a reason for hiding this comment

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

Done

dni[dni < 0] = float('nan')

# set non-zero DNI values for zenith angles >= upper_cutoff_zenith to NaN
dni[(zenith >= upper_cutoff_zenith) & (dni != 0)] = float('nan')
Copy link
Member

Choose a reason for hiding this comment

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

The doc string suggested to me that the comparison should be strictly greater than, rather than greater than or equal to.

Copy link
Contributor

Choose a reason for hiding this comment

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

Would that work and be correct English?

Zenith angle from which on nonzero DNI values are set to NaN.

Copy link
Member

Choose a reason for hiding this comment

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

Zenith angles greater than or equal to upper_cutoff_zenith will be set to NaN.

# correct DNI values for zenith angles between lower_cutoff_zenith and
# upper_cutoff_zenith that are greater than the clearsky tolerance
if clearsky_dni is not None:
dni[(zenith >= lower_cutoff_zenith) & (zenith <= upper_cutoff_zenith) &
Copy link
Member

Choose a reason for hiding this comment

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

same as above regarding comparison type

@@ -2048,3 +2048,70 @@ def _get_dirint_coeffs():
[0.743440, 0.592190, 0.603060, 0.316930, 0.794390]]

return coeffs[1:, 1:, :, :]


def dni(ghi, dhi, zenith, clearsky_dni=None, clearsky_tolerance=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 we be better off with a more descriptive name? dni_from_components? I'm not sure.

Copy link
Contributor

Choose a reason for hiding this comment

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

I find the name dni_from_components a bit confusing since components could mean different things. Maybe we just name it after what it does dni_from_ghi_and_dhi?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

According to aoi we could keep dni, but I like get_dnior calculate_dni definitely better than dni_from_something. We haven't used it so far and in my opinion the parameters should be part of the API documentation and not of the name.

Copy link
Member

Choose a reason for hiding this comment

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

Seems that dni is moderately acceptable to all.



def dni(ghi, dhi, zenith, clearsky_dni=None, clearsky_tolerance=1,
upper_cutoff_zenith=88.0, lower_cutoff_zenith=80.0, **kwargs):
Copy link
Member

Choose a reason for hiding this comment

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

no kwargs here

Copy link
Contributor

Choose a reason for hiding this comment

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

Done

-------
dni : Series
The modeled direct normal irradiance.

Copy link
Member

Choose a reason for hiding this comment

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

no extra space needed

Copy link
Contributor

Choose a reason for hiding this comment

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

Done

@@ -623,7 +623,7 @@ def effective_irradiance_model(self):
fd*self.total_irrad['poa_diffuse'])
return self

def complete_irradiance(self, times=None, weather=None):
def complete_irradiance(self, times=None, weather=None, **kwargs):
Copy link
Member

Choose a reason for hiding this comment

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

as discussed elsewhere, probably no kwargs here, too.

Copy link
Contributor

Choose a reason for hiding this comment

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

Done.
I had to set a default value for the clearsky_tolerance and chose 1.1, again according to an early comment from @adriesse. Would you suggest to have the clearsky_tolerance as an input into the complete_irradiance function, so it can be changed by the user or should the user write his/her own function in case he/she wants to change those parameters?

@adriesse
Copy link
Member

Good work. It occurs to me that the upper and lower cutoff names might not be entirely intuitive. How about: zenith_threshold_for_clearsky_limit, zenith_threshold_for_nan or something like that?

@birgits
Copy link
Contributor

birgits commented May 15, 2017

Sounds good @adriesse! How about zenith_threshold_for_zero_dni for the upper_cutoff_zenith, since the DNI for zenith angels above that value should be zero?

@birgits
Copy link
Contributor

birgits commented May 15, 2017

And I don't know what to do about the failing tests... All tests pass on my computer. Does anyone know what's going wrong?

@wholmgren
Copy link
Member

see #326 regarding the failing tests. I will merge the fix and retrigger the build.

Birgit Schachler added 3 commits May 17, 2017 13:45
…d comments

lower_cutoff_zenith is renamed to zenith_threshold_for_clearsky_limit
upper_cutoff_zenith is renamed to zenith_threshold_for_zero_dni
@birgits
Copy link
Contributor

birgits commented May 24, 2017

Thank you for fixing the failing tests @wholmgren! Do you think the PR is ready to be merged?

@wholmgren
Copy link
Member

Looks good to me. Can you add a note and your name to the 0.4.5 whats new file?

@wholmgren wholmgren merged commit f8423c9 into pvlib:master May 29, 2017
@wholmgren
Copy link
Member

thanks @birgits and @uvchik!

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.

4 participants