Requesting two additional built-in frames.

Issue #1: The positions of the bodies are returned in the GCRS and there appears to be no way to easily transform them to topocentric astrometric positions. By this I mean simply the difference in the ICRS position of the body compensated for light time and the ICRS position of the observer. This is rather important if one wishes to place the body on a background map in the ICRS. It is also the only way to compare the output of AstroPy's ephemerides with that of JPL Horizons or Brandon Rhodes' Skyfield. Why isn't there a built-in frame for doing this?

Issue #2: Currently, there is also no built-in frame for transforming the returned positions to classical apparent place. By this I mean the topocentric position with respect to the true equator and equinox of date. This is quite important to people such as myself who still like their origin of right ascension to be an actual location on the sky rather than a convenient mathematical construct. It allows one to point an equatorial mounted telescope using the local sidereal time to calculate the hour angle. And again, this is the only way to compare the output of AstroPy's ephemerides with that of JPL Horizons or Skyfield. It is also quite simple to do. Just adjust the CIRS right ascension by subtracting the equation of the origins. Note: I did find a function for doing this in solar_system.py, _apparent_position_in_true_coordinates(), but it feels sort of kludgy to use this. There should be a built-in frame for this.

Oh, I forgot to mention that this is in regard to the solar system ephemerides.

I looked into the first one and it's exactly the same as icrs_to_gcrs() without the call to atciqz(). Alternately, the transformation from GCRS is exactly the same as gcrs_to_icrs() without the addition of astrom_eb.

+1

The frame suggested in #2 would also be extremely helpful for radio astronomers who are used to thinking about source locations in hour angle.

Radio astronomer since 1975 here.

@StuartLittlefair

Why isn't there a built-in frame for doing this?

Probably simply because those who contributed to astropy have not yet needed it. Most of us are volunteers and few work on solar-system bodies (e.g., I'm not sure all transformations are safe for nearby objects).

But we very much appreciate contributions. Indeed, in #10149 we have just added a true equator, mean equinox frame following a similar request - hence, cc @StuartLittlefair.

From what you write, it sounds like the frames that might be useful for you may not be too difficult to add - though I'll readily admit that I'm not completely sure what exactly you mean. (So far, much of what we have follows the IAU/SOFA documentation; it might help to point to a formal definition.)

Marten, I understand that these concepts may be somewhat unclear to those who aren't intimately familiar with them. Astrometric coordinates are simply ICRS coordinates translated from the SSB to the location of the observer. Hence, for stars they include proper motion and parallax and for solar system bodies they include the light time to the observer's location. They are exactly the same as GCRS coordinates without the addition of aberration and gravitational light deflection. I would refer you to Figure 1 in the SOFA Astrometry Tools document. As I previously stated, the transformation from astropy's ICRS frame to astrometic is exactly the same as the icrs_to_gcrs() transformation in coordinates.builtin_frames.icrs_cirs_transforms.py except that the call to atciqz() is omitted. One other thing that I noticed when looking into this is that the light time calculation in the solar_system.py function _get_apparent_body_position() is slightly off as it calculates the light time from the target body to the Earth geocenter rather than the observer's location, but this is a very minor issue for an observer on Earth.

…

On 05/01/2020 08:32 PM, Marten van Kerkwijk wrote: Why isn't there a built-in frame for doing this? Probably simply because those who contributed to astropy have not yet needed it. Most of us are volunteers and few work on solar-system bodies (e.g., I'm not sure all transformations are safe for nearby objects). But we very much appreciate contributions. Indeed, in #10149 <#10149> we have just added a true equator, mean equinox frame following a similar request - hence, cc @StuartLittlefair <https://github.com/StuartLittlefair>. From what you write, it sounds like the frames that might be useful for you may not be too difficult to add - though I'll readily admit that I'm not completely sure what exactly you mean. (So far, much of what we have follows the IAU/SOFA documentation; it might help to point to a formal definition.) — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#10230 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACRERNBFIT3ANG4GAYK5YYDRPNS2NANCNFSM4MT4P6VQ>.

I'll chime in from the SunPy side that we'd also be interested in good implementation of this. For example, precise observations of Venus transiting the Sun require proper handling of the light travel time to an observer not necessarily sited at Earth. (We currently have our own function, get_body_heliographic_stonyhurst(), which can calculate the light-travel-time correction for an arbitrary observer, not just an Earth observer.)

