Request for day-night terminator line function

[hamiltoncj]

It would be really nice to have a skyfield function that given a UTC date and time returns a list of coordinates that represent day-night terminator line across the earth's surface. It would also be nice to calculate the different twilight lines as well. One of the parameters could be the number of coordinates to be returned. Thanks!

[JonathanMaes]

I just came across this issue, and as luck would have it I recently compiled a small script scrambled together from various sources to calculate the terminator line, see this gist: https://gist.github.com/JonathanMaes/52deb6c333c1033d5f5d22a25529eb31

You can get the different twilight lines by changing the refraction parameter, e.g. a value of approximately 6° is the civil twilight line, 12° nautical, 18° astronomical. The resulting curves look like this:

(The main pieces of code doing the calculation in the gist are the __init__ and _rotate_terminator_vertices methods. The rest is for plotting, which I included so you can easily test the result yourself. The _solar_position method can probably be replaced by more accurate calls to the skyfield.api, but I'm still new to skyfield so I did not yet use that when I wrote this script).

I just thought I'd drop this here in case if it would be helpful to this issue.

[brandon-rhodes]

@JonathanMaes — Thanks for pointing out your existing script! It looks like your script approximates the shape of the Earth as a sphere, without any flattening at the poles—do you know how much that affects the day/night line? It might not matter much if things like refraction make the terminator less certain than the effect flattening would make.

[JonathanMaes]

@brandon-rhodes I did a comparison between the terminator line calculated by the script, and the places where the sun is near the same elevation (-0.83°) calculated with skyfield and wgs84 on the other hand. Just to be sure, wgs84 does take the flattening at the poles into account, right?

It seems like, in the worst case, the line is accurate to within 2km.

This discrepancy is likely due to

• The terminator line being a bunch of line segments, so it is more accurate at a vertex than in the middle between two vertices
• The lat,lon where the sun is at zenith in the terminator script I shared above, is not calculated using skyfield but with a standalone approximation
• The oblateness of the earth?

(Some case studies:

• near the north pole (e.g. 77°N), the skyfield approach is less than 500m to the east of the terminator line calculated by the script
• near the equator (2°N) the skyfield approach is less than 200m to the west
• near the south pole (76°S), skyfield is about 1km to the west

The time of the year also does not seem to matter much.)

[brandon-rhodes]

Yes, wgs84 takes the oblateness as one of its arguments when it's built:

https://github.com/skyfielders/python-skyfield/blob/2e89cb5cc4e52644b65f4d029bc0f448e9109e45/skyfield/toposlib.py#L276

I wonder if building your own Geoid with a second argument of something like 1e10, which I think would make an almost perfect sphere, would let you determine how much of the error is from oblateness, and how much from something else.

[JonathanMaes]

From what I can see, changing it to such a Geoid does not change anything. Even the difference between a Geoid with the second argument 1e10 or just 1 is very small (about 50m at most). This means that the discrepancies in my previous reply are probably only due to the slightly different position of the sun in the simple gist script.

Could it be that place.observe(sun).apparent().altaz() does not take into account the slightly tilted horizon of the WGS84 globe?

For completeness, the core of my code that determines whether a point is in nighttime or not, is

from skyfield.api import load, wgs84
from skyfield.toposlib import Geoid
import pytz
# wgs84 = Geoid('WGS84', 6378137.0, 1e10) # Toggling this comment makes the earth round or oblate

t = load.timescale().from_datetime(datetime(2023, 6, 21, 10, 10, 10, tzinfo=pytz.UTC))
eph = load('de421.bsp')
sun, earth = eph['sun'], eph['earth']
place = earth + wgs84.latlon(LAT, LON)
elevation = place.at(t).observe(sun).apparent().altaz()[0].degrees
darkness = (el <= -0.83)
print(darkness)

[brandon-rhodes]

Just to be sure, wgs84 does take the flattening at the poles into account, right? … Could it be that place.observe(sun).apparent().altaz() does not take into account the slightly tilted horizon of the WGS84 globe?

Happily, every part of Skyfield takes into rigorous account the slightly tilted horizon of WGS84! Its results agree to very high accuracy with other astronomy libraries that use the same Earth model.

Thanks again for contributing this sample code—hopefully it will help users who need a decent day-night line, until the time comes that its code can be fully converted to using Skyfield position logic instead.

[hamiltoncj]

@JonathanMaes Jonathan, I am trying to adapt your code to work in the QGIS plugin Earth, Sun, Moon, and Planets (https://plugins.qgis.org/plugins/earthsunmoon/), at least until I can get a version that uses Skyfield. I have made one change in the code in "as_polygons". Rather than tacking on [180,-180] to the end of the array, I put [-180] and the beginning and [180] at the end. Without doing this it messes up the plots. In one case I just want the day/night terminator line and the rest of the time the polygon so I added in a check to see if I want a closed polygon generated or not. I add in extra coordinates to complete the polygon so that the drawing does not get confused in QGIS. Here are the changes for transitions.size == 1. I need to look at the other cases. transitions.size == 0 isn't working quite right. I'm not very good with NumPy so there may be a better way to make the changes.

        if transitions.size == 1: # sine across the globe (have to add two vertices at top or bottom to complete the night fill)
            p = x.argsort()
            x, y = x[p], y[p]
            y_level = -90 + 180*(solar_lat <= 0)
            x = np.append(np.insert(x, 0, [-180]), [180])
            y = np.append(np.insert(y, 0, [y_level]), [y_level])
            if return_polygon:
                if solar_lat > 0: # Sun is in the nothern hemisphere and the south will be dark
                    np.append(x, [180, 90, 0, -90, -180, -180])
                    np.append(y, [-90, -90, -90, -90, -90, y_level])
                else:
                    np.append(x, [180, 90, 0, -90, -180, -180])
                    np.append(y, [90, 90, 90, 90, 90, y_level])