sunrise and sunset from pyephem

pvlib/solarposition.py

@@ -438,7 +438,26 @@ def get_sun_rise_set_transit(time, latitude, longitude, how='numpy',
     return result
 
 
-def _ephem_setup(latitude, longitude, altitude, pressure, temperature):
+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,
+                 horizon):
     import ephem
     # initialize a PyEphem observer
     obs = ephem.Observer()
@@ -447,14 +466,76 @@ def _ephem_setup(latitude, longitude, altitude, pressure, temperature):
     obs.elevation = altitude
     obs.pressure = pressure / 100.  # convert to mBar
     obs.temp = temperature
+    obs.horizon = horizon
 
     # the PyEphem sun
     sun = ephem.Sun()
     return obs, sun
 
 
+def next_rise_set_ephem(time, latitude, longitude, altitude=0,
+                        pressure=101325, temperature=12, horizon='0:00'):
+    """
+    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.
+    horizon : string, format +/-X:YY
+        arc degrees:arc minutes from geometrical horizon for sunrise and
+        sunset, e.g., horizon='-0:34' when the sun's upper edge crosses the
+        geometrical horizon
+
+    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, horizon)
+    # 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):
+            temperature=12, horizon='+0:00'):
     """
     Calculate the solar position using the PyEphem package.
 
@@ -463,13 +544,19 @@ 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
         air temperature in degrees C.
+    horizon : string, optional, default '+0:00'
+        arc degrees:arc minutes from geometrical horizon for sunrise and
+        sunset, e.g., horizon='-0:34' when the sun's upper edge crosses the
+        geometrical horizon
 
     Returns
     -------
@@ -499,7 +586,7 @@ def pyephem(time, latitude, longitude, altitude=0, pressure=101325,
     sun_coords = pd.DataFrame(index=time)
 
     obs, sun = _ephem_setup(latitude, longitude, altitude,
-                            pressure, temperature)
+                            pressure, temperature, horizon)
 
     # make and fill lists of the sun's altitude and azimuth
     # this is the pressure and temperature corrected apparent alt/az.