That said, I'd worry that an actual astrometric frame might be a bit awkward to work with transformation-wise. It's straightforward enough to calculate the light-travel-time correction for known objects (because you can look up their location at the earlier time), but an arbitrary coordinate would need to have a velocity attached to perform the (at minimum) linear propagation of the coordinate in time. Observationally, sometimes you don't even know the distance to the object to calculate the light travel time, much less its velocity.

@mkbrewer

@mkbrewer, @ayshih - OK, yes, that really helps!

@mkbrewer - I opened #10271 for the inclusion of the location of the observer on Earth. As you note, it is a very small effect, but I agree it is missing and we should document that at the very least.

I do share @ayshih's worry about how to actually do the travel-time correction in practice. But I think the simpler part of the request is to implement the "astrometric" frame that sits in between BCRS and GCRS in Fig. 1 of SOFA's astrometric tools; with that, one could at least more easily do the calculations oneself. Is that correct? And would its standard name indeed be "Astrometric"?

p.s. The second frame also seems eminently do-able; it seems our use of the internal _apparent_position_in_true_coordinates() routine is in fact simply to facilitate comparison with skyfield... I think the main work may be in documenting when the various frames of reference are useful (I find the SOFA guide hard-going). And perhaps refer to skyfield for purposes where that is the better tool (we already rely on jplephem for the JPL ephemerides).

Light travel time is not an issue here. That is handled by get_body() for a solar system target body and an Earth based observer. For an off Earth observer, something similar would be used. I'm not sure how you are handling parallax and proper motion in these transforms as you are calling atciqz() rather than atciq(), but however you handle that for ICRS -> GCRS, it would be identical here.

"Astrometric" is appropriate. SOFA, Skyfield and Horizons all use this term. And, after all, generations of astronomers are/were familiar with Astrometric and Apparent positions in the Astronomical Almanac (Yeah, I know that's pretty obsolete at this point. Just pointing out that it's a very traditional term.)

Calling the second frame I'm asking for "Apparent", on the other hand, could lead to push back as the CIRS is supposed to be the "new" apparent, whereas true equator true equinox is the "old" apparent. Many observatories still use it though, so it's likely not going away any time soon. Perhaps taking a cue from TEME and calling it TETE would be a possibility (although that does sounds a bit French ;-} ). .

@eteq

Great, this is beginning to sound like a plan! I'm happy to go with Astrometric and TETE is probably OK too - maybe easiest done in 2 separate stages, since Astrometric sounds quite a bit easier. Let me bring in @eteq here, since he did the original implementation of the whole set of transfers in the SOFA writeup.

I do worry about your question about how we handle proper motion and RV in transformations - the machinery predates having those on coordinates, so it may have been a short-cut that we should have remedied. Again, @eteq and perhaps @adrn and/or @StuartLittlefair can help here.

TETE is quite simple also. All you need is the equation of the origins which is extracted from the nutation x precession x bias matrix calculated by pnm06a() by the function eors(). The way that I would handle this is to modify your function get_cip() to return the equation of the origins (eo) as well. Then IERS->TETE is exactly the same as your transformation iers_to_cirs() except that you subtract eo from cirs_ra to yield tete_ra.

Great! If you are up to it, a PR would be very welcome... (but do not feel obliged, I think it is now clear we should do it and we will get to it eventually).

For what it's worth, TETE lines up with what the frame is called in the USNO circular 179: https://www.usno.navy.mil/USNO/astronomical-applications/publications/Circular_179.pdf

The (unnumbered) figure on page 57 of that document has been helpful for me.

Thanks! looks like we've got a nice extra reference for in the docstring as well...

xref #10292 (comment)

It is unfortunate that you only implement geocentric position for your CIRS frame. You use apci() and apio() rather than the more rigorous apco() which builds the observer's topocentric position and velocity (if present) right in. This has two side effects. In going from ICRS to CIRS light deflection doesn't take into account a topocentric observer's position and velocity and in going from CIRS to AltAz the less rigorous simple first-order vector expression for diurnal aberration as mentioned in paragraph 4.7 of the SOFA Astrometry Tools document is used.

