[AstroPy] consistency of reference frames between SkyCoord and JPL Horizons
Michael Brewer
brewer at astro.umass.edu
Fri May 1 17:20:36 EDT 2020
Paolo,
I can tell you how to do this. First query JPL Horizons for your target
body setting the observer location to the Earth's geocenter. This will
return the astrometric position of the target body in the ICRF frame
referred to the Earth's geocenter including the proper light time
offset. Now you want these coordinates as an astropy SkyCoord object in
astropy's ICRS frame. The problem here is that astropy always refers
ICRS coordinates to the solar system barycenter. This will require you
to add the solar system barycentric position of the Earth's geocenter.
The way to do this is: First create a SkyCoord object using the
coordinates that you got from your query, call get_body_barycentric() to
get the Earth's position and then add these two together to get a proper
astropy SkyCoord object. This looks like this:
from astropy.coordinates.solar_system import get_body_barycentric
my_object_ICRS = SkyCoord( SkyCoord( JPL geocentric astrometric
coordinates ).cartesian + get_body_barycentric('earth', my_time) )
Now if you want to compare the result with JPL Horizons apparent place,
you need to realize that JPL's apparent place is referenced to the true
equator and equinox of date whereas astropy's CIRS frame is referenced
to the Celestial Intermediate Origin. Fortunately, there is a function
that converts GCRS coordinates to true apparent place:
from astropy.coordinates.solar_system import
_apparent_position_in_true_coordinates
my_object_apparent = _apparent_position_in_true_coordinates(
my_object_ICRS.transform_to('gcrs') )
The result should compare to JPL Horizons to 53 mas in RA. This 53 mas
is a rather mysterious offset that I pointed out to the JPL Horizons
people awhile ago. It took me awhile to convince them of this, but I'm
happy to see that they have now accepted it and are making a note of it
in their output files. There will be some small residual offsets due to
the fact that JPL Horizons uses a slightly different set of earth
orientation parameters than those used by astropy, but these are usually
less than 1 mas.
I have raised an issue in this regard on github (Issue #10230)
requesting two new built-in frames for converting to geocentric or
topocentric astrometric place and true apparent place. If these ever get
implemented, this process will become a lot easier.
Best Regards,
Michael Brewer
On 05/01/2020 04:01 AM, Paolo Tanga wrote:
>
> Hi Eric
>
> thank you for your other comments. Some answers below.y
>
> On 01/05/2020 05:09, Eric Jensen wrote:
>> Like Ben, I’m a little unclear on exactly what you are trying to
>> achieve beyond what you’ve demonstrated already. (I have learned
>> from your example!) You have shown that Horizons provides
>> astrometric coordinates that are ICRS, and that these are
>> straightforward to insert into a SkyCoord object as you show in your
>> example.
>
> my whole point here is that this works *only* when you query Horizons
> for astrometric coordinates with an observer at "@0" (= the Solar
> System barycenter), and this is not explicitly explained in the
> documentation. My other point is that I am not so expert in all
> features in astropy, so I wanted to have a confirmation (which
> essentially I got, I think). Eventually, somebody could be able to
> define a "custom" frame that works for all output returned by Horizons.
>
>>
>> If for some reason you want to use the apparent coordinates instead,
>> then you need to put them into a SkyCoord frame that precesses, which
>> as I understand it is what the PrecessedGeocentric frame is for. So
>> with a slight modification of your example, you can do this:
>>
>> (...)Doing this, I get coordinates that differ by about 13 arcsec
>> (mostly in RA) from the GCRS coords you create in the first part of
>> your example using the ICRS coords at the barycenter. I don’t know
>> what the source of that difference is, but perhaps some of the
>> various GR corrections.
>
> Thank you for this additional example, but in my view this shows that
> even in this approach you do not get consistent results. 13 arcsec are
> a huge discrepancy, not for most amateur use, but for today's
> standard. In fact, probably the effect is due to annual aberration
> (most of it) and light bending (a few mas). But what is the interest
> of including relativistic corrections in astropy, if then you cannot
> get an agreement at ~mas level, at least? For professional use in the
> Gaia era, this is a requirement.
>
>> There are indeed subtleties in converting between various coordinates
>> systems, as you note, but I would say that those subtleties are
>> inherent to the topic in general, not to the astropy software in
>> particular.
> True! And the whole topic is complex. But astropy already does a great
> service to the user, by providing out of the box, accurate
> transformations with all the bells and whistles. So, providing
> information on how to use (at full accuracy!) astroquery.horizons (or
> new reference frames to improve consistency!) I think goes in the same
> direction. By the way, note that the issues in the other mail thread
> by Michael Brewer also partially addresses the same problem.
>> That said, examples are quite helpful, so I’d encourage you to submit
>> an example based on your work here for inclusion on the examples page
>> at https://docs.astropy.org/en/stable/generated/examples/ .
>
> I will certainly do, including warnings about how NOT to use the
> Horizons output with SkyCoords... ;-) ...!
>
> Cheers
>
> Paolo
>
>
>>
>>
>>> On Apr 29, 2020, at 11:23 AM, Paolo Tanga <Paolo.Tanga at oca.eu
>>> <mailto:Paolo.Tanga at oca.eu>> wrote:
>>>
>>> Hi Ben
>>>
>>> thank you for taking your time to look into this question, and sorry
>>> for the delay.
>>>
>>> I would say, all that I am trying to do is to use the astroquery
>>> output in a SkyCoord object coherently. If one manages to do that,
>>> than it is possible to take profit of all the SkyCoord methods of
>>> course, to manipulate the coordinates (transform_to ) etc.
>>>
>>> As the SkyCoord is a fundamental part of the astropy suite, I would
>>> expect that there should be a "easy" method (or at least, a
>>> documented one!) to pass from an astroquery/JPL result to SkyCoord.
>>> But it turns out that this is not the case.
>>>
>>> So, what we have learned so far, if I am correct, is that:
>>>
>>> - JPL apparent coordinates. Precession - you are totally right -
>>> creates the difference between GCRS and the JPL apparent
>>> coordinates. In the documentation I would explicitly write that
>>> these cannot easily be used in a SkyCoord frame. In fact I tried
>>> playing with different reference frames (including
>>> GeocentricTrueEquator etc) but I see the same inconsistencies as
>>> with GCRS. So, I suppose that if there is a method, it is not
>>> straightforward (at least, not for me).
>>>
>>> - JPL astrometric coordinates. If computed with the observer at the
>>> barycenter of the Solar System, they can fit a SkyCoord object in
>>> ICRS. And this is the only possibility.
>>>
>>> I think that there could be other subtleties to check, but don't you
>>> think that this should be documented somewhere?
>>>
>>> Note that in the case of the ephemeris of the major planets
>>> (get_body()...) a SkyCoord object is directly returned, which is
>>> very convenient.
>>>
>>> Cheers
>>>
>>> Paolo
>>>
>>>
>>> On 25/04/2020 10:43, Benjamin Weiner wrote:
>>>> Hello Paolo,
>>>>
>>>> I am not entirely sure what end result you wish to achieve, nor am
>>>> I an expert on SkyCoord. However, I have a suggestion to
>>>> understand further these discrepancies (and thank you for including
>>>> the notebook example). In your first example:
>>>> (a) JPL's barycentric coords to build ICRS at specified time, then
>>>> transformed to GCRS at specified time.
>>>> (b) JPL's apparent coords to build GCRS at specified time.
>>>>
>>>> The separation of 16 arcmin between (a) and (b) is due to the
>>>> precession of the Earth's axis from 2000 to 2020. This can be seen
>>>> if we substitute
>>>> time = Time('2000-01-01T00:00:00', scale='utc')
>>>> in the line of your notebook that sets the time. At 2000-01-01, the
>>>> separation reduces to 8 arcseconds, which you have identified as
>>>> coming from the aberration of starlight.
>>>>
>>>> I think that the discrepancy is in building the GCRS coordinates
>>>> from the apparent RA/Dec. The JPL apparent coordinates are
>>>> referenced to the position of the Earth's axis at the specified
>>>> time. But the axes of GCRS are non-rotating with respect to the
>>>> axes of BCRS and ICRS - that is, GCRS is fixed to the Earth's
>>>> center, but not to the Earth's axis, and does not precess. I
>>>> believe that this is correct - USNO circular 179 has more gory
>>>> details on the definition of coordinate frames, see
>>>> https://www.usno.navy.mil/USNO/astronomical-applications/publications/Circular_179.pdf
>>>>
>>>> Cheers,
>>>> Ben
>>>>
>>>> On Fri, Apr 24, 2020 at 9:01 AM <astropy-request at python.org
>>>> <mailto:astropy-request at python.org>> wrote:
>>>>
>>>>
>>>> Message: 1
>>>> Date: Fri, 24 Apr 2020 15:00:22 +0200
>>>> From: Paolo Tanga <Paolo.Tanga at oca.eu <mailto:Paolo.Tanga at oca.eu>>
>>>> To: AstroPy at python.org <mailto:AstroPy at python.org>
>>>> Subject: [AstroPy] consistency of reference frames between SkyCoord
>>>> and JPL Horizons
>>>> Message-ID: <b18102bb-0e8c-8b85-4e66-8de7aecab2e4 at oca.eu
>>>> <mailto:b18102bb-0e8c-8b85-4e66-8de7aecab2e4 at oca.eu>>
>>>> Content-Type: text/plain; charset="utf-8"; Format="flowed"
>>>>
>>>> Hi all
>>>>
>>>> I am trying to plug the results by JPL Horizons into the
>>>> appropriate
>>>> reference system of a SkyCoord object, but it looks very trick.
>>>> (bewore
>>>> long message, I hope the explanation is clear for the
>>>> specialists of
>>>> reference systems)
>>>>
>>>> The Horizons server returns two types of equatorial coordinates:
>>>> "astrometric" (1) and (2) "apparent".
>>>>
>>>> Their description is the following:
>>>>
>>>> - (1): astrometric RA and Dec with respect to the observing site
>>>> (coordinate origin) in the reference frame of
>>>> the planetary ephemeris (ICRF). Compensated for down-leg
>>>> light-time
>>>> delay aberration.
>>>>
>>>> - (2): Airless apparent RA and Dec of the target with respect
>>>> to an
>>>> instantaneous reference frame defined by the Earth equator of-date
>>>> (z-axis) and meridian containing the Earth equinox of-date
>>>> (x-axis,
>>>> IAU76/80). Compensated for down-leg light-time delay,
>>>> gravitational
>>>> deflection of light, stellar aberration, precession & nutation.
>>>> Note:
>>>> equinox (RA origin) is offset -53 mas from the of-date frame
>>>> defined by
>>>> the IAU06/00a P & N system.
>>>>
>>>> My guess is none of these two systems correspond to either ICRS
>>>> or GCRS
>>>> as defined in astropy. The simplest possibility for me, would
>>>> have been
>>>> to assimilate (1), computed for the Solar System barycenter, to
>>>> ICRS.
>>>> But definitely (2) is not GCRS and the test below easily proves
>>>> that.
>>>>
>>>> Verification (see the attached notebook if you want to play with) :
>>>>
>>>> - I queried for Solar System barycentric coordinates (1) of an
>>>> asteroid
>>>> at one epoch. I inserted them in a SkyCoord object as ICRS and
>>>> use the
>>>> transform_to(GCRS) towards geocentric, at the given epoch.
>>>>
>>>> - I queried again JPL Horizons for the apparent coordinates (2)
>>>> from the
>>>> geocenter, directly.
>>>>
>>>> As expected, by comparing the results I get the patent
>>>> confirmation that
>>>> GCRS is not the frame of (2), with a discrepancy of 16 arcmin
>>>> in my
>>>> case. SO, IT DOES NOT WORK, AS EXPECTED. Let's put apparent
>>>> coordinates
>>>> (2) aside for the moment.
>>>>
>>>> Second possibility (see notebook again).
>>>>
>>>> I take coordinates (1), for the geocentric observer. I insert
>>>> them as
>>>> GCRS in a SkyCoord object. I trasform_to(ICRS) (barycentric
>>>> observer)
>>>> and compare the result to a query of (1) directly for the
>>>> barycentric
>>>> observer.
>>>>
>>>> Again this SHOULD NOT WORK accurately. In fact the definition
>>>> of (1)
>>>> does not contain annual aberration and light deflection that
>>>> are part of
>>>> GCRS. A difference of several arcsec (mostly due to annual
>>>> aberration)
>>>> SHOULD show up. And in fact IT DOES, ~8 arcsec of difference.
>>>>
>>>> The conclusions if that using SkyCoord with the output of JPL
>>>> Horizons
>>>> is far from obvious to the end user...
>>>>
>>>> ----
>>>>
>>>> So, my final questions/remarks are:
>>>>
>>>> - How to get a fully accurate consistency between a SkyCoord
>>>> frame and
>>>> JPL Horizons output, in the case of a barycentric, geocentric or
>>>> topocentric observer, for both "astrometric" and "apparent"
>>>> coordinates
>>>> as provided by JPL? Which reference SkyCoord should use and how?
>>>>
>>>> - Does astropy implement GCRS correctly in the case of Solar
>>>> System
>>>> objects, by taking into account aberration and light deflection
>>>> in the
>>>> relativistic form, for an object at finite distance (different
>>>> than for
>>>> the stars!) ?
>>>>
>>>> - ...or... Am I missing something fundamental in my interpretation?
>>>>
>>>> - Eventually, this query to Horizons turns our to be tricky to
>>>> manage as
>>>> SkyCoord for the final user. I think it would be useful to
>>>> provide the
>>>> correct recipe and/or directly a SkyCoord as output of the query...
>>>>
>>>> Best regards
>>>>
>>>> Paolo
>>>>
>>>>
>>>> --
>>>> Benjamin Weiner
>>>> Staff Scientist / Associate Astronomer, MMT / Steward Observatory
>>>> bjw at as.arizona.edu <mailto:bjw at as.arizona.edu>
>>>> http://mingus.mmto.arizona.edu/~bjw/
>>>> <http://mingus.as.arizona.edu/%7Ebjw/>
>>>>
>>>> _______________________________________________
>>>> AstroPy mailing list
>>>> AstroPy at python.org
>>>> https://mail.python.org/mailman/listinfo/astropy
>>> --
>>> ---------------------------------------------------------------------
>>> Paolo Tanga Astronomer
>>> Deputy director of Laboratoire Langrange / UMR 7293
>>> Observatoire de la Côte d'Azur Tel +33(0)492003042
>>> Bv de l'Observatoire - CS 34229 Fax +33(0)492003121
>>> 06304 Nice Cedex 4 - Francehttp://www.oca.eu/tanga
>>> https://twitter.com/ziggypao
>>> _______________________________________________
>>> AstroPy mailing list
>>> AstroPy at python.org <mailto:AstroPy at python.org>
>>> https://mail.python.org/mailman/listinfo/astropy
>>
> --
> ---------------------------------------------------------------------
> Paolo Tanga Astronomer
> Deputy director of Laboratoire Langrange / UMR 7293
> Observatoire de la Côte d'Azur Tel +33(0)492003042
> Bv de l'Observatoire - CS 34229 Fax +33(0)492003121
> 06304 Nice Cedex 4 - Francehttp://www.oca.eu/tanga
> https://twitter.com/ziggypao
>
>
> _______________________________________________
> AstroPy mailing list
> AstroPy at python.org
> https://mail.python.org/mailman/listinfo/astropy
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/astropy/attachments/20200501/763bf7e1/attachment-0001.html>
More information about the AstroPy
mailing list