Differences in GAST between astropy, skyfield, and pyephem

[aelanman]

I have some old code which uses pyephem to transform ICRS ra/dec positions to hour angle / local declination, and I'm trying to replace this with the equivalent in skyfield. I'm new to both, so I'm comparing with astropy to check results. I'm seeing slight differences in the local coordinates for each.

While tracking down the cause of those differences, I found that I get different values for the Greenwich apparent sidereal time (GAST) from each package. Between astropy and skyfield, I'm seeing a difference of a couple milliseconds. The following code does the comparison:

import ephem
from skyfield import api
import time as ptime
from datetime import datetime
from astropy import coordinates as ac, time as atime

unix_time = ptime.time()
if isinstance(unix_time, float):
    time = datetime.utcfromtimestamp(unix_time)
    time = time.replace(tzinfo=api.utc)

# Skyfield
TIMESCALE = api.load.timescale()
tt = TIMESCALE.from_datetime(time)
chime_loc = api.wgs84.latlon(
    latitude_degrees=0,
    longitude_degrees=0,
)
lst_sky = chime_loc.lst_hours_at(tt)

# Ephem
chime_eph = ephem.Observer()
chime_eph.lat = 0
chime_eph.long = 0
chime_eph.date = time
lst_eph = ac.Longitude(chime_eph.sidereal_time(), unit='rad').hour

# Astropy
t0 = atime.Time(time, scale='utc')
lst_apy = t0.sidereal_time('apparent', longitude='greenwich').hour

print(f"Astropy : {lst_apy:.7f} hours")
print(f"PyEphem : {lst_eph:.7f} hours")
print(f"Skyfield: {lst_sky:.7f} hours")

With the result:

Astropy : 8.8973182 hours
PyEphem : 8.8973459 hours
Skyfield: 8.8973119 hours

Any ideas what I might be doing wrong here?

Thanks for your help!

[brandon-rhodes]

Note that you can also ask Skyfield for the GAST with:

print(tt.gast)

To trace the difference between the two numbers, I would suggest reading through how the two different libraries compute the numbers, and seeing if they use formulae from different eras—every few years that IAU comes out with new polynomials, I believe. Glancing at the code called by the Time.gast() method in Skyfield, I'm not terribly impressed with the documentation—it looks like the references in the original code that state the formulae sources got cut off, saying things like from the reference's eq. (39) but without actually including the paper citation that the code was referencing. I'll go take a look at the USNO code that Skyfield was based on, and see if I can extend the docstrings to mention specific sources!

[brandon-rhodes]

There — I have just updated the docstrings of the functions underlying the GAST computation, so that you can track them down and compare them to the formulae used in AstroPy:

9573293

It does look like they date from the early 2000s. If you find that the results are not as accurate as your application needs, and that AstroPy's numbers are better, we could look into making the two agree better, so please do let me know what you discover!

[aelanman]

Thanks, @brandon-rhodes for the quick response, and for updating the documentation!

I've spotted a difference in fractional UT1 between skyfield and astropy. The attribute Time.ut1_fraction in skyfield differs from Time.jd2 in astropy by order 10 milliseconds. It's unclear if this is due to some difference in how they're computed or a precision limit.

Minimal example to see the difference:

from skyfield import api
from astropy.time import Time, TimeDelta
ts = api.load.timescale()
t_apy = Time.now()
t_sky = ts.from_astropy(t_apy)
t_apy.jd2 - t_sky.ut1_fraction
print(TimeDelta(t_apy.ut1.jd2 - t_sky.ut1_fraction, format='jd').sec)

This returns 0.02. This is then making a difference in ERA, which causes a difference in GMST and hence GAST.

There are some slight differences in the other terms that go into computing GMST, but those differences could just be due to numerical precision. Plugging in the same values of tt.whole and ut1_fraction in earthlib.earth_rotation_angle and in erfa.era06 yield the same result, which is good since they use the same underlying formula.

The precession, nutation, and frame bias matrices are the same between the two codes, and although they use slightly different methods to go from GMST to GAST (astropy/erfa using ERA - equation of the origin, skyfield using GMST - equation of the equinoxes), that doesn't seem to be causing a large deviation.

Could this be the main problem?

[My original message erroneously compared UTC with UT1. I've corrected this]

[brandon-rhodes]

@aelanman — Apologies for the delay in my response; travel and other obligations got in the way and only now am I catching up!

Yes, UT1 is probably the main culprit, because it's not a value produced by a formula, or not entirely: because we can't predict the Earth's motion precisely as the continents and oceans and mantle move and churn, the IERS has to simply publish numbers for where UT1 is, and occasionally issue a leap second to pull UTC back close to it.

Here's how Skyfield can be given a recent IERS table:

https://rhodesmill.org/skyfield/time.html#ut1-and-downloading-iers-data

See whether using a fresh copy of the file moves Skyfield closer to (or father away from!) the Astropy estimate for UT1.

[aelanman]

@brandon-rhodes Thanks for your reply. I've updated the IERS data table (however, I'm still having some trouble with issue #730) and I'm seeing better agreement. This plot shows the difference between astropy and skyfield for UT1 over 24 hours (using a slight variation of the code I sent before):

So it's improved by about a factor of 10. Maybe there's some difference in how skyfield and astropy interpolate the UTC-UT1 table? The discreteness does suggest that it's some numerical issue.