Astropy and JPL's Horizons query have different output

Question

I am trying to compare the position of the Earth in heliocentric system given by Astropy, with that of JPL's Horizons. I find remarkable differences and I don't know what they may be due....

I want to do this way because my purpose is to convert from geographic coordinates to heliocentric one.

    from astropy.coordinates import SkyCoord
    from astropy.time import Time
    from astropy import units as u
    from astroquery.jplhorizons import Horizons
    
    time = '2015-03-30T21:33:52.000'
    JD = Time(time).jd
    
    observing_time = Time(time, format='isot', scale='utc')
    
    #ASTROPY
    c = SkyCoord(x=0,y=0,z=0, unit='au', representation_type='cartesian', frame='gcrs', obstime=time)
    
    cc = c.transform_to("heliocentriceclipticiau76")
    print(cc.cartesian.x, cc.cartesian.y, cc.cartesian.z)
    
    cc = c.transform_to("heliocentricmeanecliptic")
    print(cc.cartesian.x, cc.cartesian.y, cc.cartesian.z)
    
    cc = c.transform_to("heliocentrictrueecliptic")
    print(cc.cartesian.x, cc.cartesian.y, cc.cartesian.z)
    
    #JPL'S HORIZONS
    obj = Horizons(id="Geocenter", location="@sun", epochs=JD, id_type='id').vectors()
    print(obj['x'], obj['y'], obj['z'])

Answer 1

Possibly this is no longer of interest after 2 years but for the benefit of future users who may get bitten by this subtlety, here's the answer. There is a subtle difference in the timescales that JPL HORIZONS, and by extension the astroquery.jplhorizons code, expects when asking for state vectors vs ephemerides. As documented under the epochs part of HorizonsClass (docs):

Epoch timescales depend on the type of query performed: UTC for ephemerides queries, TDB for element and vector queries. 

So while AstroPy is assuming UTC (you would be better passing your observing_time Time object to the SkyCoord routine's constructor so it would explicitly know what timescale it's in and be able to convert it), you are passing a Julian Date (JD) in the UTC timescale to HORIZONS which is expecting a time in the TDB timescale. At the time you have asked for, the difference TDB-UTC is ~67.184s (plus or minus 1-2 milliseconds for relativistic reasons we can ignore) so you are not asking for the position of the Earth at the same time via the two methods.

If you change your call to Horizons to obj = Horizons(id="Geocenter", location="@sun", epochs=observing_time.tdb.jd, id_type='id').vectors(delta_T=True) instead then epochs will get a correct JD in the TDB timescale corresponding to the same UTC instant. (I also included the delta_T=True in the call to .vectors() so you can see the exact value of TDB-UTC). With this change, the resulting HORIZONS position matches the Astropy one within ~2 km which is within the stated precision for the position of the Earth produced by the erfa.epv00 routine used internally by AstroPy.