Sunrise and Sunset Times¶
This tutorial computes sunrise and sunset times for every day of 2024 at a given location, then visualizes how daylight hours change over the course of the year.
satkit.sun.rise_set returns the sunrise and sunset times in UTC for a given date and ground location, using Algorithm 30 of Vallado (2013), §5.3.1, with the standard almanac definitions of sunrise (upper limb at 90°50′ from the zenith) and civil, nautical and astronomical twilight (96°, 102°, 108°). The results are converted to local time (US Eastern) for display. The input time selects the UTC calendar date (its time of day is ignored), and the sunrise and sunset returned are those of that date at the location's longitude. To get the events of a local date, pass a timezone-aware time at local noon on that date: local noon falls on the same UTC date for time zones UTC−11 to UTC+11, while local midnight falls on the previous UTC date east of Greenwich. The shaded region in the plot represents daylight hours.
The Sun's position is re-evaluated at each computed event until it converges, and it is the apparent Sun of sun.pos_mod (Meeus's solar coordinates with the largest planetary and lunar terms, good to 3.6″ in longitude) with nutation and solar parallax, so the times agree with Skyfield (JPL DE421) to within 0.5 s, typically 0.2 s, at latitudes up to 65° (a single pass of the book algorithm is off by up to 35 s there). UTC stands in for UT1 (they differ by less than 0.9 s). With use_jpl=True the Sun comes from the JPL ephemeris, with the full Earth-orientation reduction, and the times agree with Skyfield to 0.01 s. The horizon is at sea level: the observer's altitude is ignored. An observer h meters above the surrounding terrain or sea sees the horizon lowered by the dip, about 1.76′ × √h, which makes sunrise earlier and sunset later; pass sigma=90 + (50 + 1.76 * math.sqrt(h)) / 60 to account for it. On days of polar day or polar night rise_set raises an exception.
import satkit as sk
import datetime
from zoneinfo import ZoneInfo
import matplotlib.pyplot as plt
import scienceplots # noqa: F401
plt.style.use(["science", "no-latex", "../satkit.mplstyle"])
%config InlineBackend.figure_formats = ['svg']
Compute Rise and Set Times¶
For each day in 2024, compute the sunrise and sunset in UTC, then convert to US Eastern time using Python's built-in zoneinfo module.
# Array of times for each day of 2024
basetime = sk.time(2024, 1, 1)
timearr = [basetime + sk.duration(days=i) for i in range(365)]
# Coordinates of Arlington, MA
coord = sk.itrfcoord(latitude_deg=42.1514, longitude_deg=-71.1516)
# sunrise, sunset in UTC
riseset = [sk.sun.rise_set(t, coord) for t in timearr]
rise, set = zip(*riseset)
# Convert to Eastern Time
eastern = ZoneInfo("America/New_York")
drise = [r.to_datetime().astimezone(eastern) for r in rise]
dset = [s.to_datetime().astimezone(eastern) for s in set]
# Hour of day, in [0,24]
risefrac = [r.hour + r.minute / 60 + r.second / 3600 for r in drise]
setfrac = [s.hour + s.minute / 60 + s.second / 3600 for s in dset]
# Convert hour of day to a time
risetime = [datetime.time(hour=r.hour, minute=r.minute, second=r.second) for r in drise]
settime = [datetime.time(hour=s.hour, minute=s.minute, second=s.second) for s in dset]
Visualize¶
Plot sunrise and sunset times over the year. The shaded area represents daylight hours. Note the characteristic asymmetry around the solstices and the abrupt shifts from daylight saving time transitions.
def frac2str(y):
return y.strftime("%H:%M:%S")
risestring = [frac2str(r) for r in risetime]
setstring = [frac2str(s) for s in settime]
fig, ax = plt.subplots(figsize=(8, 5))
dates = [x.to_datetime() for x in timearr]
ax.plot(dates, risefrac, color='#006450', label='Sunrise')
ax.plot(dates, setfrac, color='#006450', label='Sunset')
ax.fill_between(dates, risefrac, setfrac, alpha=0.2, color='#006450')
ax.set_xlabel("Date")
ax.set_ylabel("Local Hour of Day")
ax.set_title("Sunrise and Sunset Times for 2024 in Arlington, MA")
ax.set_ylim(0, 24)
ax.set_yticks([0, 6, 12, 18, 24])
ax.set_yticklabels(["Midnight", "6 AM", "Noon", "6 PM", "Midnight"])
ax.legend()
fig.autofmt_xdate()
plt.tight_layout()
plt.show()
Moonrise and Moonset¶
satkit.moon.rise_set returns the moonrise and moonset in UTC with the same day convention as sun.rise_set: the input's UTC calendar date, over the local mean day at the location's longitude. The events are when the Moon's upper limb touches the horizon with 34′ of refraction, the USNO / Astronomical Almanac definition, with the Moon's parallax (up to 61′) and changing semidiameter computed exactly. Because the Moon rises about 50 minutes later each day, about once a month a date has no moonrise, and another has no moonset; those events are None.
The built-in analytic Moon is good to a minute or two at mid latitudes; use_jpl=True uses the JPL ephemeris instead, to within a second of other JPL-based almanacs.
# Moonrise and moonset, in UTC, for each day of 2024 in Arlington, MA
moon_riseset = [sk.moon.rise_set(t, coord) for t in timearr]
def local_hour(t):
d = t.to_datetime().astimezone(eastern)
return d.hour + d.minute / 60 + d.second / 3600
mrise = [(t.to_datetime(), local_hour(r)) for t, (r, _) in zip(timearr, moon_riseset) if r is not None]
mset = [(t.to_datetime(), local_hour(s)) for t, (_, s) in zip(timearr, moon_riseset) if s is not None]
print(f"Days without a moonrise: {sum(r is None for r, _ in moon_riseset)}")
print(f"Days without a moonset: {sum(s is None for _, s in moon_riseset)}")
fig, ax = plt.subplots(figsize=(8, 5))
ax.scatter(*zip(*mrise), s=4, color="#006450", label="Moonrise")
ax.scatter(*zip(*mset), s=4, color="#b35900", label="Moonset")
ax.set_xlabel("Date")
ax.set_ylabel("Local Hour of Day")
ax.set_title("Moonrise and Moonset Times for 2024 in Arlington, MA")
ax.set_ylim(0, 24)
ax.set_yticks([0, 6, 12, 18, 24])
ax.set_yticklabels(["Midnight", "6 AM", "Noon", "6 PM", "Midnight"])
ax.legend(loc="upper right")
fig.autofmt_xdate()
plt.tight_layout()
plt.show()
Days without a moonrise: 13 Days without a moonset: 12
Moon Phases¶
satkit.moon.phase_times lists the principal phases (New Moon, First Quarter, Full Moon, Last Quarter) in a time interval: the instants when the Moon's apparent geocentric ecliptic longitude exceeds the Sun's by 0°, 90°, 180° and 270°, as defined by the Astronomical Almanac and the USNO. With use_jpl=True the times agree with the USNO's to within their one-minute rounding; the built-in analytic model is good to about half an hour. satkit.moon.next_phase finds the next occurrence of one phase.
# Full moons of 2024, in US Eastern time, from the JPL ephemeris
phases = sk.moon.phase_times(sk.time(2024, 1, 1), sk.time(2025, 1, 1), use_jpl=True)
for phase, t in phases:
if phase == sk.moon.moonphase.FullMoon:
print(t.to_datetime().astimezone(eastern).strftime("%a %b %d %Y, %I:%M %p %Z"))
new = sk.moon.next_phase(sk.time(2024, 7, 4), sk.moon.moonphase.NewMoon, use_jpl=True)
print(f"First new moon on or after 4 July 2024: {new}")
Thu Jan 25 2024, 12:54 PM EST Sat Feb 24 2024, 07:30 AM EST Mon Mar 25 2024, 03:00 AM EDT Tue Apr 23 2024, 07:48 PM EDT Thu May 23 2024, 09:53 AM EDT Fri Jun 21 2024, 09:07 PM EDT Sun Jul 21 2024, 06:17 AM EDT Mon Aug 19 2024, 02:25 PM EDT Tue Sep 17 2024, 10:34 PM EDT Thu Oct 17 2024, 07:26 AM EDT Fri Nov 15 2024, 04:28 PM EST Sun Dec 15 2024, 04:01 AM EST First new moon on or after 4 July 2024: 2024-07-05T22:57:24.179535Z