[AstroPy] consistency of reference frames between SkyCoord and JPL Horizons
Michael Brewer
brewer at astro.umass.edu
Fri May 8 17:55:42 EDT 2020
Hi Paolo,
I forgot about time. You need to tell the GCRS frame the time of the
observation. But you must have figured that out as you got good
agreement in position, so the error in distance is odd. I tried an
example with the following script:
from astropy.time import Time
from astropy.coordinates import get_body_barycentric, SkyCoord, GCRS
from astropy.coordinates.solar_system import
_apparent_position_in_true_coordinates
import astropy.units as u
#Mars Geocentric JPL 2020-May-07 00:00:00 328.364842462
-14.682567251 328.634953692 -14.587885357 1.18472191315873
my_time = Time('2020-05-07T00:00:00')
my_astrometric = SkyCoord( 328.364842462*u.deg, -14.682567251*u.deg,
1.18472191315873*u.au )
my_apparent = SkyCoord( 328.634953692*u.deg, -14.587885357*u.deg,
1.18472191315873*u.au )
print ('JPL ICRS', my_astrometric)
my_object_ICRS = SkyCoord( my_astrometric.cartesian +
get_body_barycentric('earth', my_time) )
print ('ICRS', my_object_ICRS)
my_object_apparent = _apparent_position_in_true_coordinates(
my_object_ICRS.transform_to(GCRS(obstime=my_time)))
print ('apparent', my_object_apparent)
print ('delta ra', (my_object_apparent.ra - my_apparent.ra).to(u.mas),
'delta dec', (my_object_apparent.dec - my_apparent.dec).to(u.mas),
'delta distance', (my_object_apparent.distance -
my_apparent.distance).to(u.meter) )
And got this result:
JPL ICRS <SkyCoord (ICRS): (ra, dec, distance) in (deg, deg, AU)
(328.36484246, -14.68256725, 1.18472191)>
ICRS <SkyCoord (ICRS): (ra, dec, distance) in (deg, deg, AU)
(282.36772087, -24.4176493, 1.42516501)>
apparent <SkyCoord (GCRS: obstime=2020-05-07T00:00:00.000,
obsgeoloc=(0., 0., 0.) m, obsgeovel=(0., 0., 0.) m / s): (ra, dec,
distance) in (deg, deg, AU)
(328.63496745, -14.58788562, 1.18472191)>
delta ra 49.5348mas delta dec -0.964467mas delta distance
-6.643480019441483e-05 m
The positions are off a bit more than I'd like to see, but the distance
is right on.
By the way, thanks to an example from Stuart Littlefair, I figured out
how to do this using topocentric coordinates also. This is just a matter
of setting up an EarthLocation my_loc and then doing:
earth_loc = get_body_barycentric('earth', my_time)
obsgeoloc, obsgeovel = my_loc.get_gcrs_posvel(my_time)
earth_loc += obsgeoloc
before adding earth_loc to your topocentric astrometric coordinates from
Horizons and then you also need to tell the GCRS frame your location and
velocity:
my_object_apparent = _apparent_position_in_true_coordinates(
my_object_ICRS.transform_to(
GCRS(obsgeoloc=obsgeoloc, obsgeovel=obsgeovel,
obstime=my_time)))
Regards,
Michael Brewer
On 05/07/2020 09:59 AM, Paolo Tanga wrote:
>
> Hi Micheal
>
> I did the test following your instructions and I correctly found the
> result that you indicate: 52.5 mas of discrepancty in RA, and
> negligible in Dec (0.01 mas).
>
> However, I notice something strange. At the first step:
>
> my_object_ICRS = SkyCoord( SkyCoord( JPL geocentric astrometric
> coordinates ).cartesian + get_body_barycentric('earth', my_time) )
>
> I give as "JPL geocentric astrometric coordinates" ra, dec and
> distance. I chose first the geocentric distance (delta, in my case
> 2.0958388 au).
>
> Then I get
>
> my_object_ICRS.distance = 3.816102 au.
>
> Formally, this may be correct, as this is the result of the addition
> of barycenter-to-Earth vector. But then, even if I am now in
> geocentric ICRS, it does not represent the geocentric distance.
>
> At the end of the process, after transformation to GCRS and then to
> apparent coordinates (by _apparent_position_in_true_coordinates) I find:
>
> my_object_apparent.distance = 2.962877... au, which is off the JPL
> value for delta given above by ~0.9 au.
>
> If at the beginning I give 'r' (heliocentric distance = 2.95757466,
> instead of 'delta'), I still get:
>
> my_object_ICRS.distance = 3.810857 but now my_object_apparent.distance
>
> but then:
>
> my_object_apparent.distance =2.957567
>
> which is much closer to the geocentric distance.... So, whatever the
> choice, it looks like the distance goes completely wrong and there is
> no way to obtain the correct value directly from the transformations.
> As ra and dec are correct, the only solution seems to apply a
> correction to the distance at the end of the computation.
>
> Is this behavior consistent with what you would expect ?
>
> best regards
>
> Paolo
>
>
>
>
> On 01/05/2020 23:20, 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
> --
> ---------------------------------------------------------------------
> 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/20200508/b4ec0d46/attachment-0001.html>
More information about the AstroPy
mailing list