finding twilight times with skyfield

Question

I have a function that returns the sunrise, sunset, solar-noon and twilight times for a given location. This function uses pyephem and since it's deprecated, I'd figure I'd re-write the function to use skyfield. However, skyfield doesn't have the functions .previous_rising, .next_setting or .next_transit (at least, I can't find them), which is what I used with pyephem.

skyfield does have the function find_discrete, which will search between given times to find when a function changes, so I wrote the following to test:

from datetime import datetime, timedelta
from skyfield import api, almanac
import pytz

ts = api.load.timescale()
planets = api.load('de421.bsp')
sun, earth = planets['sun'], planets['earth']

lat, lon =  '44.59028 N', '104.71528 W'  # UFO Mooring Site
tzn, elv = 'US/Mountain', 1559

tz = pytz.timezone(tzn)
loc = api.Topos(lat, lon, elevation_m=elv)
earth_loc = earth + loc

def civil_twil(time):
    "returns true/false if sun is above -6 degrees"
    alt, az, dis = earth_loc.at(time).observe(sun).apparent().altaz()
    return alt > -6

t0 = ts.utc(datetime.now(tz))
t1 = ts.utc(tz.normalize(datetime.now(tz) + timedelta(1)))
civil_dawn = almanac.find_discrete(t0, t1, civil_twil)

But this just gives me an error that the function is missing the 'rough_period' attribute and the documentation doesn't mention what that might be. I'm going to guess that maybe it's only valid for the functions defined in the almanac class; but again, it's not mentioned.

So, how can I find/determine the twilight times via skyfield?

rfkortekaas · Accepted Answer · 2019-02-03T15:07:33.973

There is an example in the documentation to find the value for the sunrise and sunset with the find_discrete function. See the source in almanac.py for this function. The only problem here is that it calculates only when the top of the sun is apparently even with the horizon. In the following example I've changed the sunrise_sunset function to a daylength function. With this you can give the required angle to the daylength function.

from skyfield import api, almanac
from datetime import datetime, timedelta
import pytz
from skyfield.nutationlib import iau2000b

DAYLENGTH_CENTER_HORIZON = 0.0
DAYLENGTH_TOP_HORIZON = 0.26667
DAYLENGTH_TOP_HORIZON_APPARENTLY = 0.8333
DAYLENGTH_CIVIL_TWILIGHT = 6.0
DAYLENGTH_NAUTICAL_TWILIGHT = 12.0
DAYLENGTH_ASTRONOMICAL_TWILIGHT = 18.0

def daylength(ephemeris, topos, degrees):
    """Build a function of time that returns the daylength.

    The function that this returns will expect a single argument that is a 
    :class:`~skyfield.timelib.Time` and will return ``True`` if the sun is up
    or twilight has started, else ``False``.
    """
    sun = ephemeris['sun']
    topos_at = (ephemeris['earth'] + topos).at

    def is_sun_up_at(t):
        """Return `True` if the sun has risen by time `t`."""
        t._nutation_angles = iau2000b(t.tt)
        return topos_at(t).observe(sun).apparent().altaz()[0].degrees > -degrees

    is_sun_up_at.rough_period = 0.5  # twice a day
    return is_sun_up_at

ts = api.load.timescale()
planets = api.load('de421.bsp')
sun = planets['sun']
earth = planets['earth']

lat, lon = '44.59028 N', '104.71528 W'  # UFO Mooring Site
tzn, elv = 'US/Mountain', 1559

tz = pytz.timezone(tzn)
loc = api.Topos(lat, lon, elevation_m=elv)

t0 = ts.utc(datetime.now(tz))
t1 = ts.utc(tz.normalize(datetime.now(tz) + timedelta(1)))
center_time, center_up = almanac.find_discrete(t0, t1, daylength(planets, loc,
                                                    DAYLENGTH_CENTER_HORIZON))
print('Sunrise Sunset center of sun is even with horizon:')
print(center_time.utc_iso(), center_up)

apparent_top_time, apparent_top_up = almanac.find_discrete(t0, t1,
        daylength(planets, loc, DAYLENGTH_TOP_HORIZON_APPARENTLY))
print('Sunrise Sunset top of sun is apparently even with horizon:')
print(apparent_top_time.utc_iso(), apparent_top_up)


civil_twil_time, civil_twil_up = almanac.find_discrete(t0, t1,
        daylength(planets, loc, DAYLENGTH_CIVIL_TWILIGHT))
print('Civil twilight:')
print(civil_twil_time.utc_iso(), civil_twil_up)

This will print the following results:

Sunrise Sunset center of sun is even with horizon:
['2019-02-03T14:20:33Z', '2019-02-04T00:05:20Z'] [ True False]

Sunrise Sunset top of sun is apparently even with horizon:
['2019-02-03T14:15:28Z', '2019-02-04T00:10:25Z'] [ True False]

Civil twilight:
['2019-02-03T13:44:36Z', '2019-02-04T00:41:18Z'] [ True False]

The first list shows the times when a change is found and the second list shows True if the sun is up (or twilight started) and False when the sun is set.

The rough_period that was missing in your case should be a float value with how much it happens a day. The sun rises and sets once a day so the event in this function happens twice a day. This means the rough_period is 0.5. For example when you want to calculate the moon phases the rough_period can be set to 7.0 (a full moon orbit is 27.3 days which results in 6.825 days per phase). See the other examples for the season or moon phase calculation in the almanac.py source code.

finding twilight times with skyfield

1 Answers1