Examples¶
Many programmers work most efficiently when they can start with a working example and adapt it to their needs. Each of the following examples pulls together several Skyfield features to solve a general problem, that should provide readers with a basis for solving other similar problems of their own.
What time is solar noon, when the Sun transits the meridian?¶
The Earth travels fastest in its orbit when closest to the Sun in January, and slowest in July when they are farthest apart. This stretches out January days by around 10 seconds, shortens July days by the same amount, and means that your watch never reads exactly 12:00 when the Sun transits your meridian at “solar noon”.
You can compute solar noon by asking at what time the Sun transits the meridian at your location. To learn the time of solar noon today, for example, you might do this:
import datetime as dt
from pytz import timezone
from skyfield import almanac
from skyfield.api import N, W, wgs84, load
zone = timezone('America/New_York')
now = zone.localize(dt.datetime.now())
midnight = now.replace(hour=0, minute=0, second=0, microsecond=0)
next_midnight = midnight + dt.timedelta(days=1)
ts = load.timescale()
t0 = ts.from_datetime(midnight)
t1 = ts.from_datetime(next_midnight)
eph = load('de421.bsp')
bluffton = wgs84.latlon(40.8939 * N, 83.8917 * W)
f = almanac.meridian_transits(eph, eph['Sun'], bluffton)
times, events = almanac.find_discrete(t0, t1, f)
# Select transits instead of antitransits.
times = times[events == 1]
t = times[0]
tstr = str(t.astimezone(zone))[:19]
print('Solar noon:', tstr)
Solar noon: 2020-04-19 13:34:33
When will it get dark tonight?¶
Sunrise, sunset, and the several varieties of twilight are all available through the Almanac Computation module. Here’s the script I use when I want to know when it will be dark enough to see the stars — or how early I need to rise to see the morning sky:
import datetime as dt
from pytz import timezone
from skyfield import almanac
from skyfield.api import N, W, wgs84, load
# Figure out local midnight.
zone = timezone('America/New_York')
now = zone.localize(dt.datetime.now())
midnight = now.replace(hour=0, minute=0, second=0, microsecond=0)
next_midnight = midnight + dt.timedelta(days=1)
ts = load.timescale()
t0 = ts.from_datetime(midnight)
t1 = ts.from_datetime(next_midnight)
eph = load('de421.bsp')
bluffton = wgs84.latlon(40.8939 * N, 83.8917 * W)
f = almanac.dark_twilight_day(eph, bluffton)
times, events = almanac.find_discrete(t0, t1, f)
previous_e = f(t0).item()
for t, e in zip(times, events):
tstr = str(t.astimezone(zone))[:16]
if previous_e < e:
print(tstr, ' ', almanac.TWILIGHTS[e], 'starts')
else:
print(tstr, ' ', almanac.TWILIGHTS[previous_e], 'ends')
previous_e = e
2020-04-19 05:09 Astronomical twilight starts
2020-04-19 05:46 Nautical twilight starts
2020-04-19 06:20 Civil twilight starts
2020-04-19 06:49 Day starts
2020-04-19 20:20 Day ends
2020-04-19 20:48 Civil twilight ends
2020-04-19 21:23 Nautical twilight ends
2020-04-19 22:00 Astronomical twilight ends
As you can see from the above code,
if the new light level is brighter
then we say that the new level “starts”,
but if the new level is darker
then we say the previous level “ends” —
so instead of saying “astronomical twilight starts at 21:23”
we say “nautical twilight ends at 21:23.”
That’s why the code keeps up with previous_e
and compares it to the new level of twilight.
What phase is the Moon tonight?¶
The phase of the Moon is defined as the angle between the Moon and the Sun along the ecliptic. This angle is computed as the difference in the ecliptic longitude of the Moon and of the Sun. The result is an angle that is 0° for the New Moon, 90° at the First Quarter, 180° at the Full Moon, and 270° at the Last Quarter. Skyfield also has a method for computing what fraction of a spherical body is illuminated by the sun.
from skyfield.api import load
from skyfield.framelib import ecliptic_frame
ts = load.timescale()
t = ts.utc(2019, 12, 9, 15, 36)
eph = load('de421.bsp')
sun, moon, earth = eph['sun'], eph['moon'], eph['earth']
e = earth.at(t)
s = e.observe(sun).apparent()
m = e.observe(moon).apparent()
_, slon, _ = s.frame_latlon(ecliptic_frame)
_, mlon, _ = m.frame_latlon(ecliptic_frame)
phase = (mlon.degrees - slon.degrees) % 360.0
percent = 100.0 * m.fraction_illuminated(sun)
print('Phase (0°–360°): {0:.1f}'.format(phase))
print('Percent illuminated: {0:.1f}%'.format(percent))
Phase (0°–360°): 149.4
Percent illuminated: 92.9%
What is the angular diameter of a planet, given its radius?¶
Be careful to select the correct radius when predicting a planet’s angular diameter in the sky. Many web sites will quote some kind of “mean radius” that averages between a planet’s squat polar radius and its wide equatorial radius. But most astronomers instead want to know the maximum, not average, diameter across a planet’s visible face — so you will want to use the planet’s equatorial radius in your calculation.
For example, a good current estimate of Neptune’s equatorial radius is 24,764 km. We would therefore predicts its angular diameter as:
import numpy as np
from skyfield.api import Angle, load
ts = load.timescale()
time = ts.utc(2020, 12, 30)
eph = load('de421.bsp')
earth, neptune = eph['earth'], eph['neptune barycenter']
radius_km = 24764.0
astrometric = earth.at(time).observe(neptune)
ra, dec, distance = astrometric.apparent().radec()
apparent_diameter = Angle(radians=np.arcsin(radius_km / distance.km) * 2.0)
print('{:.6f} arcseconds'.format(apparent_diameter.arcseconds()))
2.257190 arcseconds
This agrees exactly with the output of the NASA HORIZONS system.
When is Venus at its greatest east and west elongations from the Sun?¶
This example illustrates the several practical steps that are often required to both find events of interest and then to learn more details about them.
The concept of “elongation from the Sun” is here explained to Skyfield with a function that for any given time
t
returns the separation between the Sun and Venus in the sky.The
find_maxima()
routine is then set loose to find the moments over the 3 years 2019–2021 at which the elongation of Venus from the Sun reaches its maximum values. Skyfield starts by computing the elongation everystep_days
= 15 days between the search’s start time and end time, then hones in everywhere it sees a local maximum: a value that’s bigger than either of the two values next to it.Finally, a
for
loop over the results not only displays each maximum but computes and displays an extra fact: whether the elongation is an east or west maximum elongation, which is defined as whether Venus’s apparent ecliptic longitude is ahead of or behind the Sun’s along the great circle of the ecliptic.
This example can serve as a template for many other kinds of custom search:
from skyfield.api import load
from skyfield.framelib import ecliptic_frame
from skyfield.searchlib import find_maxima
ts = load.timescale()
t0 = ts.utc(2019)
t1 = ts.utc(2022)
eph = load('de421.bsp')
sun, earth, venus = eph['sun'], eph['earth'], eph['venus']
def elongation_at(t):
e = earth.at(t)
s = e.observe(sun).apparent()
v = e.observe(venus).apparent()
return s.separation_from(v).degrees
elongation_at.step_days = 15.0
times, elongations = find_maxima(t0, t1, elongation_at)
for t, elongation_degrees in zip(times, elongations):
e = earth.at(t)
_, slon, _ = e.observe(sun).apparent().frame_latlon(ecliptic_frame)
_, vlon, _ = e.observe(venus).apparent().frame_latlon(ecliptic_frame)
is_east = (vlon.degrees - slon.degrees) % 360.0 < 180.0
direction = 'east' if is_east else 'west'
print('{} {:4.1f}° {} elongation'.format(
t.utc_strftime(), elongation_degrees, direction))
2019-01-06 04:53:35 UTC 47.0° west elongation
2020-03-24 22:13:32 UTC 46.1° east elongation
2020-08-13 00:14:12 UTC 45.8° west elongation
2021-10-29 20:51:56 UTC 47.0° east elongation
Are planets separated by 0° at conjunction and 180° at opposition?¶
It surprises many newcomers to astronomy that the angular separation between two planets never drops all the way to 0° at conjunction nor increases all the way to a full 180° at opposition. The reason is that the planets will still have at least a slight lingering difference in latitude, even when their longitudes have brought them together or have placed them opposite each other in the sky.
We can take as an example
the date and time of the conjunction between Mars and the Sun
computed in the Opposition and Conjunction section of the Almanac page.
How close are they in the sky at that moment?
The separation_from()
method
measures raw angular distance
between any two points in the sky:
from skyfield.api import load
from skyfield.framelib import ecliptic_frame
ts = load.timescale()
eph = load('de421.bsp')
sun, mars = eph['sun'], eph['mars']
t = ts.utc(2019, 9, 2, 10, 42, 26)
e = earth.at(t)
s = e.observe(sun).apparent()
m = e.observe(mars).apparent()
print('{:.5f}°'.format(m.separation_from(s).degrees))
1.08256°
They are more than one degree apart!
How can that be,
if their ecliptic longitudes are at that moment the same?
Let’s use Skyfield’s ecliptic_frame
to express their positions in Ecliptic coordinates:
print(' Latitude Longitude')
lat, lon, distance = s.frame_latlon(ecliptic_frame)
print('Sun {:.5f}° {:.5f}°'.format(lat.degrees, lon.degrees))
lat, lon, distance = m.frame_latlon(ecliptic_frame)
print('Mars {:.5f}° {:.5f}°'.format(lat.degrees, lon.degrees))
Latitude Longitude
Sun 0.00005° 159.68641°
Mars 1.08260° 159.68641°
While the Sun sits very close to the ecliptic —
as we would expect, since the ecliptic is defined
as the course the Sun takes around the sky each year —
the inclination of the orbit of Mars has carried it
more than one degree above the ecliptic.
That’s why the separation_from()
method
still measured an angle of more than one degree between them.
A similar situation pertains at opposition:
t = ts.utc(2020, 10, 13, 23, 25, 55)
e = earth.at(t)
s = e.observe(sun).apparent()
m = e.observe(mars).apparent()
print('Separation: {:.5f}°'.format(m.separation_from(s).degrees))
print('')
print(' Latitude Longitude')
lat, lon, distance = s.frame_latlon(ecliptic_frame)
print('Sun {:.5f}° {:.5f}°'.format(lat.degrees, lon.degrees))
lat, lon, distance = m.frame_latlon(ecliptic_frame)
print('Mars {:.5f}° {:.5f}°'.format(lat.degrees, lon.degrees))
Separation: 177.00424°
Latitude Longitude
Sun 0.00007° 201.07794°
Mars -2.99582° 21.07794°
Even though their ecliptic longitudes are 180° apart,
the fact that neither the Sun nor Mars is lying exactly on the ecliptic
means that the separation_from()
method
finds that they are not quite 180° apart.
In case you run across the term ‘elongation’ in discussions of conjunctions and oppositions, it’s shorthand for ‘the angle between a planet and the Sun’ — and so each of the angular separations printed above can, more specifically, be labeled as the ‘elongation of Mars’ on those dates.
At what angle in the sky is the crescent Moon?¶
The angle of the crescent Moon changes with the seasons. In the spring, a crescent Moon will stand high above the Sun and appear to be lit from below. In the autumn, the Moon sets farther from the Sun along the horizon and is illuminated more from the side. What if we wanted to know the exact angle?
You can find the answer by asking for the Sun’s “position angle” relative to the Moon, an angle you can compute between any two Skyfield positions. The angle will be 90° if the Sun is left of the moon, 180° if the Sun is directly below, and 270° if the Sun is to the right of the Moon.
from skyfield.api import N, W, load, wgs84
from skyfield.trigonometry import position_angle_of
ts = load.timescale()
t = ts.utc(2019, 9, 30, 23)
eph = load('de421.bsp')
sun, moon, earth = eph['sun'], eph['moon'], eph['earth']
boston = earth + wgs84.latlon(42.3583 * N, 71.0636 * W)
b = boston.at(t)
m = b.observe(moon).apparent()
s = b.observe(sun).apparent()
print(position_angle_of(m.altaz(), s.altaz()))
238deg 55' 55.3"
The position_angle_of()
routine
will not only accept
the output of altaz()
,
but also of frame_latlon()
if you want a position angle relative to the ecliptic’s north pole.
Beware, though, that radec()
produces coordinates in the opposite order
from what position_angle_of()
expects:
right ascension is like longitude, not latitude.
Try reversing the coordinates, like:
print(position_angle_of(m.radec(), s.radec()))
282deg 28' 15.7"
Drat, but this angle is backwards, because right ascension increases toward the east whereas the other angles, like azimuth, increase the other way around the circle.
When is a body or fixed coordinate above the horizon?¶
The following code will determine
when the Galactic Center is above the horizon.
The Galactic Center is an example of a fixed object,
like a star or nebula or galaxy,
whose right ascension and declination can be plugged in to a Star()
object.
The code will also work with a body from an ephemeris,
like the Sun, Moon, or one of the planets.
from skyfield.api import N, Star, W, wgs84, load
from skyfield.almanac import find_discrete, risings_and_settings
from pytz import timezone
ts = load.timescale()
t0 = ts.utc(2019, 1, 19)
t1 = ts.utc(2019, 1, 21)
moab = wgs84.latlon(38.5725 * N, 109.54972238 * W)
eph = load('de421.bsp')
gc = Star(ra_hours=(17, 45, 40.04), dec_degrees=(-29, 0, 28.1))
f = risings_and_settings(eph, gc, moab)
tz = timezone('US/Mountain')
for t, updown in zip(*find_discrete(t0, t1, f)):
print(t.astimezone(tz).strftime('%a %d %H:%M'), 'MST',
'rises' if updown else 'sets')
Sat 19 05:51 MST rises
Sat 19 14:27 MST sets
Sun 20 05:47 MST rises
Sun 20 14:23 MST sets
At what rate is a target moving across the sky?¶
Automatically-driven telescopes and antennas often need to know the rate at which an object is moving across the sky. Specifically, an instrument with a simple altazimuth mount will need to know the rates at which altitude and azimuth are changing, whereas a fancier equatorial mount will need the rates for right ascension and declination.
The solution is the same in both cases:
to call Skyfield’s
frame_latlon_and_rates()
position method
and pass it the frame of reference
in which you want the rates computed.
In either case,
you will want to start by asking Skyfield to compute an apparent position:
from skyfield.api import load, wgs84, N, E
ts = load.timescale()
t = ts.utc(2021, 2, 3, 0, 0)
planets = load('de421.bsp')
earth, mars = planets['earth'], planets['mars']
top = wgs84.latlon(35.1844866 * N, 248.347300 * E, elevation_m=2106.9128)
a = (earth + top).at(t).observe(mars).apparent()
Since every topocentric location in Skyfield
is itself a reference frame representing the location’s horizon,
we can simply pass top
to the
frame_latlon_and_rates()
method
to learn the rates at which the altitude and azimuth are changing:
alt, az, distance, alt_rate, az_rate, range_rate = (
a.frame_latlon_and_rates(top)
)
print('Alt: {:+.2f} asec/min'.format(alt_rate.arcseconds.per_minute))
print('Az: {:+.2f} asec/min'.format(az_rate.arcseconds.per_minute))
Alt: +548.66 asec/min
Az: +1586.48 asec/min
But if our instrument is on an equatorial mount,
it will need to know how fast
the right ascension and declination are changing.
In that case we import Skyfield’s reference frame
that represents the Earth’s true orientation in space,
and pass that object to the
frame_latlon_and_rates()
method:
from skyfield import framelib
eeod = framelib.true_equator_and_equinox_of_date
dec, ra, distance, dec_rate, ra_rate, range_rate = (
a.frame_latlon_and_rates(eeod)
)
print(f'RA: {ra_rate.arcseconds.per_hour:+.2f} asec/hr')
print(f'Dec: {dec_rate.arcseconds.per_hour:+.2f} asec/hr')
RA: +78.66 asec/hr
Dec: +25.61 asec/hr
Note that, contrary to Skyfield’s usual custom,
declination is returned first.
That’s why the
frame_latlon_and_rates()
method
has latlon
in its name:
it always returns the latitude-like coordinate
(angle above or below the plane) first,
then the longitude-like coordinate
(angle around the plane) second.
Finally, in case you need it,
you will notice that both of the calls above return a range_rate
that is positive if the body is moving away from you
and its range is increasing,
or is negative if the target is moving closer and the range is falling.
At the specific date and time we asked about above,
Mars is moving away:
print('Range rate: {:+.2f} km/s'.format(range_rate.km_per_s))
Range rate: +16.79 km/s
Finally,
there’s the slight complication
that some data sources and instruments
want the rate of motion around-the-sky
to be multiplied by the cosine of the body’s angle
above or below the reference plane.
By simply performing the math—and
remembering that sin()
and cos()
in Python take radian arguments—you
can produce these alternative measurements yourself:
from numpy import cos
rcos = az_rate.arcseconds.per_minute * cos(alt.radians)
print(
'Azimuth rate × cos(altitude): {:.2f} arcseconds / minute'
.format(rcos)
)
rcos = ra_rate.arcseconds.per_hour * cos(dec.radians)
print(
'RA rate × cos(declination): {:.2f} arcseconds / hour'
.format(rcos)
)
Azimuth rate × cos(altitude): 663.55 arcseconds / minute
RA rate × cos(declination): 75.16 arcseconds / hour
What is the right ascension and declination of a point in the sky?¶
An observer is often interested in the astronomical coordinates of a particular position in the sky above them. If the observer can specify the position using altitude and azimuth coordinates, then Skyfield can return its right ascension and declination.
from skyfield import api
ts = api.load.timescale()
t = ts.utc(2019, 9, 13, 20)
geographic = api.wgs84.latlon(latitude_degrees=42, longitude_degrees=-87)
observer = geographic.at(t)
pos = observer.from_altaz(alt_degrees=90, az_degrees=0)
ra, dec, distance = pos.radec()
print(ra)
print(dec)
13h 41m 14.65s
+42deg 05' 50.0"
What latitude and longitude is beneath this right ascension and declination?¶
Most Skyfield calculations,
like an observation of a planet or an Earth satellite,
directly produce a vector position centered on the Earth.
You can pass such a vector
to the subpoint()
method
of a standard geoid to compute latitude and longitude.
But sometimes the right ascension and declination of the position
are known already.
Instead of creating a Star
with those coordinates
and asking it to compute its position,
there is a simpler approach:
creating the position directly.
from skyfield.api import load, wgs84
from skyfield.positionlib import position_of_radec
ts = load.timescale()
t = ts.utc(2020, 1, 3, 12, 45)
earth = 399 # NAIF code for the Earth center of mass
ra_hours = 3.79
dec_degrees = 24.1167
pleiades = position_of_radec(ra_hours, dec_degrees, t=t, center=earth)
subpoint = wgs84.subpoint(pleiades)
print('Latitude:', subpoint.latitude)
print('Longitude:', subpoint.longitude)
Latitude: 24deg 10' 33.5"
Longitude: 123deg 16' 53.9"
Which geographic location is farther from Earth’s center?¶
After I hiked Mount Bierstadt in Colorado, a friend suggested that its 14,000 feet of elevation might have carried me farther from the Earth’s center than I had ever traveled before. It was a romantic thought: that under my own power I had hiked farther from my home planet’s core than ever before.
But there was a problem. I knew that I had once visited a city only a few degrees away from the Earth’s equator, and that the Earth’s equatorial bulge might push even modest elevations at that latitude out farther from the Earth’s center than a mountaintop in Colorado.
So I wrote a quick Skyfield script to compare the distance from the Earth’s center to both Accra, Ghana, and the top of Mount Bierstadt in Colorado.
from skyfield.api import N, W, wgs84, load
from skyfield.functions import length_of
ts = load.timescale()
t = ts.utc(2019, 1, 1)
bierstadt = wgs84.latlon(39.5828 * N, 105.6686 * W, elevation_m=4287.012)
m1 = length_of(bierstadt.at(t).position.m)
print(int(m1))
accra = wgs84.latlon(5.6037 * N, 0.1870 * W, elevation_m=61)
m2 = length_of(accra.at(t).position.m)
print(int(m2))
assert m2 > m1
print("I was", int(m2 - m1), "meters farther from the Earth's center\n"
"when I visited Accra, at nearly sea level, than atop\n"
"Mt. Bierstadt in Colorado.")
6373784
6377995
I was 4211 meters farther from the Earth's center
when I visited Accra, at nearly sea level, than atop
Mt. Bierstadt in Colorado.