Document how to find earliest and latest sunrise and sunset

[meiphoo7-Mae]

I don't know if this is the right place to ask this question, however I'm not aware of a better place. If this is not supposed to be here I apologize in advance.

I'm currently trying to integrate the marvels of the universe into my home automation system. Mostly because it is fun and I want to improve my coding skills. I managed to integrate the more common events like risings and settings of objects and the things described in the more well documented features of Skyfield. Now I'm trying to calculate less common things and I'm stuck.

For example: I want to figure out the most early and most late occurring times (and dates) of sunrise and sunset in the course of a year. I can do the calculation using the almanac.find_risings and almanac.find_settings routines and iterating over it using regular Python. But I want to know if it is faster to use find_minima and find_maxima. Somehow I simply don't understand how to do it. The documentation about these functions is not very helpful, because it is mostly about finding elongations. As a matter of fact I haven't been able to find an example on the Internet that is not related to finding elongations. I even tried a session with an AI-digital assistant which produced multiple versions of a script and each of them failing because of problems, mostly related to errors in the time(s) provided to the script. Interestingly, I do know how to perform calculations using find_discrete, but apparently using find_minima and find_maxima is more complicated. So, yes, I have to admit I can't figure this out on my own and a bit of help is very welcome.

[brandon-rhodes]

Happily, you don't need any special routines to find the answer in Python to a question like "what's the biggest number in this list?" or "what's the smallest number?" You can just ask NumPy where the biggest or smallest element is, and it will tell you. You can ignore things like find_minima(), which are for zeroing in on the place in a continuous curve that's lowest. With sunrises and sunsets, you have a simple list of 365 (or 366) numbers, and need to select the lowest.

It's important not to compare dates, because of course the "least" date in the list of sunrises will be the one on January 1, and the "maximum" will be the sunrise on December 31. Instead, you want to compare times. Happily, the Python datetime object has a .time() method with throws the year, month, and day away and just returns the hours, minutes, and seconds.

You could almost do something like this:

import numpy as np
from skyfield import almanac
from skyfield.api import N, W, load, wgs84
from pytz import timezone

ts = load.timescale()
eastern = timezone('US/Eastern')
eph = load('de421.bsp')
sun = eph['sun']
earth = eph['earth']
observer = earth + wgs84.latlon(40.8939 * N, 83.8917 * W)

t0 = ts.utc(2024, 1, 1)
t1 = ts.utc(2025, 1, 1)
t, y = almanac.find_risings(observer, sun, t0, t1)

datetimes = t.astimezone(eastern)
times = [dt.time() for dt in datetimes]

i = np.argmax(times)
print(i)
print(t[i])
print(datetimes[i])

Simple enough, right? Generate the list of sunrises as shown in the Skyfield docs, then ask which one is latest — that is, ask which sunrise has the maximum time-of-day compared to the others.

The problem is, it gives the wrong answer.

306
<Time tt=2460617.0064133313>
2024-11-02 08:08:04.927830-04:00

It says the maximum-timed = latest-timed sunrise is on November 2, which is nonsense. What's going wrong?

As always, a graph should make things clear.

import numpy as np
from skyfield import almanac
from skyfield.api import N, W, load, wgs84
from pytz import timezone

ts = load.timescale()
eastern = timezone('US/Eastern')
eph = load('de421.bsp')
sun = eph['sun']
earth = eph['earth']
observer = earth + wgs84.latlon(40.8939 * N, 83.8917 * W)

t0 = ts.utc(2024, 1, 1)
t1 = ts.utc(2025, 1, 1)
t, y = almanac.find_risings(observer, sun, t0, t1)

datetimes = t.astimezone(eastern)
times = [dt.time() for dt in datetimes]

import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.plot([t.hour + t.minute / 60.0 for t in times])
ax.grid()
fig.savefig('tmp.png')

Which generates:

It very nearly looks like a graph of sunrise times: around 8 AM in the middle of winter at the beginning of the graph, then almost at 6 AM by the middle of summer in the middle, then later again through the autumn and into winter. (Remember that this graph goes from January to December, because its 365-day x-axis is simply the days of the year.)

