How to get earth vector data from jpl horizons in python?

Question

I am trying to get the vector data for Earth using Astroquery's Horizons Class. I have the following code:

from astroquery.jplhorizons import Horizons
import numpy as np

earth = Horizons(id=399, epochs = {'start':'2005-06-20', 'stop':'2005-06-21','step':'1d'})
earthVectors = earth.vectors()
earthX = earthVectors['x'].data # X is in AU
au2km = 149_597_870.7
earthXkm = earthX * au2km # X is in km

which returns earthXkm = [-3429775.6506088143 -899299.0538429054] in kilometers. Getting this information directly from JPL Hoizons gives [-2793030.0, -2627770.0] kilometers.

There is a large discrepancy here and this is the same for all the values in the astropy table. I would also not expect the data to vary as much in one day as that from the astroquery result.

Is there an error in my code, or does the horizons vectors() method not work as intended?

Answer 1

After further analysis I have found the cause of the discrepancy. Astroquery's Horizon Class uses a default coordinate system centered at the center of the sun for vectors. The Horizon's app; however, using the Solar System barycenter as the default coordinate origin. Using the location attribute set to the solar system barycenter fixes the issue. location='@ssb' or location='500@0'

Answer 2

You could just use astropy's get_body_barycentric instead (Note that Horizons currently uses the DE441 ephemeris, so the code below will download the file for this ephemeris which is 3.3 Gb):

from astropy.coordinates import get_body_barycentric, solar_system_ephemeris
from astropy.time import Time

# set the ephemeris to use DE441
solar_system_ephemeris.set("ftp://ssd.jpl.nasa.gov/pub/eph/planets/bsp/de441.bsp")

t = Time("2005-06-20", scale="tdb")
pos = get_body_barycentric("earth", t)

print(pos.x)
<Quantity -2793031.73765342 km>

This is identical (to within a micron [probably just the numerical truncation]) to what I get from the Horizons web interface (outputting value in km rather than AU).