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

* add analytical methods from Frank Vignola's book for sunrise, sunset, and sun transit times

* ENH: remove test data file for sunrise/sunset/transit

* 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

* ENH: PERF: use numpy instead of loop to improve performance

* add one-liners to explain helpers
* remove import os from test

* ENH: use ducktyping to handle different versions of pandas

* 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

* ducktype hours if Float64index or not, convert datetime64[D] to ns before adding

* ENH: use "geometric" instead of "analytical"

* 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

* ENH: use more verbose, meaningful names for helpers _hours and _times

* change _hours to _hour_angle_to_hours
* change _times to _local_times_from_hours_since_midnite
* fix typos in warning circular orbit
* fix reference numbers
* clarify sign conventions for latitude and longitude

* remove EXPECTED_ANALYTICAL_SUNRISE_SUNSET_TRANSIT

* instead use a test pytest.fixture with a few dates like the ephem and
spa.c test

* stickler - missing whitespace after comma in list

* use np.asarray to ensure consistent return

* remove try-except, doesn't work anyway, wait for pandas >=0.15.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

* reuse expected_rise_set_spa test fixture

* 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

* stickler

* ENH: update what's new and api.rst

* 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

* stickler - line break before operator

* group sunrise/sunset enhancements together, separate iotools

* change private func name to be more descriptive

* _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

* TST: change name to test_sunrise_sunset_transit_geometric

* 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"

Author: mikofski

docs/sphinx/source/api.rst

@@ -79,6 +79,7 @@ calculations.
    solarposition.equation_of_time_spencer71
    solarposition.equation_of_time_pvcdrom
    solarposition.hour_angle
+   solarposition.sunrise_sunset_transit_geometric
 
 
 Clear sky

docs/sphinx/source/whatsnew/v0.6.1.rst

@@ -3,7 +3,7 @@
 v0.6.1 (EST October, 2018)
 --------------------------
 
-This is a minor release. We recommend all users of 0.6.0 upgrade to this
+This is a minor release. We recommend all users of v0.6.0 upgrade to this
 release.
 
 **Python 2.7 support will end on June 1, 2019**. Releases made after this
@@ -16,27 +16,35 @@ API Changes
 * Deprecated ``tmy``, ``tmy.readtmy2`` and ``tmy.readtmy3``;
   they will be removed in v0.7. Use the new :py:func:`pvlib.iotools.read_tmy2`
   and :py:func:`pvlib.iotools.read_tmy3` instead. (:issue:`261`)
-* Added kwarg `horizon` to :func:`~pvlib.solarposition.pyephem` and :func:`~pvlib.solarposition.calc_time` with default value `'+0:00'`
+* Added keyword argument ``horizon`` to :func:`~pvlib.solarposition.pyephem`
+  and :func:`~pvlib.solarposition.calc_time` with default value ``'+0:00'``
 
 
 Enhancements
 ~~~~~~~~~~~~
-* :func:`~pvlib.solarposition.rise_set_transit_ephem` returns sunrise, sunset and transit times using pyephem
-* 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`)
-  
+* :func:`~pvlib.solarposition.rise_set_transit_ephem` returns sunrise, sunset
+  and transit times using pyephem (:issue:`114`)
+* Add geometric functions for sunrise, sunset, and sun transit times,
+  :func:`~pvlib.solarposition.sunrise_sunset_transit_geometric` (:issue:`114`)
+* 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`)
+
 
 Bug fixes
 ~~~~~~~~~
 * Fix when building documentation using Matplotlib 3.0 or greater.
+* Fix and improve :func:`~pvlib.solarposition.hour_angle` (:issue:`598`)
 
 
 Testing
 ~~~~~~~
+* Add test for :func:`~pvlib.solarposition.hour_angle` (:issue:`597`)
 
 
 Contributors
 ~~~~~~~~~~~~
 * Will Holmgren (:ghuser:`wholmgren`)
 * Cliff Hansen (:ghuser:`cwhanse`)
 * Leland Boeman (:ghuser:`lboeman`)
+* Mark Mikofski (:ghuser:`mikofski`)

pvlib/solarposition.py

@@ -25,6 +25,8 @@
 from pvlib import atmosphere
 from pvlib.tools import datetime_to_djd, djd_to_datetime
 
