How to compute sunrise & sunset in Python using skyfield

The following code will compute the sunrise & sunset at a specific location & elevation using the skyfield library. Note that you need to download de413.bsp .

from skyfield import api
from skyfield import almanac
from datetime import datetime
from datetime import timedelta
import dateutil.parser
from calendar import monthrange

ts = api.load.timescale()
ephem = api.load_file('de413.bsp')

def compute_sunrise_sunset(location, year=2019, month=1, day=1):
    t0 = ts.utc(year, month, day, 0)
    # t1 = t0 plus one day
    t1 = ts.utc(t0.utc_datetime() + timedelta(days=1))
    t, y = almanac.find_discrete(t0, t1, almanac.sunrise_sunset(ephem, location))
    sunrise = None
    for time, is_sunrise in zip(t, y):
        if is_sunrise:
            sunrise = dateutil.parser.parse(time.utc_iso())
        else:
            sunset = dateutil.parser.parse(time.utc_iso())
    return sunrise, sunset

# Compute sunrise & sunset for random location near Munich
location = api.Topos('48.324777 N', '11.405610 E', elevation_m=519)
now = datetime.now()
sunrise, sunset = compute_sunrise_sunset(location, now.year, now.month, now.day)

# Print result (example)
print(f'Sunrise today: {sunrise}')
print(f'Sunset today: {sunset}')

Definition of sunrise & sunset in this context

According to the skyfield documentation:

Skyfield uses the same definition as the United States Naval Observatory: the Sun is up when its center is 0.8333 degrees below the horizon, which accounts for both its apparent radius of around 16 arcminutes and also for the 34 arcminutes by which atmospheric refraction on average lifts the image of the Sun.

Other caveats

Example output:

Sunrise today: 2022-06-19 03:12:56+00:00
Sunset today: 2022-06-19 19:18:38+00:00