ENH: geometric calculations for sunrise, sunset, and transit

pvlib/solarposition.py

@@ -1335,3 +1335,89 @@ def hour_angle(times, longitude, equation_of_time):
     # 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, hour_angle, 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 / (3600. * 1.e9) * (
+            naive_times.astype(np.int64) - times.astype(np.int64))
+    hours = (hour_angle - longitude - equation_of_time / 4.) / 15. + 12. + tzs
+    return np.asarray(hours)
+
+
+def _local_times_from_hours_since_midnite(times, hours):
+    """converts hours from an array of floats to localized times"""
+    tz_info = times.tz  # pytz timezone info
+    times = times.tz_localize(None)  # naive but still localized
+    # convert times to previous naive, local midnight
+    times = times.values.astype('datetime64[D]').astype('datetime64[ns]')
+    # add the hours until sunrise/sunset/transit
+    return pd.DatetimeIndex(
+        times + (hours * 3600. * 1.e9).astype('timedelta64[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 / (3600. * 1.e9) * (
+        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_midnite(times, sunrise_hour)
+    sunset = _local_times_from_hours_since_midnite(times, sunset_hour)
+    transit = _local_times_from_hours_since_midnite(times, transit_hour)
+    return sunrise, sunset, transit

pvlib/test/test_solarposition.py

@@ -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(
+        raises=ValueError,
+        reason = 'spa.calculate_deltat not implemented for numba yet')
     ((None, 'nrel_numba', expected_solpos_multi())),
     (67.0, 'nrel_numba', expected_solpos_multi())
     ])
@@ -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())