Earth coordinate at Sun Zenith & Day Night Boundary

[hamiltoncj]

Can Skyfield be used to find the coordinate where the sun is at its zenith for a particular date and time? The second part of this question is can it be used to calculate the day/night terminator coordinates (lat, lon) so that I can produce a polygon of the night area in QGIS. Any help would be appreciated.

[JoshPaterson]

This question nerd-sniped me pretty good, haha. Here is a script I wrote that I believe does what you need. The sun subpoint is quite easy to calculate since Skyfield has a feature for calculating subpoints. Finding points on the terminator required some custom vector math. I tried to be liberal with comments, but if you have questions about how it works let me know.

from skyfield.api import load, wgs84
import numpy as np
from skyfield.positionlib import Geocentric

eph = load('de421.bsp')
earth = eph['earth'] # vector from solar system barycenter to geocenter
sun = eph['sun'] # vector from solar system barycenter to sun
geocentric_sun = sun - earth # vector from geocenter to sun

ts = load.timescale()
t = ts.utc(2021, 5, 19, 12, 26)

sun_subpoint = wgs84.subpoint(geocentric_sun.at(t)) # subpoint method requires a geocentric position

print('subpoint latitude: ', sun_subpoint.latitude.degrees)
print('subpoint longitude: ', sun_subpoint.longitude.degrees)


def cross_product_matrix(a):
    """ computes the cross product matrix

    see https://en.wikipedia.org/wiki/Cross_product#Conversion_to_matrix_multiplication
    code from https://stackoverflow.com/questions/66707295/numpy-cross-product-matrix-function
    """
    return np.cross(a, np.identity(3) * -1)

def rotation_matrix_around_axis(axis_vector, rotation_degrees):
    """ creates a rotation matrix that rotates around a given axis

    see https://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle
    """
    rotation_radians = rotation_degrees / 180 * np.pi
    return (
        np.cos(rotation_radians) * np.identity(3)
        + np.sin(rotation_radians) * cross_product_matrix(axis_vector)
        + (1 - np.cos(rotation_radians)) * np.outer(axis_vector, axis_vector)
        )

terminator_angle_from_sun = 90.833 # 90 degrees + sun's semidiameter + refraction at horizon

sun_vec = geocentric_sun.at(t).position.au # numpy array of sun's position vector
normal_vec = np.cross(sun_vec, np.array([1,0,0]))# vector normal to sun position and x-axis
first_terminator_vec = rotation_matrix_around_axis(normal_vec, terminator_angle_from_sun) @ sun_vec # arbitrary first position on terminator

terminator_latitudes = []
terminator_longitudes = []

num_points_on_terminator = 100
for angle in np.linspace(0, 360, num_points_on_terminator):
    terminator_vector = rotation_matrix_around_axis(sun_vec, angle) @ first_terminator_vec
    terminator_position = Geocentric(terminator_vector, t=t)
    geographic_position = wgs84.subpoint(terminator_position)
    terminator_latitudes.append(geographic_position.latitude.degrees)
    terminator_longitudes.append(geographic_position.longitude.degrees)

terminator_latitudes = np.array(terminator_latitudes)
terminator_longitudes = np.array(terminator_longitudes)

import matplotlib.pyplot as plt

fig, ax = plt.subplots()
ax.scatter(terminator_longitudes, terminator_latitudes)
ax.scatter(sun_subpoint.longitude.degrees, sun_subpoint.latitude.degrees)
ax.grid(True)
plt.show()

for me this produces this plot:

[hamiltoncj]

@JoshPaterson This code only works close to properly on the summer and winter solstice. I thought I should share what I have done so far in QGIS. This displays the location of the sun and moon as well as the regions of civil twilight, nautical twilight, astronomical twilight, and night.

[brandon-rhodes]

Wow, that’s beautiful! Thanks for sharing! If you ever share the code more widely, feel free to add a link here for folks who want to read through it and see how the astronomical data is converted into QGIS-friendly coordinates.

[hamiltoncj]

Thanks! I have a few more capabilities I want to add to this before I release the code, but this will be an open source QGIS plugin and hopefully it will be added to the official QGIS plugin repo.

[JoshPaterson]

That looks great! Good job!

[hamiltoncj]

@JoshPaterson It appears I was too optimistic and was just testing the current date and time. If I run this for March 21, 2021 at 1300 UTM time, I get these results.

What should be displayed is something like the lightest boundary outer edge of this picture. See: https://www.timeanddate.com/worldclock/sunearth.html?month=3&day=21&year=2021&hour=13&min=0&sec=0&n=&ntxt=&earth=0

Do you have any ideas what is wrong?

[hamiltoncj]

@JoshPaterson Hi Josh, It has been a while since I took a look at this code. I don't know if your are still around and interested. This code is close to correct at the summer and winter solstice, but not any other time. I think the rotation of the earth's axis is not being taken into consideration some how. If you are still around, do you think this problem could be solved. If I can get it working correctly, then I will add a day-night, twilight algorithm into QGIS. Thanks.