The problem, obviously, is Daylight Savings Time. Twice a year, it messes up the graph with a sharp spike, once when clocks "spring forward", and once when they "fall back". And just before they "fall back" in November, the DST scheme allows a sunrise that's so late in the day — at 8:08 AM, as we saw in the printout above — that it beats all of the mid-winter sunsets by being later in the day than any of them.

So, I guess this answer is accurate if you really mean to ask, "which sunrise occurs when my clock says the latest time?"

But usually when people ask about the "latest sunrise", they mean relative to Standard Time, not DST. Happily, pytz has a set of timezones that are fixed relative to UTC instead of having Daylight Savings Time. So:

import numpy as np
from skyfield import almanac
from skyfield.api import N, W, load, wgs84
from pytz import timezone

ts = load.timescale()
eastern = timezone('Etc/GMT-5')
eph = load('de421.bsp')
sun = eph['sun']
earth = eph['earth']
observer = earth + wgs84.latlon(40.8939 * N, 83.8917 * W)

t0 = ts.utc(2024, 1, 1)
t1 = ts.utc(2025, 1, 1)
t, y = almanac.find_risings(observer, sun, t0, t1)

datetimes = t.astimezone(eastern)
times = [dt.time() for dt in datetimes]

i = np.argmax(times)
print(i)
print(t[i])
print(datetimes[i])

The output indeed names a plausible date for latest-sunrise:

3
<Time tt=2460314.0426957025>
2024-01-04 18:00:19.724704+05:00

So, to summarize: finding the biggest or smallest item in a list is easy in Python, either by writing a little loop yourself that remembers the smallest or biggest item so far, or by calling a pre-written loop of the same sort like argmax() or argmin(). The only subtle thing here is how to measure time-of-day in a way that doesn't jump around and ruin your answer.

Try this out yourself, and see if it works!

[meiphoo7-Mae]

Thanks for your response. Your example code is about the same as what I did. I asked about find_maxima and find_minima mainly because of bad performance of my code. Unfortunately there were two problems: 1) find_maxima and find_minima only work on a continuous range of values, as you've pointed out. 2) In my code I accidently put the almanac formula within the main for-loop, which is very bad for performance of the code. I only noticed that mistake after your reply.

Remains the observation that I still don't know how to use find_minima and find_maxima. But I've to find a use case for that first.

[brandon-rhodes]

I am glad you figured out how to speed up your code, and that you understand now why those two routines won't help you find an earliest or latest sunrise or sunset.

I think I should add this example to the Skyfield "examples" docs, since getting a good answer turned out to be a bit more subtle than I expected. I've tagged this issue a "documentation request", and I'll close it once I get my discussion above pasted over into the docs somewhere.

[meiphoo7-Mae]

I'm glad to know that this question serves a certain purpose and that it does not end up as a waste of time.

In the meantime I continued programming and I found out that a lot of things can be calculated using plain Python or the find_discrete routine. However I believe I've found a use case for find_minima or find_maxima that is not related to elongation. I want to figure out when (one of) the outer planets starts and ends the 'opposition loop'. I don't know if this is the correct terminology, but it is a translation of the terminology used in the language of my country. What I mean is the phenomenon that an outer planet starts to move retrograde relative to the stars and some time later resumes its usual motion.

[brandon-rhodes]

However I believe I've found a use case for find_minima or find_maxima that is not related to elongation. I want to figure out when (one of) the outer planets starts and ends the 'opposition loop'. I don't know if this is the correct terminology, but it is a translation of the terminology used in the language of my country.

Oh, yes, good idea — those routines would be perfect for finding the moments when the planet turns around at each end of its retrograde loop.

I would have thought that there would be a great technical term for those moments, but looking at Guy Ottewell's Astronomical Calendar, which always gives the dates of those events, I see that they are simply described using words — "stationary in longitude" and "stationary in right ascension":

(This is from the 2001 edition.) Maybe I should someday think of adding routines for this to Skyfield's almanac module.