ENH: geometric calculations for sunrise, sunset, and transit #583
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -1210,3 +1210,69 @@ 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): | ||
| 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): | ||
| 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)]) | ||
|
There was a problem hiding this comment. This loop using pure python is likely to be painfully slow and defeat the purpose of an analytical approach. There was a problem hiding this comment. I think I can use numpy arrays with There was a problem hiding this comment. Nice catch @wholmgren , I get about a 70x speed improvement (from 43.5 ms to 605us) to by using numpy 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) |
||
|
|
||
cwhanse marked this conversation as resolved.
Show resolved
Hide resolved
|
||
|
|
||
| def sunrise_sunset_transit_analytical(times, latitude, longitude, declination, | ||
| equation_of_time): | ||
| """ | ||
| Analytical calculation of solar sunrise, sunset, and transit. | ||
|
There was a problem hiding this comment. I don't care for the term There was a problem hiding this comment. 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., There was a problem hiding this comment. "geometric" might be an alternative to "analytical" There was a problem hiding this comment. I really like "geometric" but would we change zenith and azimuth to match? There was a problem hiding this comment. I'd say yes, to follow the naming pattern, tacitly labeling the collection of spherical trigonometry methods as the 'geometric' model. Deprecate There was a problem hiding this comment. 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! There was a problem hiding this comment. 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. | ||
|
There was a problem hiding this comment. Probably "geometric" instead of analytic. Also guessing it neglects more than refraction, but haven't looked at the sources. There was a problem hiding this comment. Also spherical earth, sun is represented by a point at its center, ... There was a problem hiding this comment. 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." There was a problem hiding this comment. from the old test the max differences are
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. | ||
|
There was a problem hiding this comment. The description of |
||
| latitude : float | ||
| Latitude in degrees | ||
| longitude : float | ||
cwhanse marked this conversation as resolved.
Show resolved
Hide resolved
|
||
| 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) | ||
|
|
||
| """ | ||
cwhanse marked this conversation as resolved.
Show resolved
Hide resolved
cwhanse marked this conversation as resolved.
Show resolved
Hide resolved
|
||
| 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 | ||
There was a problem hiding this comment.
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
_hoursand_times.