+NS_PER_HR = 1.e9 * 3600.  # nanoseconds per hour
+
 
 def get_solarposition(time, latitude, longitude,
                       altitude=None, pressure=None,
@@ -1175,7 +1177,7 @@ def declination_cooper69(dayofyear):
     return dec
 
 
-def solar_azimuth_analytical(latitude, hour_angle, declination, zenith):
+def solar_azimuth_analytical(latitude, hourangle, declination, zenith):
     """
     Analytical expression of solar azimuth angle based on spherical
     trigonometry.
@@ -1184,7 +1186,7 @@ def solar_azimuth_analytical(latitude, hour_angle, declination, zenith):
     ----------
     latitude : numeric
         Latitude of location in radians.
-    hour_angle : numeric
+    hourangle : numeric
         Hour angle in the local solar time in radians.
     declination : numeric
         Declination of the sun in radians.
@@ -1240,12 +1242,12 @@ def solar_azimuth_analytical(latitude, hour_angle, declination, zenith):
 
     # when NaN values occur in input, ignore and pass to output
     with np.errstate(invalid='ignore'):
-        sign_ha = np.sign(hour_angle)
+        sign_ha = np.sign(hourangle)
 
-    return (sign_ha * np.arccos(cos_azi) + np.pi)
+    return sign_ha * np.arccos(cos_azi) + np.pi
 
 
-def solar_zenith_analytical(latitude, hour_angle, declination):
+def solar_zenith_analytical(latitude, hourangle, declination):
     """
     Analytical expression of solar zenith angle based on spherical
     trigonometry.
@@ -1257,7 +1259,7 @@ def solar_zenith_analytical(latitude, hour_angle, declination):
     ----------
     latitude : numeric
         Latitude of location in radians.
-    hour_angle : numeric
+    hourangle : numeric
         Hour angle in the local solar time in radians.
     declination : numeric
         Declination of the sun in radians.
@@ -1288,7 +1290,7 @@ def solar_zenith_analytical(latitude, hour_angle, declination):
     declination_spencer71 declination_cooper69 hour_angle
     """
     return np.arccos(
-        np.cos(declination) * np.cos(latitude) * np.cos(hour_angle) +
+        np.cos(declination) * np.cos(latitude) * np.cos(hourangle) +
         np.sin(declination) * np.sin(latitude)
     )
 
@@ -1300,7 +1302,8 @@ def hour_angle(times, longitude, equation_of_time):
     Parameters
     ----------
     times : :class:`pandas.DatetimeIndex`
-        Corresponding timestamps, must be timezone aware.
+        Corresponding timestamps, must be localized to the timezone for the
+        ``longitude``.
     longitude : numeric
         Longitude in degrees
     equation_of_time : numeric
@@ -1329,9 +1332,99 @@ def hour_angle(times, longitude, equation_of_time):
     """
     naive_times = times.tz_localize(None)  # naive but still localized
     # hours - timezone = (times - normalized_times) - (naive_times - times)
-    hrs_minus_tzs = 1 / (3600. * 1.e9) * (
+    hrs_minus_tzs = 1 / NS_PER_HR * (
         2 * times.astype(np.int64) - times.normalize().astype(np.int64) -
         naive_times.astype(np.int64))
     # ensure array return instead of a version-dependent pandas <T>Index
     return np.asarray(
         15. * (hrs_minus_tzs - 12.) + longitude + equation_of_time / 4.)
+
+
+def _hour_angle_to_hours(times, hourangle, longitude, equation_of_time):
+    """converts hour angles in degrees to hours as a numpy array"""
+    naive_times = times.tz_localize(None)  # naive but still localized
+    tzs = 1 / NS_PER_HR * (
+        naive_times.astype(np.int64) - times.astype(np.int64))
+    hours = (hourangle - longitude - equation_of_time / 4.) / 15. + 12. + tzs
+    return np.asarray(hours)
+
+
+def _local_times_from_hours_since_midnight(times, hours):
+    """
+    converts hours since midnight from an array of floats to localized times
+    """
+    tz_info = times.tz  # pytz timezone info
+    naive_times = times.tz_localize(None)  # naive but still localized
+    # normalize local, naive times to previous midnight and add the hours until
+    # sunrise, sunset, and transit
+    return pd.DatetimeIndex(
+        (naive_times.normalize().astype(np.int64) +
+         (hours * NS_PER_HR).astype(np.int64)).astype('datetime64[ns]'),
+        tz=tz_info)
+
+
+def _times_to_hours_after_local_midnight(times):
+    """convert local pandas datetime indices to array of hours as floats"""
+    times = times.tz_localize(None)
+    hrs = 1 / NS_PER_HR * (
+        times.astype(np.int64) - times.normalize().astype(np.int64))
+    return np.array(hrs)
+
+
+def sunrise_sunset_transit_geometric(times, latitude, longitude, declination,
+                                     equation_of_time):
+    """
+    Geometric calculation of solar sunrise, sunset, and transit.
+
+    .. warning:: The geometric calculation assumes a circular earth orbit with
+        the sun as a point source at its center, and neglects the effect of
+        atmospheric refraction on zenith. The error depends on location and
+        time of year but is of order 10 minutes.
+
+    Parameters
+    ----------
+    times : pandas.DatetimeIndex
+        Corresponding timestamps, must be localized to the timezone for the
+        ``latitude`` and ``longitude``.
+    latitude : float
+        Latitude in degrees, positive north of equator, negative to south
+    longitude : float
+        Longitude in degrees, positive east of prime meridian, negative to west
+    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
+    ----------
+    [1] J. A. Duffie and W. A. Beckman,  "Solar Engineering of Thermal
+    Processes, 3rd Edition," J. Wiley and Sons, New York (2006)
+
+    [2] Frank Vignola et al., "Solar And Infrared Radiation Measurements,"
+    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 = _hour_angle_to_hours(
+        times, sunrise_angle, longitude, equation_of_time)
+    sunset_hour = _hour_angle_to_hours(
+        times, sunset_angle, longitude, equation_of_time)
+    transit_hour = _hour_angle_to_hours(times, 0, longitude, equation_of_time)
+    sunrise = _local_times_from_hours_since_midnight(times, sunrise_hour)
+    sunset = _local_times_from_hours_since_midnight(times, sunset_hour)
+    transit = _local_times_from_hours_since_midnight(times, transit_hour)
+    return sunrise, sunset, transit

pvlib/test/test_solarposition.py

@@ -529,11 +529,13 @@ def test_get_solarposition_altitude(altitude, expected, golden):
 
 @pytest.mark.parametrize(
     "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')
-    ((None, 'nrel_numba', expected_solpos_multi())),
-    (67.0, 'nrel_numba', expected_solpos_multi())
+        (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')
+        ((None, 'nrel_numba', expected_solpos_multi())),
+        (67.0, 'nrel_numba', expected_solpos_multi())
     ])
 def test_get_solarposition_deltat(delta_t, method, expected, golden):
     times = pd.date_range(datetime.datetime(2003,10,17,13,30,30),
@@ -690,8 +692,9 @@ def test_analytical_azimuth():
                     [ -30., -180.,    5.],
                     [ -30.,  180.,   10.]]))
 
-    zeniths  = solarposition.solar_zenith_analytical(*test_angles.T)
-    azimuths = solarposition.solar_azimuth_analytical(*test_angles.T, zenith=zeniths)
+    zeniths = solarposition.solar_zenith_analytical(*test_angles.T)
+    azimuths = solarposition.solar_azimuth_analytical(*test_angles.T,
+                                                      zenith=zeniths)
 
     assert not np.isnan(azimuths).any()
 
@@ -718,3 +721,48 @@ 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_sunrise_sunset_transit_geometric(expected_rise_set_spa, golden_mst):
+    """Test geometric calculations for sunrise, sunset, and transit times"""
+    times = expected_rise_set_spa.index
+    latitude = golden_mst.latitude
+    longitude = golden_mst.longitude
+    eot = solarposition.equation_of_time_spencer71(times.dayofyear)  # minutes
+    decl = solarposition.declination_spencer71(times.dayofyear)  # radians
+    sr, ss, st = solarposition.sunrise_sunset_transit_geometric(
+        times, latitude=latitude, longitude=longitude, declination=decl,
+        equation_of_time=eot)
+    # sunrise: 2015-01-02 07:26:39.763224487, 2015-08-02 05:04:35.688533801
+    # sunset:  2015-01-02 16:41:29.951096777, 2015-08-02 19:09:46.597355085
+    # transit: 2015-01-02 12:04:04.857160632, 2015-08-02 12:07:11.142944443
+    test_sunrise = solarposition._times_to_hours_after_local_midnight(sr)
+    test_sunset = solarposition._times_to_hours_after_local_midnight(ss)
+    test_transit = solarposition._times_to_hours_after_local_midnight(st)
+    # convert expected SPA sunrise, sunset, transit to local datetime indices
+    expected_sunrise = pd.DatetimeIndex(expected_rise_set_spa.sunrise.values,
+                                        tz='UTC').tz_convert(golden_mst.tz)
+    expected_sunset = pd.DatetimeIndex(expected_rise_set_spa.sunset.values,
+                                       tz='UTC').tz_convert(golden_mst.tz)
+    expected_transit = pd.DatetimeIndex(expected_rise_set_spa.transit.values,
+                                        tz='UTC').tz_convert(golden_mst.tz)
+    # convert expected times to hours since midnight as arrays of floats
+    expected_sunrise = solarposition._times_to_hours_after_local_midnight(
+        expected_sunrise)
+    expected_sunset = solarposition._times_to_hours_after_local_midnight(
+        expected_sunset)
+    expected_transit = solarposition._times_to_hours_after_local_midnight(
+        expected_transit)
+    # geometric time has about 4-6 minute error compared to SPA sunset/sunrise
+    expected_sunrise_error = np.array(
+        [0.07910089555555544, 0.06908014805555496])  # 4.8[min], 4.2[min]
+    expected_sunset_error = np.array(
+        [-0.1036246955555562, -0.06983406805555603])  # -6.2[min], -4.2[min]
+    expected_transit_error = np.array(
+        [-0.011150788888889096, 0.0036508177777765383])  # -40[sec], 13.3[sec]
+    assert np.allclose(test_sunrise, expected_sunrise,
+                       atol=np.abs(expected_sunrise_error).max())
+    assert np.allclose(test_sunset, expected_sunset,
+                       atol=np.abs(expected_sunset_error).max())
+    assert np.allclose(test_transit, expected_transit,
+                       atol=np.abs(expected_transit_error).max())