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.

Sign up for GitHub

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
Merged

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

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
19 commits
Select commit Hold shift + click to select a range
3b0e597
add analytical methods from Frank Vignola's book for sunrise, sunset,…
mikofski Sep 19, 2018
af206f2
ENH: remove test data file for sunrise/sunset/transit
mikofski Sep 23, 2018
e6319fe
ENH: PERF: use numpy instead of loop to improve performance
mikofski Sep 24, 2018
59323d3
ENH: use ducktyping to handle different versions of pandas
mikofski Sep 24, 2018
44aed34
ducktype hours if Float64index or not, convert datetime64[D] to ns be…
mikofski Sep 24, 2018
f1e8970
ENH: use "geometric" instead of "analytical"
mikofski Sep 24, 2018
7fa6953
ENH: use more verbose, meaningful names for helpers _hours and _times
mikofski Sep 25, 2018
199a5b5
remove EXPECTED_ANALYTICAL_SUNRISE_SUNSET_TRANSIT
mikofski Sep 30, 2018
a8661f5
stickler - missing whitespace after comma in list
mikofski Sep 30, 2018
5d4012c
use np.asarray to ensure consistent return
mikofski Oct 5, 2018
8b257de
remove try-except, doesn't work anyway, wait for pandas >=0.15.0
mikofski Oct 6, 2018
c15f19b
Merge branch 'master' into analytical_sunrise_sunset_transit
mikofski Oct 6, 2018
a74f944
reuse expected_rise_set_spa test fixture
mikofski Oct 7, 2018
2ace87e
stickler
mikofski Oct 7, 2018
344f679
ENH: update what's new and api.rst
mikofski Oct 8, 2018
bf492cc
stickler - line break before operator
mikofski Oct 9, 2018
887fe29
group sunrise/sunset enhancements together, separate iotools
mikofski Oct 9, 2018
7c90b9c
change private func name to be more descriptive
mikofski Oct 9, 2018
3077923
TST: change name to test_sunrise_sunset_transit_geometric
mikofski Oct 9, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
81 changes: 81 additions & 0 deletions pvlib/solarposition.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1210,3 +1210,84 @@ 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.

"""helper that converts hour angles in degrees into hours"""
tz_info = times.tz # pytz timezone info
tz = tz_info.utcoffset(times).total_seconds() / 3600.
return (hour_angle - longitude - equation_of_time / 4.) / 15. + 12. + tz


def _times(times, hours):
"""helper converts hours from an array of floats to localized times"""
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.

pass
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 •
edited
Loading

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 •
edited
Loading

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 •
edited
Loading

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.

# NEW pandas, "AttributeError: Cannot directly set timezone. Use
# tz_localize() or tz_convert() as appropriate"
times = times.tz_localize(None)
return pd.DatetimeIndex(
times.values.astype('datetime64[D]').astype('datetime64[ns]') +
(hours * 3600. * 1.e9).astype('timedelta64[ns]'),
tz=tz_info)


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 •
edited
Loading

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.


.. 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.


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.

latitude : float
Latitude in degrees
longitude : float
Longitude in degrees
declination : numeric
declination angle in radians at ``times``
equation_of_time : numeric
difference in time between solar time and mean solar time in minutes

Returns
-------
sunrise : datetime
localized sunrise time
sunset : datetime
localized sunset time
transit : datetime
localized sun transit time

References
----------
[2] J. A. Duffie and W. A. Beckman, "Solar Engineering of Thermal
Processes, 3rd Edition" pp. 9-11, J. Wiley and Sons, New York (2006)

[3] Frank Vignola et al., "Solar And Infrared Radiation Measurements",
p. 13, CRC Press (2012)

