Skip to content

ENH: geometric calculations for sunrise, sunset, and transit #583

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 19 commits into from
Oct 9, 2018

Conversation

mikofski
Copy link
Member

@mikofski mikofski commented Sep 19, 2018

… and sun transit times

pvlib python pull request guidelines

Thank you for your contribution to pvlib python! You may delete all of these instructions except for the list below.

You may submit a pull request with your code at any stage of completion.

The following items must be addressed before the code can be merged. Please don't hesitate to ask for help if you're unsure of how to accomplish any of the items below:

  • Closes issue add sunrise / sunset to Location #114
  • 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.

Brief description of the problem and proposed solution (if not already fully described in the issue linked to above):


def _hours(times, hour_angle, longitude, equation_of_time):
timezone = times.tz.utcoffset(times).total_seconds() / 3600.
return (hour_angle - longitude - equation_of_time/4.)/15. + 12. + timezone
Copy link

Choose a reason for hiding this comment

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

E226 missing whitespace around arithmetic operator

Copy link
Member Author

Choose a reason for hiding this comment

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

from pep8: whitespace in expressions:

If operators with different priorities are used, consider adding whitespace around the operators with the lowest priority(ies). Use your own judgment; however, never use more than one space, and always have the same amount of whitespace on both sides of a binary operator.

Copy link
Member

Choose a reason for hiding this comment

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

fine with me here and below.

sunrise = (sr.time() for sr in sst[0])
sunset = (ss.time() for ss in sst[1])
transit = (tr.time() for tr in sst[2])
sunrise_hours = [sr.hour + (sr.minute + sr.second/60.)/60.
Copy link

Choose a reason for hiding this comment

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

E226 missing whitespace around arithmetic operator

Copy link
Member Author

Choose a reason for hiding this comment

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

from pep8: whitespace in expressions:

If operators with different priorities are used, consider adding whitespace around the operators with the lowest priority(ies). Use your own judgment; however, never use more than one space, and always have the same amount of whitespace on both sides of a binary operator.

transit = (tr.time() for tr in sst[2])
sunrise_hours = [sr.hour + (sr.minute + sr.second/60.)/60.
for sr in sunrise]
sunset_hours = [ss.hour + (ss.minute + ss.second/60.)/60. for ss in sunset]
Copy link

Choose a reason for hiding this comment

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

E226 missing whitespace around arithmetic operator

Copy link
Member Author

Choose a reason for hiding this comment

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

from pep8: whitespace in expressions:

If operators with different priorities are used, consider adding whitespace around the operators with the lowest priority(ies). Use your own judgment; however, never use more than one space, and always have the same amount of whitespace on both sides of a binary operator.

sunrise_hours = [sr.hour + (sr.minute + sr.second/60.)/60.
for sr in sunrise]
sunset_hours = [ss.hour + (ss.minute + ss.second/60.)/60. for ss in sunset]
transit_hours = [tr.hour + (tr.minute + tr.second/60.)/60.
Copy link

Choose a reason for hiding this comment

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

E226 missing whitespace around arithmetic operator

Copy link
Member Author

Choose a reason for hiding this comment

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

from pep8: whitespace in expressions:

If operators with different priorities are used, consider adding whitespace around the operators with the lowest priority(ies). Use your own judgment; however, never use more than one space, and always have the same amount of whitespace on both sides of a binary operator.

@mikofski
Copy link
Member Author

@thunderfish24 ask an ye shall receive

@mikofski
Copy link
Member Author

note: test failure for py3.6 is unrelated

@adriesse
Copy link
Member

adriesse commented Sep 20, 2018

By the way, the same formula is in Duffie and Beckman and in probably many other solar textbooks!

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.

thanks mark!

@@ -1210,3 +1210,50 @@ def hour_angle(times, longitude, equation_of_time):
)).total_seconds() / 3600. for t in times])
timezone = times.tz.utcoffset(times).total_seconds() / 3600.
return 15. * (hours - 12. - timezone) + longitude + equation_of_time / 4.


def _hours(times, hour_angle, longitude, equation_of_time):
Copy link
Member

Choose a reason for hiding this comment

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

I'd like a more descriptive function name or a one line documentation for _hours and _times.


def _hours(times, hour_angle, longitude, equation_of_time):
timezone = times.tz.utcoffset(times).total_seconds() / 3600.
return (hour_angle - longitude - equation_of_time/4.)/15. + 12. + timezone
Copy link
Member

Choose a reason for hiding this comment

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

fine with me here and below.

def _times(times, hours):
return np.array([(dt.timedelta(hours=h) + t.tz.localize(
dt.datetime(t.year, t.month, t.day)
)) for h, t in zip(hours, times)])
Copy link
Member

Choose a reason for hiding this comment

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

This loop using pure python is likely to be painfully slow and defeat the purpose of an analytical approach.

Copy link
Member Author

Choose a reason for hiding this comment

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

I think I can use numpy arrays with timedelta64[s] here instead

Copy link
Member Author

Choose a reason for hiding this comment

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

Nice catch @wholmgren , I get about a 70x speed improvement (from 43.5 ms to 605us) to by using numpy datetime64 and timedelta64[ns]. I wonder where else we could use these? There is an oddity that hours is a Float64Index from pandas, not just floats, which bothers me, but ok for now.

See this script that shows improvements with numpy over looping.

In [1]: import pandas as pd

In [2]: import datetime as dt

In [3]: import numpy as np

In [4]: d = pd.DatetimeIndex(start='2013-01-01 00:00:00', end='2013-12-31 23:59:59', freq='H', tz='Etc/GMT+8')

