keflavich
commented
3 years ago @mkelley can you help with this?
Closed nkelhos closed 3 years ago
keflavich
commented
3 years ago @mkelley can you help with this?
mkelley
commented
3 years ago Hi @nkelhos , one way to get more information on the results is to examine the raw response from Horizons.
from astroquery.jplhorizons import Horizons
from astropy.time import Time
obj = Horizons( id='6', id_type='majorbody', location=None, epochs=Time(['2021-01-01']).jd)
raw = obj.vectors(get_raw_response=True)
print(raw)*******************************************************************************
Revised: April 12, 2021 Saturn Barycenter 6
Dynamical point:
---------------
The location of the center-of-mass of the Saturn system. The point about which
Saturn and its satellites revolve. See 699 for Saturn center.
*******************************************************************************
*******************************************************************************
Ephemeris / WWW_USER Tue Aug 3 07:54:21 2021 Pasadena, USA / Horizons
*******************************************************************************
Target body name: Saturn Barycenter (6) {source: DE441}
Center body name: Sun (10) {source: DE441}
Center-site name: BODY CENTER
*******************************************************************************
Start time : A.D. 2021-Jan-01 00:00:00.0000 TDB
Stop time : A.D. 2021-Jan-01 00:00:00.5000 TDB
Step-size : 0 steps
*******************************************************************************
Center geodetic : 0.00000000,0.00000000,0.0000000 {E-lon(deg),Lat(deg),Alt(km)}
Center cylindric: 0.00000000,0.00000000,0.0000000 {E-lon(deg),Dxy(km),Dz(km)}
Center radii : 696000.0 x 696000.0 x 696000.0 k{Equator, meridian, pole}
Output units : AU-D
Output type : GEOMETRIC cartesian states
Output format : 3 (position, velocity, LT, range, range-rate)
Reference frame : Ecliptic of J2000.0
*******************************************************************************
JDTDB, Calendar Date (TDB), X, Y, Z, VX, VY, VZ, LT, RG, RR,
**************************************************************************************************************************************************************************************************************************************************************************
$$SOE
2459215.500000000, A.D. 2021-Jan-01 00:00:00.0000, 5.490358280342060E+00, -8.342122639745645E+00, -7.347716850247014E-02, 4.356569070507491E-03, 3.057983937039585E-03, -2.265013009596026E-04, 5.768018082597512E-02, 9.987013721697220E+00, -1.576356389034122E-04,
$$EOE
**************************************************************************************************************************************************************************************************************************************************************************
Coordinate system description:
Ecliptic at the standard reference epoch
Reference epoch: J2000.0
X-Y plane: adopted Earth orbital plane at the reference epoch
Note: obliquity of 84381.448 arcseconds (IAU76) wrt ICRF equator
X-axis : ICRF
Z-axis : perpendicular to the X-Y plane in the directional (+ or -) sense
of Earth's north pole at the reference epoch.
Symbol meaning [1 au= 149597870.700 km, 1 day= 86400.0 s]:
JDTDB Julian Day Number, Barycentric Dynamical Time
X X-component of position vector (au)
Y Y-component of position vector (au)
Z Z-component of position vector (au)
VX X-component of velocity vector (au/day)
VY Y-component of velocity vector (au/day)
VZ Z-component of velocity vector (au/day)
LT One-way down-leg Newtonian light-time (day)
RG Range; distance from coordinate center (au)
RR Range-rate; radial velocity wrt coord. center (au/day)
Geometric states/elements have no aberrations applied.
Computations by ...
Solar System Dynamics Group, Horizons On-Line Ephemeris System
4800 Oak Grove Drive, Jet Propulsion Laboratory
Pasadena, CA 91109 USA
Information : https://ssd.jpl.nasa.gov/
Documentation: https://ssd.jpl.nasa.gov/?horizons_doc
Connect : https://ssd.jpl.nasa.gov/?horizons (browser)
telnet ssd.jpl.nasa.gov 6775 (command-line)
e-mail command interface available
Script and CGI interfaces available
Author : Jon.D.Giorgini@jpl.nasa.gov
*******************************************************************************
!$$SOF
TABLE_TYPE = VECTORS
OUT_UNITS = AU-D
COMMAND = "6"
CENTER = '500@10'
CSV_FORMAT = "YES"
REF_PLANE = ECLIPTIC
REF_SYSTEM = J2000
TP_TYPE = ABSOLUTE
LABELS = YES
VECT_CORR = "NONE"
VEC_DELTA_T = NO
OBJ_DATA = YES
TLIST = 2459215.5It does appear to be using the Sun and an obliquity of the ecliptic of 84381.448 arcseconds. See astropy/astropy#6461 for some notes on that reference plane. Does this address your problem?
BTW, I thought that there might instead be a problem with the time scale. In my example above, Horizons took my time input to be in the TDB scale. But, if one was expecting UTC (as is the default for astropy Time objects), the difference is quite small, 0.01 Saturn radii:
obj_utc_input = Horizons( id='6', id_type='majorbody', location=None, epochs=Time(['2021-01-01']).jd)
result_utc_input = obj_utc_input.vectors()
obj_tdb_input = Horizons( id='6', id_type='majorbody', location=None, epochs=Time(['2021-01-01']).tdb.jd)
result_tdb_input = obj_tdb_input.vectors()
d = np.sqrt((result_tdb_input['x'] - result_utc_input['x'])**2 + (result_tdb_input['y'] - result_utc_input['y'])**2 + (result_tdb_input['z'] - result_utc_input['z'])**2).data
print(d * u.au.to('km') / 58232)The documentation on the time scale not correct and I will fix that now: https://github.com/astropy/astroquery/blob/32b0af5e37775af09d1f780f1ecdb6759215b933/astroquery/jplhorizons/core.py#L59
nkelhos
commented
3 years ago @mkelley, thanks for pointing out raw_response and the tdb time difference. I think I'm now pretty convinced jpl uses the barycenter of the solarsystem. However, from this:
Ecliptic at the standard reference epoch
Reference epoch: J2000.0
X-Y plane: adopted Earth orbital plane at the reference epoch
Note: obliquity of 84381.448 arcseconds (IAU76) wrt ICRF equator
X-axis : ICRF
Z-axis : perpendicular to the X-Y plane in the directional (+ or -) sense
of Earth's north pole at the reference epoch.its not quite clear if they use the Mean ecliptic or True ecliptic. Wikipedia points out:
Use of the "mean" locations means that nutation is averaged out or omitted.On a whim, I went and plotted the differences between the astropy and jpl coordinates. This resulted in a few plots, but generally:
The above plot is the XY coordinates, and the below plot is the YZ coordinates.
The orange arrows change direction slowly throughout one saturn year. I went and plotted the position of jupiter, and jupiter passing by seems to coincide with the strong inward-pointing around 2001 and 2019.
The outward-pointing differences around 2011 seem odd, maybe due to an outer planet, so I plotted more planets:
but I realized its probably because jupiter is furthest from saturn at that time.
After the fact, it seems kind of obvious that nearby planets will squeeze orbits, but it seems odd the perturbations are handled differently between the two coordinates.
What struck me as weird is the orange arrows point from jpl to astropy, so it seems like de432s has stronger perturbations than de441, even though its 7 years older. Maybe they both modelled the perturbations, but de432 used a cruder model or something.
These plots were made with a python script, which I've included here if anyone wants to see how I did this. planets_for_github.log (change extension to .py)
I guess I'm going to try getting de440 or de441 loaded into astropy, to see if that reduces the differences between astropy and jpl, to make sure the coordinate systems are converting properly.
mkelley
commented
3 years ago Interesting plots, @nkelhos .
@mkelley, thanks for pointing out raw_response. I think I'm now pretty convinced jpl uses the barycenter of the solarsystem.
But that isn't consistent with the Horizons output for location=None, which contains:
Center body name: Sun (10) {source: DE441}
Center-site name: BODY CENTERUsing location='500@0' gets the barycenter:
Center body name: Solar System Barycenter (0) {source: DE441}
Center-site name: BODY CENTERI did this once for the Sun and once for the SSB and the vectors are different:
# Sun
2459215.500000000, A.D. 2021-Jan-01 00:00:00.0000, 5.490358280342060E+00, -8.342122639745645E+00, -7.347716850247014E-02, 4.356569070507491E-03, 3.057983937039585E-03, -2.265013009596026E-04, 5.768018082597512E-02, 9.987013721697220E+00, -1.576356389034122E-04,
# Solar System Barycenter
2459215.500000000, A.D. 2021-Jan-01 00:00:00.0000, 5.483706347577082E+00, -8.336119049840741E+00, -7.337139832615820E-02, 4.349729206314884E-03, 3.052221555062912E-03, -2.262854082984947E-04, 5.763009436375799E-02, 9.978341519594688E+00, -1.577859116832880E-04,The difference between the two is of order 0.01 au == 26 Saturn radii. Does using '500@0' solve the problem? I ran your code with this and it looks a lot better:
# location=None
difference vectors:
minimum : 4.672e+08 meters = 4.01 saturn diameters
mean : 6.734e+08 meters = 5.78 saturn diameters
standev : 1.118e+08 meters = 0.96 saturn diameters
maximum : 8.103e+08 meters = 6.96 saturn diameters
# location='500@0'
difference vectors:
minimum : 9.674e+07 meters = 0.83 saturn diameters
mean : 9.853e+07 meters = 0.85 saturn diameters
standev : 1.129e+06 meters = 0.01 saturn diameters
maximum : 1.003e+08 meters = 0.86 saturn diametersBut that isn't consistent with the Horizons output for location=NoneOof, thanks for catching that @mkelley . I think I didn't read the docpage here closely enough:
The default is location=None, which uses a geocentric location for ephemerides queries and the Sun as central body for orbital elements and state vector queries.The Horizons web browser describes more location choices, and it looks like the 500 part isn't necessary, unless I want RA/Dec from a particular observatory. Using @0 doesn't seem to have any effect. In retrospect, the 'perturbations from Jupiter' is probably mathematically similar to movement in the solar system barycenter, since Jupiter is the second-largest mass in the area.
The location=@0 plot looks much better:
Even with the arrows still being scaled by 100x, they are now much smaller.
However, all the arrows now point along the orbital path which is suspicious in its own right. This could be a timing issue, and since they're all the same length, it could be due to using a wrong equinox. When I switch astropy to use BarycentricMeanEcliptic (instead of BarycentricTrueEcliptic), the difference between jpl and astropy shrinks from 0.85 saturn diameters to 0.17 Saturn diameters:
location='@0', barycentrictrueecliptic
difference vectors:
minimum : 9.326e+07 meters = 0.80 saturn diameters
mean : 9.858e+07 meters = 0.85 saturn diameters
standev : 3.688e+06 meters = 0.03 saturn diameters
maximum : 1.037e+08 meters = 0.89 saturn diameters
location='@0', barycentricmeanecliptic
difference vectors:
minimum : 1.918e+07 meters = 0.16 saturn diameters
mean : 2.035e+07 meters = 0.17 saturn diameters
standev : 8.066e+05 meters = 0.01 saturn diameters
maximum : 2.147e+07 meters = 0.18 saturn diametersNow when I look at the difference vectors, they're small enough they need to be scaled by 10000x, and you can see they're all pointing inwards now:
I'm not sure what this could be caused by. I'm vaguely aware that in ephemeris observations, the RA/Dec error is usually much smaller than the range error. Updating to de440 might improve this, since it has 7 more years of observations than de432s.
nkelhos
commented
3 years ago I've expanded the plots to include all the planets, and you can see the orange difference arrows between astropy(de440) and jpl(de440) point almost directly at the origin.
I've also noticed that the length of the orange difference vectors scale (by a fixed constant) based on the orbital radius of the planet, which is also suspicious.
A) Mean Orange Difference Vector Lengths:
mercury : 8.4352e+05 meters = 0.1729 mercury diameters
venus : 1.5398e+06 meters = 0.1272 venus diameters
earth : 2.1307e+06 meters = 0.1672 earth diameters
mars : 3.2602e+06 meters = 0.4809 mars diameters
jupiter : 1.1075e+07 meters = 0.0792 jupiter diameters
saturn : 2.0360e+07 meters = 0.1748 saturn diameters
uranus : 4.0906e+07 meters = 0.8064 uranus diameters
neptune : 6.4039e+07 meters = 1.3004 neptune diameters
pluto : 8.6706e+07 meters = 36.4834 pluto diametersB) Mean Planet Orbit Radius
mercury : 5.9260628e+10 meters
venus : 1.0817808e+11 meters
earth : 1.4968840e+11 meters
mars : 2.2904337e+11 meters
jupiter : 7.7804092e+11 meters
saturn : 1.4303519e+12 meters
uranus : 2.8738029e+12 meters
neptune : 4.4989899e+12 meters
pluto : 6.0914655e+12 metersRatio A/B :
mercury : 1.4234027e-05
venus : 1.4234031e-05
earth : 1.4234021e-05
mars : 1.4234027e-05
jupiter : 1.4234017e-05
saturn : 1.4234044e-05
uranus : 1.4234022e-05
neptune : 1.4234029e-05
pluto : 1.4234094e-05The constant is 1.42340e-05. Another way to say this is that, you could tell me the radius of some planet, and by multiplying the radius by that constant, I could tell you how far apart astropy and jpl are, which means something predictable is wrong.
jpl's Horizons().vectors() likes to return values in AU, and I believe de440 also stores distances in AU, so I was thinking of checking the AU-to-meter conversion I was using. In my script I was using code like:
AU2meters = 1.496e11
obj = Horizons( id='6', id_type='majorbody', location='@0', epochs=times.tdb.jd )
jpl = obj.vectors()
jpl_x = list( [ x * AU2meters for x in jpl['x'] ] ) # meters
jpl_y = list( [ y * AU2meters for y in jpl['y'] ] ) # meters
jpl_z = list( [ z * AU2meters for z in jpl['z'] ] ) # meters so I tried using the more precise 1.49597870700e11, but that didnt affect anything.
Then I tried using the astropy.units module as much as possible:
AU2meters = 1.49597870700e11 # astropy 4.2 code (iau2015.py)
obj = Horizons( id='6', id_type='majorbody', location='@0', epochs=times.tdb.jd )
jpl = obj.vectors()
#jpl_x = list( [ x * AU2meters for x in jpl['x'] ] ) # off by a large factor
#jpl_y = list( [ y * AU2meters for y in jpl['y'] ] ) # off by a large factor
#jpl_z = list( [ z * AU2meters for z in jpl['z'] ] ) # off by a large factor
jpl_x = list( ( [ x for x in jpl['x'] ] * units.au).to_value(units.m) ) # 100x improvement!
jpl_y = list( ( [ y for y in jpl['y'] ] * units.au).to_value(units.m) ) # 100x improvement!
jpl_z = list( ( [ z for z in jpl['z'] ] * units.au).to_value(units.m) ) # 100x improvement!and with that, the difference vectors shrank by a factor of 100x, which is great!
mercury : 7.9603e+03 meters = 0.0016 mercury diameters
venus : 1.4598e+04 meters = 0.0012 venus diameters
earth : 1.9893e+04 meters = 0.0016 earth diameters
mars : 3.0642e+04 meters = 0.0045 mars diameters
jupiter : 1.0325e+05 meters = 0.0007 jupiter diameters
saturn : 1.9534e+05 meters = 0.0017 saturn diameters
uranus : 3.8294e+05 meters = 0.0075 uranus diameters
neptune : 6.0469e+05 meters = 0.0123 neptune diameters
pluto : 8.9051e+05 meters = 0.3747 pluto diameterse.g. saturn went from differences of 0.17 diameters to 0.0017 diameters. We're now at the point where the astropy and jpl barycenter coordinates are only different by 8-891 km.
However, this makes me pretty uneasy, since going from AU to meters should just be multiplying by one constant, but when I do that manually, I get a different result. It could be that both astropy and jpl now match better because they've both had the same wrong conversion applied within the python machinery of the script :( .
I've attached the latest version of the script here: for_github_v2.log (change extension to .py)
mkelley
commented
3 years ago Great, it looks like astroquery is working as it should be, so perhaps this issue can be closed? Or, perhaps we can update the documentation and name the astropy reference frame that the Horizons vectors are returned in?
so I tried using the more precise 1.49597870700e11, but that didnt affect anything.
That sounds very strange, so I ran your updated script once with lines 101-103 enabled, and again with lines 98-100 enabled. But the results were the same.
nkelhos
commented
3 years ago The jpl coords seem to match the ones in the get_raw_response text, but I'm still investigating whats going on.
perhaps we can update the documentation and name the astropy reference frame that the Horizons vectors are returned in?
I'd planned on using this conversion in a project, but if it would be appropriate, I can also add this as a possible output of Horizons().vectors(). Maybe have a separate kwarg for returning sky coordinates, like vectors(get_skycoords=True) and have it return an astropy.coordinates.SkyCoord() object with the coordinates already in the right system. I'm not quite sure if that fits into astroquery's design schema, though it does look like astroquery already depends on astropy, so it might be ok?
nkelhos
commented
3 years ago I ran your updated script once with lines 101-103 enabled, and again with lines 98-100 enabled. But the results were the same.
Yes, I get the same behavior. I think I was wrong, the improvement was from the AU2meters. With lines 98-100 enabled, and setting
AU2meters = 1.496e11 # meters/AUyou get a difference of 0.17 saturn diameters. With
AU2meters = 1.49597870700e11 # astropy 4.2 code (iau2015.py)you get 0.0017 saturn diameters.
Now that this is fixed, the difference vectors have shrunk, and they can be seen below:
I was thinking this was another equinox problem, but mercury and pluto's vectors are asymmetric, unlike the other planets, so I'm not sure. I'm noticing in the raw response that jpl is using de441, and from the publication on de440 & de441's release:
Therefore, DE441 is less accurate than DE440 for the current century, but covers a much longer duration of years −13,200 to +17,191, compared to DE440 covering years 1550–2650.
So my next step is probably to try switching to de441, to see if that makes up the difference.
mkelley
commented
3 years ago I do want to see the output use SkyCoord (and Time) objects. It is a significant change in the output, and using a keyword to control the output is OK, but adds complexity, so I'm more in favor of putting the SkyCoord in a separate column or replacing the RA, Dec columns altogether. @keflavich do you have any recommendations?
nkelhos
commented
3 years ago The angles of the difference vectors are aligned along the orbit, and while they're different lengths for different planets, they're all roughly the same 7.6e-6 deg from the solar system barycenter.
Mean Difference Vectors:
mercury 7.9603e+03 meters = 0.0016 mercury diameters = 7.69652e-06 deg
venus 1.4598e+04 meters = 0.0012 venus diameters = 7.73162e-06 deg
earth 1.9893e+04 meters = 0.0016 earth diameters = 7.61445e-06 deg
mars 3.0642e+04 meters = 0.0045 mars diameters = 7.66519e-06 deg
jupiter 1.0325e+05 meters = 0.0007 jupiter diameters = 7.60387e-06 deg
saturn 1.9534e+05 meters = 0.0017 saturn diameters = 7.82468e-06 deg
uranus 3.8294e+05 meters = 0.0075 uranus diameters = 7.63481e-06 deg
neptune 6.0469e+05 meters = 0.0123 neptune diameters = 7.70104e-06 deg
pluto 8.9051e+05 meters = 0.3747 pluto diameters = 8.37621e-06 degI tried switching to de441_part-2.bsp, but it offered no improvement. It might be a difference in how many significant figures were used in the equinox specification, or just a slightly different mean equinox time.
mkelley
commented
3 years ago After referring to astropy/astropy#8394 and several related threads...
It looks like you can get past that last 30 mas by using the CustomBarycentricEcliptic frame and giving it the same obliquity that Horizons uses:
frame = CustomBarycentricEcliptic(obliquity=84381.448 * units.arcsec)
apy_s = apy_sc.transform_to(frame)I also got good results for a heliocentric computation, but in this case, it is important to use DE440 and the observation times for the Sun's ephemeris by promoting a few lines into the solar_system_ephemeris state context:
with solar_system_ephemeris.set('de440.bsp'):
apy = get_body_barycentric(pastr, times)
apy_sc = SkyCoord(apy, representation_type='cartesian', frame='icrs', obstime=times)
apy_s = apy_sc.transform_to('heliocentriceclipticiau76')
bsipocz
commented
3 years ago so I'm more in favor of putting the SkyCoord in a separate column or replacing the RA, Dec columns altogether.
Replacing the columns sounds the most favourable to me. In fact doing something like this is on our wishlist for a while. I suppose we can provide a few documented examples of how to get the previous raw strings out of such a column to help users switch.
nkelhos
commented
3 years ago I switched to:
with solar_system_ephemeris.set('de440.bsp') :
apy = get_body_barycentric( planet, times )
apy_sc = SkyCoord( apy, representation_type='cartesian', frame='icrs' )
frame = CustomBarycentricEcliptic( obliquity=84381.448 * units.arcsec )
apy_s = apy_sc.transform_to(frame) Promoting the bottom three lines didn't change anything in the following differences. The difference between the astropy and jpl coordinate systems is down to 10s of cm over each planets' year, and 10^10 deg from the origin:
mercury : 49.828 cm : 4.818e-10 deg
venus : 36.720 cm : 1.945e-10 deg
earth : 30.286 cm : 1.159e-10 deg
mars : 24.376 cm : 6.098e-11 deg
jupiter : 12.810 cm : 9.434e-12 deg
saturn : 10.321 cm : 4.134e-12 deg
uranus : 6.447 cm : 1.285e-12 deg
neptune : 5.777 cm : 7.357e-13 deg
pluto : 4.783 cm : 4.499e-13 degThe difference vectors now oscillate forward and backwards along the orbital paths:
mkelley
commented
3 years ago Promoting the bottom three lines didn't change anything in the following differences. The difference between the astropy and jpl coordinate systems is down to 10s of cm over each planets' year, and 10^10 deg from the origin:
Right, that shouldn't make a difference here. Sorry if I confused you. I was trying to figure out how to get heliocentric vectors. That was more about making notes for myself.
Now that the difference is down to the micro-arcsec level, I think the issue is nearly resolved. I'll open a PR that documents the frames to use, and that should finish it.
Thanks for the discussion and work here. I definitely learned a few things!
I'm trying to get the coordinates of an asteroid via astroquery.jplhorizons, and load them into astropy, but I'm a little uncertain about the coordinate system jplhorizons outputs.
If I get an object's xyz coordinates via:
The astroquery description for these x, y, and z coordinates is shown here.
My question is: What coordinate system are these x, y, and z values in?
In the astroquery doc page above, I can see references to "REF_PLANE" as ecliptic, and with
location=Nonesetting it to use the sun it might start using the sun's ecliptic plane, but I'm not sure if its using a True or Mean ecliptic, or if its sun-geocentric, sun-barycentric, or sun-jupiter-barycentric, solar-system-barycentric, or what.Astropy has a large number of coordinate systems. Through comparison with Saturn's coordinates in astropy and jplhorizons, I've narrowed it down to something like BarycentricTrueEcliptic or BarycentricMeanEcliptic.
I think BarycentricMean/TrueEcliptic coordinate systems use the barycenter of the entire solar system, but on the astroquery.jplhorizons page it says
location=None, which uses the Sun as central body for orbital elements and state vector queries.Neither BarycentricTrueEcliptic or BarycentricMeanEcliptic systems match up perfectly.
de430orjplsolar system ephemerides didn't change the errors.In summary, BarycentricTrueEcliptic gets the position of saturn to within 7.6 diameters, but I should be getting within 2.25 diameters, so something is still not quite right.
And, if there's an easier way to get jpl horizons asteroid coordinates into astropy, please let me know. Thanks for your help.