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.
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.