ENH: geometric calculations for sunrise, sunset, and transit

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`)
+* 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,
+  :func:`~pvlib.solarposition.sunrise_sunset_transit_geometric` (:issue:`114`)
+
 
 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,
@@ -1329,9 +1331,96 @@ 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 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(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 timezone aware.
+    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,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):
+    """Test analytical 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(sr)
+    test_sunset = solarposition._times_to_hours(ss)
+    test_transit = solarposition._times_to_hours(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(expected_sunrise)
+    expected_sunset = solarposition._times_to_hours(expected_sunset)
+    expected_transit = solarposition._times_to_hours(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())