sunrise and sunset from pyephem

pvlib/solarposition.py

@@ -438,6 +438,24 @@ def get_sun_rise_set_transit(time, latitude, longitude, how='numpy',
     return result
 
 
+def _ephem_convert_to_seconds_and_microseconds(date):
+    # utility from unreleased PyEphem 3.6.7.1
+    """Converts a PyEphem date into seconds"""
+    microseconds = int(round(24 * 60 * 60 * 1000000 * date))
+    seconds, microseconds = divmod(microseconds, 1000000)
+    seconds -= 2209032000  # difference between epoch 1900 and epoch 1970
+    return seconds, microseconds
+
+
+def _ephem_to_timezone(date, tzinfo):
+    # utility from unreleased PyEphem 3.6.7.1
+    """"Convert a PyEphem Date into a timezone aware python datetime"""
+    seconds, microseconds = _ephem_convert_to_seconds_and_microseconds(date)
+    date = dt.datetime.fromtimestamp(seconds, tzinfo)
+    date = date.replace(microsecond=microseconds)
+    return date
+
+
 def _ephem_setup(latitude, longitude, altitude, pressure, temperature):
     import ephem
     # initialize a PyEphem observer
@@ -453,6 +471,63 @@ def _ephem_setup(latitude, longitude, altitude, pressure, temperature):
     return obs, sun
 
 
+def ephem_next_rise_set(time, latitude, longitude, altitude=0,
+                        pressure=101325, temperature=12):
+    """
+    Calculate the next sunrise and sunset times using the PyEphem package.
+
+    Parameters
+    ----------
+    time : pandas.DatetimeIndex
+        Localized or UTC.
+    latitude : float
+        positive is north of 0
+    longitude : float
+        positive is east of 0
+    altitude : float, default 0
+        distance above sea level in meters.
+    pressure : int or float, optional, default 101325
+        air pressure in Pascals.
+    temperature : int or float, optional, default 12
+        air temperature in degrees C.
+
+    Returns
+    -------
+    pandas.DataFrame
+        index is the same as input time.index
+        columns are 'sunrise' and 'sunset', times are localized to the
+        timezone time.tz
+
+    See also
+    --------
+    pyephem
+    """
+
+    try:
+        import ephem
+    except ImportError:
+        raise ImportError('PyEphem must be installed')
+
+    # if localized, convert to UTC. otherwise, assume UTC.
+    try:
+        time_utc = time.tz_convert('UTC')
+    except TypeError:
+        time_utc = time
+
+    obs, sun = _ephem_setup(latitude, longitude, altitude,
+                            pressure, temperature)
+    # create lists of the next sunrise and sunset time localized to time.tz
+    next_sunrise = []
+    next_sunset = []
+    for thetime in time_utc:
+        obs.date = ephem.Date(thetime)
+        next_sunrise.append(_ephem_to_timezone(obs.next_rising(sun), time.tz))
+        next_sunset.append(_ephem_to_timezone(obs.next_setting(sun), time.tz))
+
+    return pd.DataFrame(index=time, data={'sunrise': next_sunrise,
+                                          'sunset': next_sunset})
+
+
 def pyephem(time, latitude, longitude, altitude=0, pressure=101325,
             temperature=12):
     """
@@ -463,9 +538,11 @@ def pyephem(time, latitude, longitude, altitude=0, pressure=101325,
     time : pandas.DatetimeIndex
         Localized or UTC.
     latitude : float
+        positive is north of 0
     longitude : float
+        positive is east of 0
     altitude : float, default 0
-        distance above sea level.
+        distance above sea level in meters.
     pressure : int or float, optional, default 101325
         air pressure in Pascals.
     temperature : int or float, optional, default 12

pvlib/test/test_solarposition.py

@@ -29,6 +29,7 @@
 
 tol = 5
 
+
 @pytest.fixture()
 def expected_solpos():
     return pd.DataFrame({'elevation': 39.872046,
@@ -37,6 +38,7 @@ def expected_solpos():
                          'apparent_elevation': 39.888378},
                         index=['2003-10-17T12:30:30Z'])
 
+
 @pytest.fixture()
 def expected_solpos_multi():
     return pd.DataFrame({'elevation': [39.872046, 39.505196],
@@ -45,6 +47,27 @@ def expected_solpos_multi():
                          'apparent_elevation': [39.888378, 39.521740]},
                         index=[['2003-10-17T12:30:30Z', '2003-10-18T12:30:30Z']])
 
+
+@pytest.fixture()
+def expected_rise_set():
+    # for Golden, CO, from USNO website
+    times = pd.DatetimeIndex([datetime.datetime(2015, 1, 2),
+                              datetime.datetime(2015, 8, 2),
+                              ]).tz_localize('MST')
+    sunrise = pd.DatetimeIndex([datetime.datetime(2015, 1, 2, 7, 22, 0),
+                                datetime.datetime(2015, 8, 2, 5, 0, 0)
+                                ]).tz_localize('MST').tolist()
+    sunset = pd.DatetimeIndex([datetime.datetime(2015, 1, 2, 16, 48, 0),
+                               datetime.datetime(2015, 8, 2, 19, 13, 0)
+                               ]).tz_localize('MST').tolist()
+    transit = pd.DatetimeIndex([datetime.datetime(2015, 1, 2, 12, 5, 0),
+                                datetime.datetime(2015, 8, 2, 12, 7, 0)
+                               ]).tz_localize('MST').tolist()
+    return pd.DataFrame({'sunrise': sunrise,
+                         'sunset': sunset,
+                         'transit': transit}, index=times)
+
+
 # the physical tests are run at the same time as the NREL SPA test.
 # pyephem reproduces the NREL result to 2 decimal places.
 # this doesn't mean that one code is better than the other.
@@ -125,7 +148,7 @@ def test_spa_python_numba_physical_dst(expected_solpos):
 
 
 @needs_pandas_0_17
-def test_get_sun_rise_set_transit():
+def test_get_sun_rise_set_transit(expected_rise_set):
     south = Location(-35.0, 0.0, tz='UTC')
     times = pd.DatetimeIndex([datetime.datetime(1996, 7, 5, 0),
                               datetime.datetime(2004, 12, 4, 0)]
@@ -136,47 +159,50 @@ def test_get_sun_rise_set_transit():
     sunset = pd.DatetimeIndex([datetime.datetime(1996, 7, 5, 17, 1, 4),
                                datetime.datetime(2004, 12, 4, 19, 2, 2)]
                               ).tz_localize('UTC').tolist()
+    frame = pd.DataFrame({'sunrise':sunrise, 'sunset':sunset}, index=times)
+
     result = solarposition.get_sun_rise_set_transit(times, south.latitude,
                                                     south.longitude,
                                                     delta_t=64.0)
-    frame = pd.DataFrame({'sunrise':sunrise, 'sunset':sunset}, index=times)
     result_rounded = pd.DataFrame(index=result.index)
     # need to iterate because to_datetime does not accept 2D data
     # the rounding fails on pandas < 0.17
     for col, data in result.iteritems():
-        result_rounded[col] = pd.to_datetime(
-            np.floor(data.values.astype(np.int64) / 1e9)*1e9, utc=True)
+        result_rounded[col] = data.dt.round('1s')
 
     del result_rounded['transit']
     assert_frame_equal(frame, result_rounded)
 
-
-    # tests from USNO
-    # Golden
-    golden = Location(39.0, -105.0, tz='MST')
-    times = pd.DatetimeIndex([datetime.datetime(2015, 1, 2),
-                              datetime.datetime(2015, 8, 2),]
-                             ).tz_localize('MST')
-    sunrise = pd.DatetimeIndex([datetime.datetime(2015, 1, 2, 7, 19, 2),
-                                datetime.datetime(2015, 8, 2, 5, 1, 26)
-                                ]).tz_localize('MST').tolist()
-    sunset = pd.DatetimeIndex([datetime.datetime(2015, 1, 2, 16, 49, 10),
-                               datetime.datetime(2015, 8, 2, 19, 11, 31)
-                               ]).tz_localize('MST').tolist()
-    result = solarposition.get_sun_rise_set_transit(times, golden.latitude,
+    # test for Golden, CO compare to USNO
+    result = solarposition.get_sun_rise_set_transit(expected_rise_set.index,
+                                                    golden.latitude,
                                                     golden.longitude,
                                                     delta_t=64.0)
-    frame = pd.DataFrame({'sunrise':sunrise, 'sunset':sunset}, index=times)
+
+    # round to nearest minute
     result_rounded = pd.DataFrame(index=result.index)
     # need to iterate because to_datetime does not accept 2D data
-    # the rounding fails on pandas < 0.17
     for col, data in result.iteritems():
-        result_rounded[col] = (pd.to_datetime(
-            np.floor(data.values.astype(np.int64) / 1e9)*1e9, utc=True)
-            .tz_convert('MST'))
+        result_rounded[col] = data.dt.round('min').tz_convert('MST')
 
-    del result_rounded['transit']
-    assert_frame_equal(frame, result_rounded)
+    assert_frame_equal(expected_rise_set, result_rounded)
+
+
+@requires_ephem
+def test_next_rise_set_ephem(expected_rise_set):
+    # test for Golden, CO compare to USNO
+    result = solarposition.ephem_next_rise_set(expected_rise_set.index,
+                                               golden.latitude,
+                                               golden.longitude,
+                                               golden.altitude,
+                                               pressure=101325,
+                                               temperature=12)
+    # round to nearest minute
+    result_rounded = pd.DataFrame(index=result.index)
+    for col, data in result.iteritems():
+        result_rounded[col] = data.dt.round('min').tz_convert('MST')
+    del expected_rise_set('transit')
+    assert_frame_equal(expected_rise_set, result_rounded)
 
 
 @requires_ephem
@@ -190,6 +216,7 @@ def test_pyephem_physical(expected_solpos):
     assert_frame_equal(expected_solpos.round(2),
                        ephem_data[expected_solpos.columns].round(2))
 
+
 @requires_ephem
 def test_pyephem_physical_dst(expected_solpos):
     times = pd.date_range(datetime.datetime(2003,10,17,13,30,30), periods=1,
@@ -201,6 +228,7 @@ def test_pyephem_physical_dst(expected_solpos):
     assert_frame_equal(expected_solpos.round(2),
                        ephem_data[expected_solpos.columns].round(2))
 
+
 @requires_ephem
 def test_calc_time():
     import pytz
@@ -227,6 +255,7 @@ def test_calc_time():
     assert_allclose((az.replace(second=0, microsecond=0) -
                           epoch_dt).total_seconds(), actual_timestamp)
 
+
 @requires_ephem
 def test_earthsun_distance():
     times = pd.date_range(datetime.datetime(2003,10,17,13,30,30),