[AstroPy] consistency of reference frames between SkyCoord and JPL Horizons
Michael Brewer
brewer at astro.umass.edu
Fri May 1 18:09:29 EDT 2020
Oh, I forgot. This accuracy assumes that you are using the same DE430
ephemeris that JPL Horizon uses.
solar_system_ephemeris.set('de430')
On 05/01/2020 05:20 PM, Michael Brewer wrote:
> 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
>
>
>
> _______________________________________________
> 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/80151e30/attachment-0001.html>
More information about the AstroPy
mailing list