"""
latitude_rad = np.radians(latitude) # radians
sunset_angle_rad = np.arccos(-np.tan(declination) * np.tan(latitude_rad))
sunset_angle = np.degrees(sunset_angle_rad) # degrees
# solar noon is at hour angle zero
# so sunrise is just negative of sunset
sunrise_angle = -sunset_angle
sunrise_hour = _hours(times, sunrise_angle, longitude, equation_of_time)
sunset_hour = _hours(times, sunset_angle, longitude, equation_of_time)
transit_hour = _hours(times, 0, longitude, equation_of_time)
sunrise = _times(times, sunrise_hour)
sunset = _times(times, sunset_hour)
transit = _times(times, transit_hour)
return sunrise, sunset, transit
58 changes: 57 additions & 1 deletion pvlib/test/test_solarposition.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,6 @@
from conftest import (requires_ephem, needs_pandas_0_17,
requires_spa_c, requires_numba)


# setup times and locations to be tested.
times = pd.date_range(start=datetime.datetime(2014,6,24),
end=datetime.datetime(2014,6,26), freq='15Min')
Expand Down Expand Up @@ -499,3 +498,60 @@ def test_analytical_azimuth():
azimuths = solarposition.solar_azimuth_analytical(*test_angles.T, zenith=zeniths)

assert not np.isnan(azimuths).any()


EXPECTED_ANALYTICAL_SUNRISE_SUNSET_TRANSIT = [
(7.455277777777778, 16.691666666666666, 12.073611111111111),
(7.3975, 16.98138888888889, 12.189444444444444),
(7.161388888888889, 17.338333333333335, 12.25),
(6.79, 17.69472222222222, 12.2425),
(6.307777777777778, 18.0375, 12.1725),
(5.821666666666666, 18.336111111111112, 12.078888888888889),
(5.3613888888888885, 18.629722222222224, 11.995555555555555),
(4.978888888888889, 18.925555555555555, 11.952222222222222),
(4.709722222222222, 19.213333333333335, 11.961666666666666),
(4.617222222222222, 19.405833333333334, 12.011388888888888),
(4.686944444444444, 19.458333333333332, 12.072777777777778),
(4.8825, 19.340555555555557, 12.111666666666666),
(5.158333333333333, 19.043333333333333, 12.100833333333334),
(5.433333333333334, 18.6375, 12.035277777777777),
(5.704166666666667, 18.160833333333333, 11.9325),
(5.9847222222222225, 17.666666666666668, 11.825555555555555),
(6.315277777777778, 17.185, 11.75),
(6.666944444444445, 16.82, 11.743611111111111),
(7.0225, 16.593055555555555, 11.807777777777778),
(7.311111111111111, 16.54277777777778, 11.926944444444445)]


def test_analytical_sunrise_sunset_transit():
"""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

times = pd.DatetimeIndex(start='2018-01-01 0:00:00',
end='2018-12-31 23:59:59',
freq='H').tz_localize('Etc/GMT+7')
lat, lon = 39.743, -105.178 # degrees
eot = solarposition.equation_of_time_pvcdrom(times.dayofyear) # minutes
decl = solarposition.declination_spencer71(times.dayofyear) # radians
sr, ss, st = solarposition.sunrise_sunset_transit_analytical(
times, latitude=lat, longitude=lon, declination=decl,
equation_of_time=eot)
sst = list(
zip(*[(sr[idx], ss[idx], st[idx]) for idx in range(0, 8760, 438)]))
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.
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.
for tr in transit]
test_data_file = EXPECTED_ANALYTICAL_SUNRISE_SUNSET_TRANSIT
test_data_type = np.dtype(
[('sunrise', float), ('sunset', float), ('transit', float)])
test_data = np.array(test_data_file, dtype=test_data_type)
# test data created using NREL SPA-C algorithm at following conditions:
# year=2018, time_zone=-7, longitude=-105.178, latitude=39.743,
# elevation=1830.14, pressure=820, temperature=11, delta_t=67
assert np.allclose(sunrise_hours, test_data['sunrise'])
assert np.allclose(sunset_hours, test_data['sunset'])
assert np.allclose(transit_hours, test_data['transit'])