In [5]: h = np.loadtxt('../sunrise_sunset_transit.txt', dtype=np.dtype([('sunrise', float), ('sunset', float), ('transit', float)]))

In [6]: sr = pd.Float64Index(h['sunrise'])

In [7]: %paste
def _times(times, hours):
    return np.array([(dt.timedelta(hours=h) + t.tz.localize(
        dt.datetime(t.year, t.month, t.day)
    )) for h, t in zip(hours, times)])

## -- End pasted text --

In [8]: %timeit _times(d, sr)
43.5 ms ± 354 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [9]: del _times

In [10]: %paste
def _times(times, hours):
    """helper converts hours from an array of floats to localized times"""
    tz_info = times.tz
    return pd.DatetimeIndex(
        times.tz_localize(None).values.astype('datetime64[D]')
        + (hours.values * 3600. * 1.e9).astype('timedelta64[ns]'), tz=tz_info)

## -- End pasted text --

In [11]: %timeit _times(d, sr)
605 µs ± 5.07 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

@@ -15,6 +16,9 @@
from conftest import (requires_ephem, needs_pandas_0_17,
requires_spa_c, requires_numba)

DIRNAME = os.path.dirname(__file__)
PROJDIR = os.path.dirname(DIRNAME)
DATADIR = os.path.join(PROJDIR, 'data')
Copy link
Member

Choose a reason for hiding this comment

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

probably better to put these in the function if they're not going to be reused elsewhere.

transit_hours = [tr.hour + (tr.minute + tr.second/60.)/60.
for tr in transit]
test_data_file = os.path.join(
DATADIR, 'sunrise_sunset_transit.txt')
Copy link
Member

Choose a reason for hiding this comment

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

do we need such a large data file to comprehensively test the algorithm? I am guessing 10-20 well chosen points is sufficient.

Copy link
Member Author

Choose a reason for hiding this comment

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

I guess not, but I thought others might want to use this, it's a full year of hourly sunrise/sunset/transit times from NREL SPA. It's about 600kb, is that too much? Whatever you think is best, maybe ditch it completely, and just test 5-10 values hardcoded in the test?

sunrise = _times(times, sunrise_hour)
sunset = _times(times, sunset_hour)
transit = _times(times, transit_hour)
return (sunrise, sunset, transit)
Copy link
Member

Choose a reason for hiding this comment

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

We should try to make the API consistent with the existing function. So this should be changed to support dict or DataFrame output or we can discuss changing the API of the other function elsewhere.

Copy link
Member Author

Choose a reason for hiding this comment

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

I was thinking that this was a lower level helper function, and that a higher level api function like get_sunrise would call it and based on the how argument, and then package up the return values into a DataFrame as needed. But happy to change this to output a DataFrame for now if that's what the consensus is.

Copy link
Member

Choose a reason for hiding this comment

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

I'm in favor of returning a pandas object rather than the tuple, and also, the intermediate function get_next_sunrise with a how argument that wraps the geometric, pyephem and spa methods. The input is a tz-aware datetime and location. get_sunrise is somewhat ambiguous: assume the datetime is solar noon at the location. Is the desired output the sunrise of the current day, or the next sunrise (the following day)? get_sunrise(date) clarifies things but requires the function to localize date (for pyephem) or determine the corresponding date in UTC time (for spa). For these reasons I'm suggesting the next approach.

To complete the function set we'd add get_previous_sunrise, get_next_sunset, etc.



def _hours(times, hour_angle, longitude, equation_of_time):
timezone = times.tz.utcoffset(times).total_seconds() / 3600.
Copy link
Member

Choose a reason for hiding this comment

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

hmm odd to see times as both the object and explicit argument of its own method.

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 is weird, but I think it's the only way, I broke it up into two lines that explains a bit more:

tz_info = times.tz  # this is the pytz timezone info object that has the utcoffset method
tz = tz_info.utcoffset(times)  # this method requires a datetime-like argument to operate on


Parameters
----------
times : :class:`pandas.DatetimeIndex`
Copy link
Member

Choose a reason for hiding this comment

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

please remove the rst markup here. it detracts from readability when inspecting the function signature outside of the rendered docs.