On the other hand, your GCRS frame does implement topocentric position rigorously. Therefore, I find it curious that you don't derive your CIRS frame from GCRS. After all, with GCRS you are just one step of applying the CIO based BPN matrix from true topocentric CIRS intermediate (apparent) place. The call to apcs() in this transform properly sets up the ASTROM data structure for rigorously calculating gravitational light deflection, diurnal aberration and diurnal parallax including the observer's position and velocity (including taking the observer's velocity into account in beta - 1). I particularly like the way that diurnal aberration is naturally included in Ab() this way, so that the transformation from topocentric CIRS to AltAz can be done without resorting to the use of astrom->diurab.

Now you may think that only a cranky old celestial mechanic like myself would complain about such small errors, but SOFA is supposed to be definitive and capable of sub-mas accuracy, so astropy really should reflect that.

Since TETE should definitely offer both geocentric and topocentric position, given the situation with your implementation of CIRS, it would be best to derive TETE from GCRS using the equinox based BPN provided by Pnm06a() to do it correctly.

@mkbrewer

@mkbrewer - we've only done as best as we can - I remember @eteq, who did most of the work, very much struggling with it... And I think I can speak for all of us saying that we really appreciate being told how to do it even better! I think we might as well get CIRS right as well, so opened a new issue, #10356

That said, I'd worry that an actual astrometric frame might be a bit awkward to work with transformation-wise. It's straightforward enough to calculate the light-travel-time correction for known objects (because you can look up their location at the earlier time), but an arbitrary coordinate would need to have a velocity attached to perform the (at minimum) linear propagation of the coordinate in time. Observationally, sometimes you don't even know the distance to the object to calculate the light travel time, much less its velocity.

Addressing this comment, there are basically two cases.

For stars we have the published values of proper motion, radial velocity and annual parallax at epoch. The astrometric position of the star is calculated from these by Pmpx() called by Atciq() (but not Atciqz()).

For solar system bodies, things are a bit more complicated. Any object moving through the solar system is following an elliptical, hyperbolic or parabolic trajectory. It is very simple to calculate this trajectory from the orbital elements provided by the Minor Planet Center by solving Kepler's equation. In fact, I have a program to do this that consists of only a few hundred lines of C code.Since much of this code is involved in calculating the eccentric anomaly which is a simple minimization problem, this code could be condensed in Python by using SciPy's minimization algorithms. Since small bodies' orbits are constantly being perturbed, however, this simple solution should be treated like the solutions for satellites given TLEs (i.e. only accurate for a few weeks or even days). That said, we always found this solution perfectly adequate for comet and asteroid observations so long as the orbital elements were updated frequently enough. For longer term ephemerides, the perturbations must be taken into account using numerical integration. One such code is offered by Bill Gray in his lunar repository, but that is likely well beyond the scope of astropy.

For solar system bodies, things are a bit more complicated. Any object moving through the solar system is following an elliptical, hyperbolic or parabolic trajectory. ...

My concern is actually about how any trajectory-calculation machinery would be integrated into the coordination-transformation framework. As it currently stands, an Astropy coordinate is best thought of as attached to a location in inertial space rather than a (moving) body. (Aside from when using SkyCoord.apply_space_motion(), but that's not really integrated into the framework.) Adding on the functionality of tracking a solar-system body raises questions such as:

What new sorts of frame attributes would be needed?
- The name of the body, as well as where to retrieve the orbital elements given that name?
- The full orbital elements for a body, if the body can't be looked up?
What if the user can't provide orbital elements, and only has the instantaneous velocity?
- Do we propagate it as if the body was moving linearly, or do we make a creative assumption about the nature of its orbit?
- What if the user knows the 3D location but only the projected 2D velocity?
- What if the user knows location and velocity only in 2D?

Of course, the transformation framework can raise errors in some of these situations, but what it can and should support needs to be carefully considered.

@astrojuanlu

My sense would be that for astropy we'd stick with trying to do the coordinates as well as possible with just a proper motion for time-of-flight effects - beyond that it would become too specialized and more suited for an affiliated package focussed on the solar system. Indeed, maybe poliastro? Which reminds me to ping @astrojuanlu - who may be interested in this!

As it currently stands, an Astropy coordinate is best thought of as attached to a location in inertial space rather than a (moving) body

This isn't true for planets, so any solar system body would be treated the same way. To do this properly, you'd need to provide a way for the user to define a solar system body (say a comet or asteroid) by providing a name and a set of orbital elements in order to calculate its position and velocity at any given time. Once you have those, the light time delay can be calculated exactly as is done for the planets.

My sense would be that for astropy we'd stick with trying to do the coordinates as well as possible with just a proper motion for time-of-flight effects.

Proper motion doesn't actually have anything to do with time-of-flight. The catalog coordinates of a star are as measured from Earth. They already include time-of-flight. Proper motion is just derived from the change in measured position over time.

None of this is to suggest that you should implement this. I'm just describing how it could be done is all.

@mhvk

Thanks for the ping @mhvk! Trying to give an organized stream of thoughts from what I read:

The questions that @ayshih raises in his last comment are very relevant. At the moment, the way we deal with this in poliastro is by having:

Functions that convert between cartesian elements and Keplerian classical elements
Several propagators that take these elements and, under some customizable assumptions (i.e. pure two body motion, selected perturbations) return the position of a body in a future instant

And the different cases are:

For Solar System bodies, we query for the ephemerides (position + velocity) at the desired instant directly using astropy + jplephem, therefore using only (1)
For major non planetary bodies we do the same but using astroquery
For some minor bodies, we use sbpy and then propagate the given orbit to the desired state, therefore using (1) and (2)

In my opinion baking all this complexity in Astropy to make a magical solution that just works it's quite impractical (which in any case is not what @mkbrewer is suggesting anyway, as per his last comment)

So, perhaps a good solution from an API point of view would be to allow the user to pass optional (arrays of) position and velocity, and use them if given to compute the necessary astrometric corrections? That way, the user can decide how to compute such extra parameters - for example using the methods I outlined above.

@ayshih

Although it is really quite simple to solve Kepler's equation for a two body (the target body and the Sun) trajectory given a set of orbital elements, I'd like to re-emphasize that I was just responding to the original comment by @ayshih and this one by @StuartLittlefair in issue #9362.

Thinking about this issue, I'm not sure we can do the "right thing" for fast moving solar system objects.

What is missing in our transformations is the planetary aberration term, or correcting for the fact that the light that reaches us at obstime was actually emitted earlier, when the body had a different position.

For solar system objects we have ephemerides for we correct for this since we can predict their location at an arbitrary time. For other fast moving objects, even if the user provides instantaneous velocity information we can't correct unless we assume motion in a straight line.

@mkbrewer

Although it is really quite simple to solve Kepler's equation for a two body (the target body and the Sun) trajectory given a set of orbital elements

I beg to disagree @mkbrewer - sure, it's quite simple to solve Kepler's equation, but to do it in a performant and robust way in all eccentricity regimes is still a research topic on its own. In particular, the near parabolic region (eccentricity around 1.0) is numerically unstable and it's not at all trivial how to deal with it. Many asteroids and comets have this property so it's not a weird corner case. I suffered this in poliastro very recently.

But anyway, I am not sure if this is still on topic :)

That's interesting. I guess I've never had to deal with eccentricities that close to 1.00. I am curious as to how far off the solution would be if one just used a simple parabolic orbit in that case to propagate the trajectory for a few days until new orbital elements became available. And yes, we are getting very OT here.

@mkbrewer

@mkbrewer I am currently wrangling with implementing a TETE frame, which has been relatively simple to do but I would like to test my implementation. Could you provide your code output for the position of the moon over a day? I'm thinking of the plots you shared in #10666. I have some files you shared but cannot replicate those plots...

Navigation Menu

Repository navigation

Requesting two additional built-in frames. #10230

1 remaining item

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Requesting two additional built-in frames. #10230

Description

Activity

mkbrewer commented on Apr 29, 2020

mkbrewer commented on Apr 29, 2020

bhazelton commented on May 1, 2020

mkbrewer commented on May 1, 2020

mhvk commented on May 2, 2020

mkbrewer commented on May 2, 2020

ayshih commented on May 2, 2020

ayshih commented on May 2, 2020

mhvk commented on May 2, 2020

mkbrewer commented on May 2, 2020

mhvk commented on May 2, 2020

mkbrewer commented on May 2, 2020

mhvk commented on May 2, 2020

1 remaining item

bhazelton commented on May 3, 2020

mhvk commented on May 3, 2020

pllim commented on May 11, 2020

mkbrewer commented on May 15, 2020

mhvk commented on May 15, 2020

mkbrewer commented on May 17, 2020

ayshih commented on May 18, 2020

mhvk commented on May 18, 2020

mkbrewer commented on May 18, 2020

astrojuanlu commented on May 18, 2020

mkbrewer commented on May 18, 2020

astrojuanlu commented on May 18, 2020

mkbrewer commented on May 18, 2020

StuartLittlefair commented on Oct 15, 2020

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions