Discrepancy between PyEphem and Skyfield HA/dec results

Question

I'm refactoring some old code that used PyEphem to use Skyfield, and I'm getting a slight difference in the results of the GHA/dec of a body.

def sf(year):
    from skyfield.api import N, W, wgs84
    from skyfield.api import load
    from skyfield.units import Angle

    ts = load.timescale()
    eph = load('de421.bsp')
    greenwich = wgs84.latlon(58.47722 * N, 0.0 * W)
    t = ts.utc(year)
    e = eph['earth'] + greenwich
    v = eph['venus']
    ha = e.at(t).observe(v).hadec()
    gha = Angle(degrees=360.0 + ha[0]._degrees)

    print(f'{gha}, {ha[1]}')


def pe(year):
    import ephem

    greenwich = ephem.Observer()
    greenwich.lon = 0.0
    greenwich.lat = ephem.degrees('51:28:38')
    greenwich.pressure = 0.0
    greenwich.horizon = '-0:34'

    t = ephem.date(str(year))
    v = ephem.Venus()
    greenwich.date = t
    greenwich.epoch = t
    v.compute(greenwich)
    gha = ephem.twopi - v.g_ra + greenwich.sidereal_time()

    print(f'{ephem.degrees(gha)}, {ephem.degrees(v.g_dec)}')


if __name__ == '__main__':
    sf(2016)
    pe(2016)

Yields the results:

219deg 41' 57.5", -18deg 37' 03.9"
219:42:15.7, -18:36:56.0

The old PyEphem code agrees with the Nautical Almanac I have in front of me.

This is definitely a PEBKAC, but I'm scratching my head to find what temporal or spatial transformation I've missed.

Answer 1

It looks like you are asking Skyfield for astrometric positions, but PyEphem for apparent positions. According to the PyEphem documentation:

"g_ra and ra — Apparent geocentric right ascension for the epoch-of-date"

https://rhodesmill.org/pyephem/quick.html#body-compute-date

Whereas with Skyfield, you have to call the .apparent() method on a position to learn the corresponding apparent position; it does not happen automatically:

https://rhodesmill.org/skyfield/positions.html#barycentric-astrometric-apparent

See if that change eliminates most of the difference between coordinates.