def sunrise_sunset_transit_analytical(times, latitude, longitude, declination,
equation_of_time):
"""
Analytical calculation of solar sunrise, sunset, and transit.
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 care for the term analytical but I don't have a better one to suggest. Perhaps the references mentioned have a better name for this?

Copy link
Member

Choose a reason for hiding this comment

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

We may want to clean up the naming convention in this module. We have prefixes 'spa' and 'nrel' for functions using SPA, we sometimes have 'pyephem' but not always, and use the suffix 'analytical' (I'm Ok with this) for some purely geometric methods.

I prefer _, e.g., declination_cooper69. where there are multiple options to calculate quantity.

Copy link
Member

Choose a reason for hiding this comment

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

"geometric" might be an alternative to "analytical"

Copy link
Member Author

Choose a reason for hiding this comment

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

I really like "geometric" but would we change zenith and azimuth to match?

Copy link
Member

@cwhanse cwhanse Sep 24, 2018

Choose a reason for hiding this comment

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

I'd say yes, to follow the naming pattern, tacitly labeling the collection of spherical trigonometry methods as the 'geometric' model.

Deprecate _analytical to 0.7, I'd say.

Copy link
Member Author

Choose a reason for hiding this comment

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

OK, I like it! Everyone else in agreement? It'll have to wait until after SPI, but I'll try to get back to this asap. Thanks!

Copy link
Member

Choose a reason for hiding this comment

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

I like "geometric". We should use "geometric" for the new function in this PR, but let's take on the deprecation of the existing functions in a separate PR.

* use hardcoded values for 10-20 expected values
* remove DIRNAME, etc
* make stickler happy
* no rst markup in numpydoc parameters
* add references to duffie/beckman and Frank Vignola
* get pytz timezone info first so it's not weird
* add returns to docstring
@@ -1,5 +1,6 @@
import calendar
import datetime
import os
Copy link

Choose a reason for hiding this comment

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

F401 'os' imported but unused

* add one-liners to explain helpers
* remove import os from test
tz_info = times.tz
return pd.DatetimeIndex(
times.tz_localize(None).values.astype('datetime64[D]')
+ (hours.values * 3600. * 1.e9).astype('timedelta64[ns]'), tz=tz_info)
Copy link

Choose a reason for hiding this comment

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

W503 line break before binary operator

Copy link
Member Author

Choose a reason for hiding this comment

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

pep8: Should a line break before or after a binary operator?

To solve this readability problem, mathematicians and their publishers follow the opposite convention. Donald Knuth explains the traditional rule in his Computers and Typesetting series:

“Although formulas within a paragraph always break after binary operations and relations, displayed formulas always break before binary operations”

In Python code, it is permissible to break before or after a binary operator, as long as the convention is consistent locally. For new code Knuth's style is suggested.

emphasis mine

@mikofski
Copy link
Member Author

Notes:

  • travis py35 failure was due to bad http request that timed out
  • the py27-min fails because in pandas-0.14, the minimum requirement, you can't reset a tz-aware datetime

do new features have to support the minimum requirements? I can see that patching old code should be backwards compatible, but if you are installing pvlib now to use an analytical sunrise/sunset equation, then you will have the newest pandas

@mikofski
Copy link
Member Author

mikofski commented Sep 24, 2018

And FYI: pandas-0.14.0 isn't even compatible with numpy-1.10.1 according to conda

  • numpy-1.8.2 is the newest version that works with pandas-0.14.0
  • the oldest version of pandas that works on numpy-1.10.1 is pandas-0.17.0

* if old pandas change tz directly, and hours is already a numpy array
* if new pandas then use tz_localize to strip tz, and change hours from
Float64Index to numpy array
* move operator to end to match existing style in pvlib
@mikofski
Copy link
Member Author

@wholmgren I am seeing some backwards compatibility issues on Travis with numpy and pandas minimum requirements

Since there are so many API changes in v0.6 can we consider maybe updating some of the minimum requirements now? Pandas is currently v0.23.4 and NumPy is at 1.15, keeping really old requirements seems limiting to me a little. Thanks.

@mikofski
Copy link
Member Author

Also, can we maybe consider dropping support for Python-3.4?

@wholmgren
Copy link
Member

Since there are so many API changes in v0.6 can we consider maybe updating some of the minimum requirements now? Pandas is currently v0.23.4 and NumPy is at 1.15, keeping really old requirements seems limiting to me a little.

I am all in favor of raising the minimums if there is a reason to do so. But I take these incompatibilities as sign that the approach might be too complicated and will be difficult to maintain in the future. Let's think about if this is really the implementation that we want.

Somewhere you commented something to the effect of the spec'ed minimum numpy/pandas versions are incompatible -- is it that they are incompatible or that Continuum doesn't supply pre-built binaries for them?

Also, can we maybe consider dropping support for Python-3.4?

We can consider this too but not in this pull request.

@wholmgren
Copy link
Member

My guess is that much of the logic can be simplified if we do some combination of the following

  1. assume the input time is localized or in UTC
  2. convert to ints instead of datetime64.
  3. utilize the same helper functions that irradiance.get_extra_radiation uses.

@mikofski
Copy link
Member Author

it seems to be working now, only had to make 2 minor changes to satisfy the CI

  • ducktype hours for different versions of pandas that return either float or Float64index
  • cast the date back to nano-seconds for different versions of numpy that require safe casting of inputs to ufunc add

Somewhere you commented something to the effect of the spec'ed minimum numpy/pandas versions are incompatible -- is it that they are incompatible or that Continuum doesn't supply pre-built binaries for them?

I was wrong, you were right, according to setup.py pandas-0.14.0 min npy ver is 1.6.1 for py2 or 1.7.0 for py3

assume the input time is localized or in UTC

I am already assuming that the times must be tz-aware and that they have a tz or tz_info field

convert to ints instead of datetime64.

good idea, like Linux time, can we pursue this as a course of last resort? The calculations seem to be working now, I believe that internally that is what timedelta64 already does, since it is int already.

utilize the same helper functions that irradiance.get_extra_radiation uses.

Do you mean _handle_extra_radiation_types, I think this could be useful for handling the input if it's not a tz-aware pandas DatetimeIndex, but my view of the API is that would be handled by the get_* function for all how args, assuming that sunrise_sunset_transit_analytical is a low level function called by a get_sunrise... high level function. ok?

otherwise, (minus what's new) I think this is ready for review

@mikofski
Copy link
Member Author

mikofski commented Sep 24, 2018

I am all in favor of raising the minimums if there is a reason to do so.

BTW: pandas-0.14.0 is from 2014, IMO we should not support this because it's not a v1.0 release, so a version that is 4 years old may have more issues than just incompatibility.

But I take these incompatibilities as sign that the approach might be too complicated and will be difficult to maintain in the future. Let's think about if this is really the implementation that we want.

did you see this comment, the simpler approach was 70x slower, I don't want to create a maintenance headache, but I don't think this approach is too complicated.

Analytical calculation of solar sunrise, sunset, and transit.

.. warning:: The analytic form neglects the effect of atmospheric
refraction.
Copy link
Member

Choose a reason for hiding this comment

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

Probably "geometric" instead of analytic. Also guessing it neglects more than refraction, but haven't looked at the sources.

Copy link
Member

Choose a reason for hiding this comment

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

Also spherical earth, sun is represented by a point at its center, ...

Copy link
Member

Choose a reason for hiding this comment

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

Here and in the other "analytical"/"geometrical" functions should we give an estimate of the typical magnitude of the errors? Something like "The error depends on location and time of year but is of order 10 minutes."

Copy link
Member Author

Choose a reason for hiding this comment

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

from the old test the max differences are

  • sunrise 0.11 hours
  • sunset 0.11 hours
  • transit 0.02 hours

this is a rough estimate, only to 2 decimals, I manually iterated the 2nd decimal to see what tolerance I could get.

Also the errors are one sided, ie: for sunrise the geometric estimate is always earlier up to 0.11 hours, and down to some positive, non-zero number, so there's always some sunrise error that is a minute or so earlier than the apparent sunrise, and ie: for sunset it's always later up to 0.11 hours, and down to a positive, non-zero number, so there's always some sunset error that is a minute or so later than the apparent sunset.

tz_info = times.tz
try:
hours = hours.values # make a new reference
except AttributeError:
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 understand the point of this try/except. Is hours not always the same type? If not, suggest we change _hours() to do the right thing with something like

hours = # your eqn
# enforce Series
hours = pd.Series(hours, index=times)
# OR enforce array
hours = np.array(hours)
return hours

Then we can remove this try except.

Copy link
Member Author

Choose a reason for hiding this comment

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

no apparently depending on the pandas version it changes, sometimes it's a numpy array of floats and sometimes it's a pandas series of pandas.Float64Index, I didn't bother to track down which versions are which, but it seems to be in the 0.teen versions, I think >0.20 yields series. I'm sure I could enforce the type somewhere else, as an output of some calculation, but I didn't bother to

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 would simplify things to enforce it elsewhere. I am sorry for being slow/stubborn but these functions are hard for me to follow. I assume it will be even harder for future me to understand them when something inevitably breaks.

try:
# OLD pandas, "TypeError: Already tz-aware, use tz_convert to convert."
times.tz = None
except AttributeError:
Copy link
Member

@wholmgren wholmgren Sep 24, 2018

Choose a reason for hiding this comment

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

I also don't understand this try/except. Your API declaration requires a localized DatetimeIndex, so times.tz_convert is guaranteed to work. This is true for pandas 0.14 and 0.23.

Copy link
Member Author

@mikofski mikofski Sep 24, 2018

Choose a reason for hiding this comment

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

tz_convert changes the timezone from one to another, ie: it applies the tz_utcoffset so eg: from Etc/GMT+8 to UTC would add 8 hours to all of the times.

But that's not what we want. We just want to make the tz-aware time a local naive time, so we want to effectively replace tz_info with None.

Again in some older version of pandas, you could just do this by setting tz=None, and in newer ones you can call tz_localize(None) to convert a tz-aware time to a naive one, without applying a utc offset.

EG: 2017-01-01T12:30:00-0800 just becomes 2017-01-01T12:30:00 same time, but no more timezone

Copy link
Member

Choose a reason for hiding this comment

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

I understand now, thanks.

newer ones you can call tz_localize(None) to convert a tz-aware time to a naive one, without applying a utc offset.

Is this a documented behavior of tz_localize? I can't find official support for it. If not I recommend a different implementation.

Copy link
Member Author

Choose a reason for hiding this comment

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

I believe it is documented, I've used this before, I don't recall where, but I will look for it.

@cwhanse and @wholmgren thanks so much for being patient and providing really great feedback, most of which I believe I have addressed and incorporated.

  • enforce hours as numpy array in _hours() instead of in _times()
  • break up calculation into separate steps so they are more clear
  • use "geometric" instead of "analytical"
  • in warning, explain circular orbit with sun as center point, and state max errors of 6 minutes or so

Copy link
Member Author

@mikofski mikofski Sep 24, 2018

Choose a reason for hiding this comment

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

@wholmgren it's in the pandas documentation for DatetimeIndex.tz_localize

Passing None will remove the time zone information preserving local time.

It's and also here:

Passing None will remove the time zone information preserving local time.

and there

None will remove timezone holding local time.

but neither here nor there

@cwhanse normalize seems very interesting, and perhaps what we're both looking for? Since v0.16, see what's new.

Convert times to midnight.
The time component of the date-timeise converted to midnight i.e. 00:00:00. This is useful in cases, when the time does not matter. Length is unaltered. The timezones are unaffected.

* change name in solarposition.py and test_solarposition.py
* enforce that hours is a numpy array in _hours instead of in _times
* add more useful comments, separate into steps to explain process
@mikofski
Copy link
Member Author

thanks for feedback!, tests are passing! all good now?

@wholmgren
Copy link
Member

wholmgren commented Sep 24, 2018 via email

@wholmgren
Copy link
Member

If I understand correctly, pandas 0.16 adds normalize for Timestamps and Series. 0.14 does have it for DatetimeIndex:

>>> times = pd.DatetimeIndex(start='20180924 1600', end='20180924 1700', freq='1H', tz='America/Phoenix')
>>> times
<class 'pandas.tseries.index.DatetimeIndex'>
[2018-09-24 16:00:00-07:00, 2018-09-24 17:00:00-07:00]
Length: 2, Freq: H, Timezone: America/Phoenix
>>> times.normalize()
<class 'pandas.tseries.index.DatetimeIndex'>
[2018-09-24 00:00:00-07:00, 2018-09-24 00:00:00-07:00]
Length: 2, Freq: None, Timezone: America/Phoenix
>>> pd.__version__
'0.14.0'

* also remove utcoffset(), unreliable with pandas dataframes, only really
works with scalar datetimes and timestamps
* instead use tz_localise(None) to get naive, but still local time, and
subtract the localized time, to get a sequence of timezones
"""converts hour angles in degrees to hours as a numpy array"""
naive_times = times.tz_localize(None) # naive but still localized
tzs = 1 / (3600. * 1.e9) * (
naive_times.astype(np.int64) - times.astype(np.int64))
Copy link

Choose a reason for hiding this comment

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

E126 continuation line over-indented for hanging indent

@@ -718,3 +718,56 @@ def test_hour_angle():
# sunrise: 4 seconds, sunset: 48 seconds, transit: 0.2 seconds
# but the differences may be due to other SPA input parameters
assert np.allclose(hours, expected)
@pytest.fixture()
Copy link

Choose a reason for hiding this comment

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

E302 expected 2 blank lines, found 0

* remove geometric expected sunrise/sunset/transite test fixture
* check abs tolerance between geometric and spa are less than max absolute
error difference which is typically 4-6 minutes for sunrise/sunset, and
about 10-40 seconds for transite
* move _times_to_hours to solarposition, and use the same approach as
other hour/time conversions, use pandas>= tz_localize(None) or normalize
to get the local, naive times, and time at previous midnight, then convert
to np.int64, and then convert to hours from nanoseconds
* since fixture is a dataframe, get column, convert to MST timezone, before
calling _times_to_hours(<pd.DatetimeIndex>)
* get times from expected_rise_set_spa.index
* get lat/lon from golden_mst test_fixture
@@ -531,7 +531,9 @@ def test_get_solarposition_altitude(altitude, expected, golden):
"delta_t, method, expected", [
(None, 'nrel_numpy', expected_solpos_multi()),
(67.0, 'nrel_numpy', expected_solpos_multi()),
pytest.mark.xfail(raises=ValueError, reason = 'spa.calculate_deltat not implemented for numba yet')
pytest.mark.xfail(
Copy link

Choose a reason for hiding this comment

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

E122 continuation line missing indentation or outdented

pytest.mark.xfail(raises=ValueError, reason = 'spa.calculate_deltat not implemented for numba yet')
pytest.mark.xfail(
raises=ValueError,
reason = 'spa.calculate_deltat not implemented for numba yet')
Copy link

Choose a reason for hiding this comment

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

E251 unexpected spaces around keyword / parameter equals

@mikofski
Copy link
Member Author

mikofski commented Oct 7, 2018

@wholmgren , @cwhanse , ready to merge I think,

  • I reused the SPA fixture, and set atol to be roughly 6-minutes for sunset/sunrise, and about 1-minute for transit,
  • also moved _times_to_hours from test_solarposition.py to solarposition.py and changed to the approach Will suggested (pandas>=0.15, normalize, and astype(np.int64)) which seems to work best
  • had to finagle the columns in expected_rise_set_spa from numpy datetimes at UTC to the pandas DatetimeIndex localized to golden_MST
  • still seeing the build error in py27-min which is still set to use pandas==0.14

@mikofski
Copy link
Member Author

mikofski commented Oct 7, 2018

question: should we define

NS_PER_HR = 1.e9 * 3600.

somewhere, to save time? it's calculated in a few places now

@wholmgren
Copy link
Member

Thanks @mikofski. Ok with me if you want to define NS_PER_HR somewhere in solarposition.py.

I'd like to see this PR tested against the pandas 0.15 changes in #595. Probably easiest to wait until that goes in. I want to test that one myself, probably on Wednesday. Sorry for the backlog.

midnite or midnight? Seems like we should prefer the standard spelling.

I'm less confident working with the datetime64 types than with integers, but ok with me if you think it should work reliably.

Can you revise the name for this PR?

Also a few items on the checklist to be completed.

@mikofski mikofski changed the title add analytical methods from Frank Vignola's book for sunrise, sunset,… ENH: geometric calculations for sunrise, sunset, and transit Oct 8, 2018
* wrap long lines, link to issues
* add constant NS_PER_HR and substiture for (3600. * 1.e9)
* change arg name to hourangle avoid shadowing function hour_angle()
* fix spelling in _local_times_from_hours_since_midnight()
* use normalize() and astype(np.int64) instead of
naive_times.astype('datetime64[D]').astype('datetim64[ns]') and
midnight_times + (hours*NS_PER_HR).astype('timedelta64[ns]') for
stability and consistency with other hour/times conversion functions
# sunrise, sunset, and transit
return pd.DatetimeIndex(
(naive_times.normalize().astype(np.int64)
+ (hours * NS_PER_HR).astype(np.int64)).astype('datetime64[ns]'),
Copy link

Choose a reason for hiding this comment

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

W503 line break before binary operator

@mikofski
Copy link
Member Author

mikofski commented Oct 9, 2018

@wholmgren updated api and what's new docs, and removed numpy.timedelta64[ns] in favor of your suggested integer approach for stability with numpy and consistency with other hour-time conversion functions.

All tests pass in both environments with pandas>=0.15.0 & numpy>=1.10.1:

py27-min

virtualenv

[master] mikm@PTOD702634:~/Projects/pvlib-python$ virtualenv -p python2 venv2
Running virtualenv with interpreter /usr/bin/python2
New python executable in /home/mikm/Projects/pvlib-python/venv2/bin/python2
Also creating executable in /home/mikm/Projects/pvlib-python/venv2/bin/python
Installing setuptools, pkg_resources, pip, wheel...done.
[master] mikm@PTOD702634:~/Projects/pvlib-python$ . venv2/bin/activate
(venv2) [master] mikm@PTOD702634:~/Projects/pvlib-python$ python --version
Python 2.7.12
(venv2) [master] mikm@PTOD702634:~/Projects/pvlib-python$ pip install numpy==1.10.1 pandas==0.15.0 pytz six pytest pytest-cov pytest-mock nose scipy tables numba siphon netcdf4 ephem cython ipython sphinx sphinx_rtd_theme numpydoc matplotlib

tests

(venv2) [analytical_sunrise_sunset_transit] mikm@PTOD702634:~/Projects/pvlib-python$ pytest -v --disable-warnings pvlib/test/test_solarposition.py::test_geometric_sunrise_sunset_transit pvlib/test/test_solarposition.py::test_hour_angle
============================================================================================================ test session starts ============================================================================================================
platform linux2 -- Python 2.7.12, pytest-3.8.2, py-1.6.0, pluggy-0.7.1 -- /home/mikm/Projects/pvlib-python/venv2/bin/python2
cachedir: .pytest_cache
rootdir: /home/mikm/Projects/pvlib-python, inifile:
plugins: cov-2.6.0, mock-1.10.0
collected 2 items

pvlib/test/test_solarposition.py::test_geometric_sunrise_sunset_transit PASSED                                                                                                                                                        [ 50%]
pvlib/test/test_solarposition.py::test_hour_angle PASSED

pip freeze

(venv2) [analytical_sunrise_sunset_transit] mikm@PTOD702634:~/Projects/pvlib-python$ pip freeze
alabaster==0.7.12
atomicwrites==1.2.1
attrs==18.2.0
Babel==2.6.0
backports.functools-lru-cache==1.5
backports.shutil-get-terminal-size==1.0.0
beautifulsoup4==4.6.3
certifi==2018.8.24
cftime==1.0.1
chardet==3.0.4
coverage==4.5.1
cycler==0.10.0
Cython==0.28.5
decorator==4.3.0
docutils==0.14
enum34==1.1.6
ephem==3.7.6.0
funcsigs==1.0.2
idna==2.7
imagesize==1.1.0
ipython==5.8.0
ipython-genutils==0.2.0
Jinja2==2.10
kiwisolver==1.0.1
llvmlite==0.25.0
MarkupSafe==1.0
matplotlib==2.2.3
mock==2.0.0
more-itertools==4.3.0
netCDF4==1.4.1
nose==1.3.7
numba==0.40.0
numexpr==2.6.8
numpy==1.10.1     <-- *** numpy==1.10.1 ***
numpydoc==0.8.0
packaging==18.0
pandas==0.15.0     <-- *** numpy==0.15.0 ***
pathlib2==2.3.2
pbr==4.3.0
pexpect==4.6.0
pickleshare==0.7.5
pkg-resources==0.0.0
pluggy==0.7.1
prompt-toolkit==1.0.15
protobuf==3.6.1
ptyprocess==0.6.0
-e [email protected]:mikofski/pvlib-python.git@bf492cc33dd65653c29bfa20e020e8bce503f799#egg=pvlib
py==1.6.0
Pygments==2.2.0
pyparsing==2.2.2
pytest==3.8.2
pytest-cov==2.6.0
pytest-mock==1.10.0
python-dateutil==2.7.3
pytz==2018.5
requests==2.19.1
scandir==1.9.0
scipy==1.1.0
simplegeneric==0.8.1
singledispatch==3.4.0.3
siphon==0.8.0
six==1.11.0
snowballstemmer==1.2.1
Sphinx==1.8.1
sphinx-rtd-theme==0.4.2
sphinxcontrib-websupport==1.1.0
subprocess32==3.5.2
tables==3.4.4
traitlets==4.3.2
typing==3.6.6
urllib3==1.23
wcwidth==0.1.7

py3x

tests

(py37) C:\Users\mikm\Projects\pvlib-python>pytest -v pvlib\test\test_solarposition.py
============================= test session starts =============================
platform win32 -- Python 3.7.0, pytest-3.8.1, py-1.6.0, pluggy-0.7.1 -- C:\Users\mikm\AppData\Local\Continuum\miniconda3\envs\py37\python.exe
cachedir: .pytest_cache
rootdir: C:\Users\mikm\Projects\pvlib-python, inifile:
collected 35 items

pvlib/test/test_solarposition.py::test_spa_c_physical PASSED             [  2%]
pvlib/test/test_solarposition.py::test_spa_c_physical_dst PASSED         [  5%]
pvlib/test/test_solarposition.py::test_spa_python_numpy_physical PASSED  [  8%]
pvlib/test/test_solarposition.py::test_spa_python_numpy_physical_dst PASSED [ 11%]
pvlib/test/test_solarposition.py::test_spa_python_numba_physical PASSED  [ 14%]
pvlib/test/test_solarposition.py::test_spa_python_numba_physical_dst PASSED [ 17%]
pvlib/test/test_solarposition.py::test_get_sun_rise_set_transit PASSED   [ 20%]
pvlib/test/test_solarposition.py::test_rise_set_transit_ephem PASSED     [ 22%]
pvlib/test/test_solarposition.py::test_rise_set_transit_ephem_error PASSED [ 25%]
pvlib/test/test_solarposition.py::test_rise_set_transit_ephem_horizon PASSED [ 28%]
pvlib/test/test_solarposition.py::test_pyephem_physical PASSED           [ 31%]
pvlib/test/test_solarposition.py::test_pyephem_physical_dst PASSED       [ 34%]
pvlib/test/test_solarposition.py::test_calc_time PASSED                  [ 37%]
pvlib/test/test_solarposition.py::test_earthsun_distance PASSED          [ 40%]
pvlib/test/test_solarposition.py::test_ephemeris_physical PASSED         [ 42%]
pvlib/test/test_solarposition.py::test_ephemeris_physical_dst PASSED     [ 45%]
pvlib/test/test_solarposition.py::test_ephemeris_physical_no_tz PASSED   [ 48%]
pvlib/test/test_solarposition.py::test_get_solarposition_error PASSED    [ 51%]
pvlib/test/test_solarposition.py::test_get_solarposition_pressure[82000-expected0] PASSED [ 54%]
pvlib/test/test_solarposition.py::test_get_solarposition_pressure[90000-expected1] PASSED [ 57%]
pvlib/test/test_solarposition.py::test_get_solarposition_altitude[1830.14-expected0] PASSED [ 60%]
pvlib/test/test_solarposition.py::test_get_solarposition_altitude[2000-expected1] PASSED [ 62%]
pvlib/test/test_solarposition.py::test_get_solarposition_deltat[None-nrel_numpy-expected0] PASSED [ 65%]
pvlib/test/test_solarposition.py::test_get_solarposition_deltat[67.0-nrel_numpy-expected1] PASSED [ 68%]
pvlib/test/test_solarposition.py::test_get_solarposition_deltat[None-nrel_numba-expected2] xfail [ 71%]
pvlib/test/test_solarposition.py::test_get_solarposition_deltat[67.0-nrel_numba-expected3] PASSED [ 74%]
pvlib/test/test_solarposition.py::test_get_solarposition_no_kwargs PASSED [ 77%]
pvlib/test/test_solarposition.py::test_get_solarposition_method_pyephem PASSED [ 80%]
pvlib/test/test_solarposition.py::test_nrel_earthsun_distance PASSED     [ 82%]
pvlib/test/test_solarposition.py::test_equation_of_time PASSED           [ 85%]
pvlib/test/test_solarposition.py::test_declination PASSED                [ 88%]
pvlib/test/test_solarposition.py::test_analytical_zenith PASSED          [ 91%]
pvlib/test/test_solarposition.py::test_analytical_azimuth PASSED         [ 94%]
pvlib/test/test_solarposition.py::test_hour_angle PASSED                 [ 97%]
pvlib/test/test_solarposition.py::test_geometric_sunrise_sunset_transit PASSED [100%]

============================== warnings summary ===============================
<omitted>

============== 34 passed, 1 xfailed, 13 warnings in 8.12 seconds ==============

conda list

(py37) C:\Users\mikm\Projects\pvlib-python>conda list
# packages in environment at C:\Users\mikm\AppData\Local\Continuum\miniconda3\envs\py37:
#
# Name                    Version                   Build  Channel
asn1crypto                0.24.0                   py37_0
astroid                   2.0.4                    py37_0
atomicwrites              1.2.1                    py37_0
attrs                     18.2.0           py37h28b3542_0
backcall                  0.1.0                    py37_0
beautifulsoup4            4.6.3                    py37_0
blas                      1.0                         mkl
bzip2                     1.0.6                hfa6e2cd_5
ca-certificates           2018.03.07                    0
certifi                   2018.8.24                py37_1
cffi                      1.11.5           py37h74b6da3_1
cftime                    1.0.0b1          py37h452e1ab_0
chardet                   3.0.4                    py37_1
colorama                  0.3.9                    py37_0
cryptography              2.3.1            py37h74b6da3_0
curl                      7.61.1               h7602738_0
cycler                    0.10.0                   py37_0
cython                    0.28.5           py37h6538335_0
decorator                 4.3.0                    py37_0
ephem                     3.7.6.0          py37hfa6e2cd_1
freetype                  2.9.1                ha9979f8_1
hdf4                      4.2.13               h712560f_2
hdf5                      1.10.2               hac2f561_1
icc_rt                    2017.0.4             h97af966_0
icu                       58.2                 ha66f8fd_1
idna                      2.7                      py37_0
intel-openmp              2019.0                      118
ipython                   7.0.1            py37h39e3cac_0
ipython_genutils          0.2.0                    py37_0
isort                     4.3.4                    py37_0
jedi                      0.12.1                   py37_0
jinja2                    2.10                     py37_0
jpeg                      9b                   hb83a4c4_2
kiwisolver                1.0.1            py37h6538335_0
lazy-object-proxy         1.3.1            py37hfa6e2cd_2
libcurl                   7.61.1               h7602738_0
libnetcdf                 4.6.1                h70ac11d_1
libpng                    1.6.34               h79bbb47_0
libssh2                   1.8.0                hd619d38_4
llvmlite                  0.25.0                   py37_0
markupsafe                1.0              py37hfa6e2cd_1
matplotlib                2.2.3            py37hd159220_0
mccabe                    0.6.1                    py37_1
mkl                       2019.0                      118
mkl_fft                   1.0.6            py37hdbbee80_0
mkl_random                1.0.1            py37h77b88f5_1
more-itertools            4.3.0                    py37_0
netcdf4                   1.4.1            py37hbfe741f_0
numba                     0.40.0           py37hf9181ef_0
numpy                     1.15.2           py37ha559c80_0
numpy-base                1.15.2           py37h8128ebf_0
openssl                   1.0.2p               hfa6e2cd_0
pandas                    0.23.4           py37h830ac7b_0
parso                     0.3.1                    py37_0
pickleshare               0.7.5                    py37_0
pip                       10.0.1                   py37_0
pluggy                    0.7.1            py37h28b3542_0
prompt_toolkit            2.0.5                    py37_0
pvlib                     0.6.0+10.g2908875.dirty           <pip>
py                        1.6.0                    py37_0
pycparser                 2.19                     py37_0
pygments                  2.2.0                    py37_0
pylint                    2.1.1                    py37_0
pyopenssl                 18.0.0                   py37_0
pyparsing                 2.2.1                    py37_0
pyqt                      5.9.2            py37ha878b3d_0
pysocks                   1.6.8                    py37_0
pytest                    3.8.1                    py37_0
python                    3.7.0                hea74fb7_0
python-dateutil           2.7.3                    py37_0
pytz                      2018.5                   py37_0
qt                        5.9.6            vc14h1e9a669_2
requests                  2.19.1                   py37_0
scipy                     1.1.0            py37h4f6bf74_1
setuptools                40.4.3                   py37_0
simplegeneric             0.8.1                    py37_2
sip                       4.19.12          py37h6538335_0
siphon                    0.8.0                     <pip>
six                       1.11.0                   py37_1
sqlite                    3.25.2               hfa6e2cd_0
tornado                   5.1.1            py37hfa6e2cd_0
traitlets                 4.3.2                    py37_0
urllib3                   1.23                     py37_0
vc                        14.1                 h0510ff6_4
vs2015_runtime            14.15.26706          h3a45250_0
wcwidth                   0.1.7                    py37_0
wheel                     0.32.0                   py37_0
win_inet_pton             1.0.1                    py37_1
wincertstore              0.2                      py37_0
wrapt                     1.10.11          py37hfa6e2cd_2
zlib                      1.2.11               h8395fce_2

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.

Great, thanks for putting in the extra effort for the manual pandas 0.15 testing. In that case I think it's ready to go. I'd like to give @cwhanse a chance to review it again given his detailed comments.

@@ -718,3 +721,45 @@ def test_hour_angle():
# sunrise: 4 seconds, sunset: 48 seconds, transit: 0.2 seconds
# but the differences may be due to other SPA input parameters
assert np.allclose(hours, expected)


def test_geometric_sunrise_sunset_transit(expected_rise_set_spa, golden_mst):
Copy link
Member

Choose a reason for hiding this comment

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

Change function name to test_sunrise_sunset_transit_geometric



def test_geometric_sunrise_sunset_transit(expected_rise_set_spa, golden_mst):
"""Test analytical calculations for sunrise, sunset, and transit times"""
Copy link
Member

Choose a reason for hiding this comment

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

geometric rather than analytical

tz=tz_info)


def _times_to_hours(times):
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 the function name should be more specific, perhaps _hours_after_local_midnight

The function description should also be more specific.



def _local_times_from_hours_since_midnight(times, hours):
"""converts hours from an array of floats to localized times"""
Copy link
Member

Choose a reason for hiding this comment

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

This implies (I think) that hours is hours since local (?) midnight? I think the function description should define what is expected for times and hours. times only serves to communicate a date. hours appears to be hours since midnight.

Parameters
----------
times : pandas.DatetimeIndex
Corresponding timestamps, must be timezone aware.
Copy link
Member

Choose a reason for hiding this comment

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

The description of times needs to specify localized to the timezone for latitude and longitude, it would help to state that sunrise, sunset and transit are calculated for the date part of times.

* Created :py:func:`pvlib.iotools.read_srml` and
:py:func:`pvlib.iotools.read_srml_month_from_solardat` to read University of
Oregon Solar Radiation Monitoring Laboratory data. (:issue:`589`)
* Add geometric functions for sunrise, sunset, and sun transit times,
Copy link
Member

Choose a reason for hiding this comment

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

Wasn't going to pick this nit but since I have other comments, let's sort this list to group the sunrise etc. changes together and the iotools items together.

* _times_to_hours -> _times_to_hours_after_local_midnight
* change corresponding usage in test_solarposition.py and  wrap some
long lines
* in docstring specify that times should be localized to the timzeone
of the lat/lon args in hour_angle() and sunrise_sunset_transit_geometric()
* change arg in solar_azimuth_analytical() and solar_zenith_analytical()
from hour_angle -> hourangle to avoid shadowing func name hour_angle()
from outer scope and remove redundant parentheses from return
* change description to say "geometric" vs. analytical
* improve description of _local_times_from_hours_since_midnight to be
more specific by adding "hours since midnight"
@mikofski
Copy link
Member Author

mikofski commented Oct 9, 2018

@cwhanse thanks for your feedback, I believe I've addressed your comments.

@mikofski
Copy link
Member Author

mikofski commented Oct 9, 2018

FYI: the appveyor error is unrelated, caused by http error

An HTTP error occurred when trying to retrieve this URL.
HTTP errors are often intermittent, and a simple retry will get you on your way.

@cwhanse
Copy link
Member

cwhanse commented Oct 9, 2018

@wholmgren OK to merge?

@wholmgren
Copy link
Member

Thanks Mark!

@wholmgren wholmgren merged commit 4538bef into pvlib:master Oct 9, 2018
@mikofski mikofski deleted the analytical_sunrise_sunset_transit branch November 5, 2018 01:22
@wholmgren wholmgren added this to the 0.6.1 milestone Jan 21, 2019
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