@@ -711,7 +798,8 @@ def ephemeris(time, latitude, longitude, pressure=101325, temperature=12):
 
 
 def calc_time(lower_bound, upper_bound, latitude, longitude, attribute, value,
-              altitude=0, pressure=101325, temperature=12, xtol=1.0e-12):
+              altitude=0, pressure=101325, temperature=12, horizon='+0:00',
+              xtol=1.0e-12):
     """
     Calculate the time between lower_bound and upper_bound
     where the attribute is equal to value. Uses PyEphem for
@@ -736,6 +824,10 @@ def calc_time(lower_bound, upper_bound, latitude, longitude, attribute, value,
         atmospheric correction.
     temperature : int or float, optional, default 12
         Air temperature in degrees C.
+    horizon : string, optional, default '+0:00'
+        arc degrees:arc minutes from geometrical horizon for sunrise and
+        sunset, e.g., horizon='-0:34' when the sun's upper edge crosses the
+        geometrical horizon
     xtol : float, optional, default 1.0e-12
         The allowed error in the result from value
 
@@ -758,7 +850,7 @@ def calc_time(lower_bound, upper_bound, latitude, longitude, attribute, value,
         raise ImportError('The calc_time function requires scipy')
 
     obs, sun = _ephem_setup(latitude, longitude, altitude,
-                            pressure, temperature)
+                            pressure, temperature, horizon)
 
     def compute_attr(thetime, target, attr):
         obs.date = thetime

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,51 @@ 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_spa():
+    # for Golden, CO, from NREL SPA 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, 21, 55),
+                                datetime.datetime(2015, 8, 2, 5, 0, 27)
+                                ]).tz_localize('MST').tolist()
+    sunset = pd.DatetimeIndex([datetime.datetime(2015, 1, 2, 16, 47, 43),
+                               datetime.datetime(2015, 8, 2, 19, 13, 58)
+                               ]).tz_localize('MST').tolist()
+    transit = pd.DatetimeIndex([datetime.datetime(2015, 1, 2, 12, 4, 45),
+                                datetime.datetime(2015, 8, 2, 12, 6, 58)
+                                ]).tz_localize('MST').tolist()
+    return pd.DataFrame({'transit': transit,
+                         'sunrise': sunrise,
+                         'sunset': sunset},
+                        index=times)
+
+
+@pytest.fixture()
+def expected_rise_set_ephem():
+    # for Golden, CO, from USNO websites
+    times = pd.DatetimeIndex([datetime.datetime(2015, 1, 2),
+                              datetime.datetime(2015, 1, 3),
+                              datetime.datetime(2015, 8, 2),
+                              ]).tz_localize('MST')
+    sunrise = pd.DatetimeIndex([datetime.datetime(2015, 1, 2, 7, 22, 0),
+                                datetime.datetime(2015, 1, 3, 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, 1, 3, 16, 49, 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, 1, 3, 12, 5, 0),
+                                datetime.datetime(2015, 8, 2, 12, 7, 0)
+                                ]).tz_localize('MST').tolist()
+    return pd.DataFrame({'transit': transit,
+                         'sunrise': sunrise,
+                         'sunset': sunset},
+                        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 +172,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_spa):
     south = Location(-35.0, 0.0, tz='UTC')
     times = pd.DatetimeIndex([datetime.datetime(1996, 7, 5, 0),
                               datetime.datetime(2004, 12, 4, 0)]
@@ -134,49 +181,103 @@ def test_get_sun_rise_set_transit():
                                 datetime.datetime(2004, 12, 4, 4, 38, 57)]
                                ).tz_localize('UTC').tolist()
     sunset = pd.DatetimeIndex([datetime.datetime(1996, 7, 5, 17, 1, 4),
-                               datetime.datetime(2004, 12, 4, 19, 2, 2)]
+                               datetime.datetime(2004, 12, 4, 19, 2, 3)]
                               ).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)
+                                                    delta_t=65.0)
     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)
 
+    # test for Golden, CO compare to NREL SPA
+    result = solarposition.get_sun_rise_set_transit(
+            expected_rise_set_spa.index, golden.latitude, golden.longitude,
+            delta_t=65.0)
 
-    # 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,
-                                                    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('s').tz_convert('MST')
 
-    del result_rounded['transit']
-    assert_frame_equal(frame, result_rounded)
+    assert_frame_equal(expected_rise_set_spa, result_rounded)
+
+
+@requires_ephem
+def test_next_rise_set_ephem(expected_rise_set_ephem):
+    # test for Golden, CO compare to USNO
+    result = solarposition.next_rise_set_ephem(expected_rise_set_ephem.index,
+                                               golden.latitude,
+                                               golden.longitude,
+                                               golden.altitude,
+                                               pressure=0,
+                                               temperature=11,
+                                               horizon='-0:34')
+    # 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_ephem['transit']
+    assert_frame_equal(expected_rise_set_ephem, result_rounded)
+
+    # test with times
+    times = pd.DatetimeIndex([datetime.datetime(2015, 1, 2, 3, 0, 0),
+                              datetime.datetime(2015, 1, 2, 10, 15, 0),
+                              datetime.datetime(2015, 1, 2, 15, 3, 0),
+                              datetime.datetime(2015, 1, 2, 21, 6, 7)
+                              ]).tz_localize('MST')
+    expected = pd.DataFrame(index=times,
+                            columns=['sunrise', 'sunset'],
+                            dtype='datetime64[ns]')
+    expected['sunrise'] = pd.Series(index=times, data =
+       [expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 2), 'sunrise'],
+        expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 3), 'sunrise'],
+        expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 3), 'sunrise'],
+        expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 3), 'sunrise']])
+    expected['sunset'] = pd.Series(index=times, data =
+       [expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 2), 'sunset'],
+        expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 2), 'sunset'],
+        expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 2), 'sunset'],
+        expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 3), 'sunset']])
+
+    result = solarposition.next_rise_set_ephem(times,
+                                               golden.latitude,
+                                               golden.longitude,
+                                               golden.altitude,
+                                               pressure=0,
+                                               temperature=11,
+                                               horizon='-0:34')
+    # 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')
+    assert_frame_equal(expected, result_rounded)
+
+    # test with different timezone
+    times = times.tz_convert('UTC')
+    expected = expected.tz_convert('UTC')
+    for col, data in expected.iteritems():
+        expected[col] = data.dt.tz_convert('UTC')
+    result = solarposition.next_rise_set_ephem(times,
+                                               golden.latitude,
+                                               golden.longitude,
+                                               golden.altitude,
+                                               pressure=0,
+                                               temperature=11,
+                                               horizon='-0:34')
+    # 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(times.tz)
+    assert_frame_equal(expected, result_rounded)
 
 
 @requires_ephem
@@ -190,6 +291,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 +303,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 +330,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),