Not sure exactly what you want as the origin and what corrections for aberrations are wanted and whether you want Mercury itself or its barycenter but there didn't appear to be much/any difference in this case (no Mercurian moons...) but the following code defines a mapping for both and also for Parker Solar Probe (PSP) to the correct HORIZONS id. You can find these HORIZONS SPK ids out by typing the names into the web version of HORIZONS.
This example produces a daily position for 2020; change the start, end and ephem_step_size to taste but don't abuse the JPL service by asking for giant amounts of output calculations.
from datetime import datetime
from astropy.time import Time
from astroquery.jplhorizons import Horizons
object_id_map = { 'MercuryBody' : 199,
'MercuryBary' : 1,
'PSP' : -96
}
start = datetime(2020,1,1)
end = datetime(2021,1,1)
ephem_step_size='1d'
origin = '@ssb'
obj_name = 'MercuryBody'
# Setup HORIZONS query for the desired body by SPK id
eph = Horizons(id=object_id_map.get(obj_name, 1), id_type='id', epochs={'start' : start.strftime("%Y-%m-%d %H:%M"),
'stop' : end.strftime("%Y-%m-%d %H:%M"), 'step' : ephem_step_size}, location=origin)
# Query HORIZONS for the state vectors, set the reference plane to
# default 'ecliptic' (ecliptic and mean equinox of reference epoch),
# no aberrations (light travel time and stellar aberration not included),
# include time difference of TDB-UTC in output for possible further use (although
# astropy.Time can do this just as well)
table = eph.vectors(refplane='ecliptic', aberrations='geometric'. delta_T=True)
# Add a `datetime` column to the table for easier filtering. Needs to be done by making a
# Time array not a Column due to numpy issue (https://github.com/astropy/astropy/issues/9374) May not apply in Astropy 4.0+
dates = Time([datetime.strptime(d, "%Y-%b-%d %H:%M") for d in table['datetime_str']])
tables.add_column(dates, name='datetime')
print(table.colnames)
['targetname', 'datetime_jd', 'datetime_str', 'delta_T', 'x', 'y', 'z', 'vx', 'vy', 'vz', 'lighttime', 'range', 'range_rate']
x = table['x']
print(x.shape)
(367,)
The resulting table contains the X, Y, Z components of the vector (in AU) in the x, y, z columns and the velocities (in AU/day) in vz, vy, vz. AstroPy Columns are built on top of numpy arrays and provide the same functionality with additional functionality to handle Quantitys (Astropy Table documentation reference) as illustrated in the last line where we extract the X component.
