Propagating objects in time as part of the transformation framework - #10404

(Now posting in the right thread... Sorry)

@ayshih - this looks very promising indeed! I'll try to bring in @eteq and @adrn again, as I realize now we did in fact have a discussion about some of this when the proper motions were first implemented, and how to handle transformations to a different obstime was explicitly left as a to-do item - see #6280. Not dealing with it has left us in a fine mess...

Your approach of having a context manager may well be the most sensible, and I certainly like in this PR that it seems to properly identify the cause by making explicit the previously implicit ICRS->BCRS transformation.

You are definitely on the right track here. I would suggest that you iterate on the light time delay in SolarSystemLinearPropagator so that applying progagate() followed by depropagate() gets you back to the original position for a fast moving target.

DEFAULT_REF_EPOCH = Time('J2000', scale='tt')

The time scale here should properly be 'tdb'.

c2_V2 = speed_of_light ** 2 - vel.dot(vel)

Again, you may achieve greater accuracy by expressing this as:

c2_V2 = (speed_of_light + np.sqrt(vel.dot(vel))) * (speed_of_light - np.sqrt(vel.dot(vel)))

Shouldn't the SkyCoord itself carry information on whether it is "solar system" or not?

I agree that it might be easier for the user to figure out how to do this if the relevant information could be specified when instantiating an ICRS SkyCoord object. See: #10372 (comment)

I would suggest that you iterate on the light time delay in SolarSystemLinearPropagator so that applying progagate() followed by depropagate() gets you back to the original position for a fast moving target.

I shouldn't need to. The propagate() call calculates the exact solution for linear motion (barring the potential small loss of accuracy you've noted for c2_V2, which I put off fixing), so depropagate() should exactly reverse it. That is, it solves for the time where the location retarded by that time has that same light travel time to the observer. The propagate() code is complex to intentionally avoid an iterated calculation.

If it doesn't exactly reverse, then I've messed up the math.

Thanks @ayshih this does look very promising and necessary! I'll ping @eteq on Slack to see if he can dig in to the details here, as he is more of an expert on this (I understand the general concepts but never have to deal with these things in practice in my own work).

One minor comment: The context manager approach is clever, but feels a little awkward to me. However, I haven't thought of another idea that would work with the current transformation framework, so fine to ignore this comment for now!

Nice work @ayshih !

I very much like the inclusion of an explicit BCRS frame and hope to find time to look at the detailed code. In the meantime though, there's one minor suggestion, and an awkward point raised by your inclusion.

First - minor issue - the docstring for the BCRS class should probably explain that BCRS is identical to ICRS, but objects are corrected for space motion to the specified obstime.

The awkward point is that, according to the sofa docs the transform from ICRS to all other frames should go through BCRS. This would be a breaking API and would also "hide" the space motion correction from users converting a coordinate to/from ICRS. I'm not sure yet what I feel the correct solution is, but thought I should raise the issue...

@StuartLittlefair - our transformations do go through BCRS for any code that relies on SOFA. The problem is that when we included proper motions, a decision was made to not propagate those by default, i.e., they are explicitly left out in the transformation chain (which you pointed out already then was a mistake, but at the time the result was that we raised an issue - #6280 - that we ended up just forgetting about...).

My feeling is still that somehow we have to be able to let the SkyCoord itself be the place which stores how the coordinate is transformed (see #10404 (review) above). If that becomes some kind of option, then the context manager could perhaps determine the default for that option.

@mhvk I don't think they do, and they won't under this PR either. If I define a SkyCoord (frame is ICRS) with proper motion, parallax etc, then the transform to (e.g.) GCRS should be ICRS->BCRS->GCRS.

However, this PR leaves the direct ICRS->GCRS transform in the transform graph (using the function icrs_to_gcrs. That will be the shortest route from ICRS to GCRS and so the following code would bypass BCRS entirely.

# coo is an ICRS SkyCoord with proper motion, parallax etc
coo.transform_to(GCRS(obstime=Time("J2015.5")))

At the moment, and still under this PR as it stands, if you want true GCRS you'd have to do

t = Time("J2015.5")
coo.apply_space_motion(obstime=t).transform_to(GCRS(obstime=t))

@StuartLittlefair - sorry, I was unclear - I didn't mean we could just use our current transformation, but that since they rely on sofa anyway, we can change them such that they take proper motion into account as needed. In the example function you pointed to, it would mean changing the call to atciqz to atciq.

Though obviously how to give a choice is the crux of the problem... Do we need an ICRS_no_pm_propagation frame for getting the current behaviour?!

@mhvk I see what you mean. We could keep the direct ICRS->GCRS step in place, but update all transforms out of ICRS to include space motion.

It seems simpler to have one transform that implements the space motion correction (ICRS->BCRS) and make BCRS the only route in and out of ICRS. This might be slower though.

As for how to give a choice about how this is handled to the user; this is really hard. Remember that frames can also use these transforms, so an implementation should work for frames and skycoords.

I also think it is desirable to allow SkyCoords to have a mixture of solar system objects and stellar objects, and transform them to a single BCRS object correctly. The use case I have in mind is my own robotic telescope scheduler; we have a vectorised SkyCoord object containing all the targets for a night, and we mix and match SS objects with stars!

@StuartLittlefair - nice telescope, it sounds like! For that, though, I think you'd want any proper motions on stars properly propagated (just in case Barnard's star is on the list...). I.e., would it be a problem if all entries were propagated?

No, it would be perfect.

I'm not sure why anyone would not want to apply space motion as the result would be the incorrect position of the star on the date of observation.

I've been thinking a fair bit about this challenge of how to choose which propagation method to apply to ICRS - BCRS transforms. It seems to be the crux of implementing the BCRS frame correctly, and the BCRS frame seems to be the key to nicely handling the implementation of an Astrometric frame as requested in #10230.

I keep coming back to @mhvk 's comment that the SkyCoord itself should hold the information about what type of object is being represented. In a sense we already do that, since SkyCoords have proper motion and radial velocity attributes - i.e we are currently assuming that the SkyCoord represents a star.

I believe the solution to this is to add attributes to SkyCoord that allows us to recognise that it is a Solar System body. My proposal is an attribute solar_system_body_name, or just solar_system_body which follows the requirements of coordinates.get_body - i.e a string that is one of coord.solar_system_ephemeris.bodies, or a list of 2-tuples that specify a path through the JPL kernels.

The presence of stellar motion attributes (proper motion, rvs) would indicate a ICRS <--> BCRS transform using apply_space_motion, whereas the presence of solar_system_body_name would use the solar system machinery.

Support for currently non-supported solar system bodies could be added to coordinates.get_body

(aside:) a consistent implementation of BCRS would also require some changes to apply_space_motion.

SkyCoord.apply_space_motion when applied on ICRS coordinates should return coordinates in the BCRS frame for a start...

Since in a way we're discussing different classes of objects which behave slightly differently, another option would be to introduce a different class to represent them, i.e., define a SolarSystemCoord class or so. I think this would be particularly relevant if we were to link to actual ephemerides rather than just have a near/far flag. A different flag also has the advantage that we can separate concerns a bit: even for stars we do need the option (ideally with a default of True) to propagate proper motion!

I'm keener on retaining SkyCoord for both - Although it's never been clearly defined I've always felt that SkyCoord represents a thing at a particular point and time. It might be a lot of work to introduce a new class, and they would behave identically apart from the method used to progopate coords between obstimes I guess?

For stars I think we should always apply space motion for transforms out of ICRS. Arguably it is a bug that we do not. If we allow this as an option, the default should definitely be true...

Yes, agreed it is a bug not to apply motion...

On the solar-system vs star case, I think the main thing I would like to avoid is to have a decision be made just on distance (though fine if that informs the default!). I slightly prefer the subclass since it is for different objects, and while for stars only one type of behaviour seems to make sense, for solar system objects it seems less clear whether there is only one type (but I realize in part I still do not fully understand what is needed!). Anyway, at some level these are just implementation details.

One more concrete suggestion, though, would be to break up the work in smaller pieces. E.g., first correct the star case (which I think is relatively easy) and only then add the solar-system pieces.

Coming back to this (sorry for the sporadic involvement). I think we need to resolve the sub-class/attribute question as a first step. An attribute would (in principle) allow a vector SkyCoord to contain a mix of object types, so it would be my preference.

Once that is decided, we could split the work into the following steps?

1. Add the attribute to SkyCoord, and update the solar system routines (e.g get_body) to add that attribute automatically
2. Create a private version of apply_space_motion which uses the stellar or planetary routine as appropriate
3. Add the BCRS frame, using the private routine above.
4. Proceed to add an Astrometric frame

The list sounds good. I'm still slightly in favour of a subclass as I think it cleaner and I suspect there may be other behaviour for which some special methods would make a lot of sense for solar-system objects and not so much for sky coordinates. But it is only a mild preference (i.e., as far as I'm concerned, whoever does the work gets to decide!).

For an attribute, it does seem to be asking for trouble if we were to allow a mix in a single SkyCoord. Is there a specific example you were thinking of where that would be useful?

I was thinking of someone trying to do - for example - a lot of Alt/Az calculations with a mix of targets (stellar and planetary). Even for that use case though you could always split your targets into two batches. As you say, allowing that mix sounds problematic.

Thinking of the issues you've seen with Time and TimeDelta I guess much of the current SkyCoord code would move to a CoordBase class, with SkyCoord and SolarSystemCoord inheriting from that?

@mkbrewer - Isn't this PR trying to solve the problem that we don't have a frame to represent solar system objects at an absolute time? As I mentioned, I think "ICRS" for solar system bodies effectively becomes that given the transformation at https://github.com/astropy/astropy/pull/10404/files#diff-86311e5590ca8c380eae0e346fdc47e0736e57601f5e4decc5ab6f5fbd2e7d83R125-R147.

More generally, wouldn't it be useful to have this option as part of the transformation network instead of having the custom _get_apparent_body_position() used by get_body()? If not, then I think we can proceed much more simply by just implementing the ICRS->BCRS transformation for stars! (Actually, I think that would be good regardless - I'm very much in favour of splitting things up in small pieces and doing those that are clearly OK!).

@mhvk There are several aspects to this. First we need to settle on exactly what we mean by a BCRS frame. I would suggest that this represents the ICRS aligned position of the object referred to the SSB, adjusted for space motion (not parallax as that is within the purview of the astrometric frame) and light time to the observer. Second, we need to settle on how to realize this frame. In that regard, it seems to me that there are three classes of objects:

• Extra solar objects defined by having a catalog position at epoch J2000, with or without a distance, components of proper motion and parallax.

• Solar system objects with positions read from an ephemeris at epoch obstime, without any velocity information.

• Solar system objects not read from an ephemeris with both position and velocity information given at a reference epoch within a few days of obstime.

For the third class of objects, I say that the reference epoch must be within a few days of obstime as any solar system object with a velocity less than solar system escape velocity is actually traveling in an elliptical orbit so that using a linear propagator to advance or retard its position using its velocity at the given reference epoch will not be valid for long.

Now the problem becomes simply how to properly identify each class of object in order to realize its position in the BCRS frame at obstime. For this, there are several options. I can think of three (there may be more):

• Assign an attribute specifying the reference epoch similar to the equinox attribute used in the FK5 class. If this is J2000, then the object is a member of the first class. If this is obstime and there is no velocity information, then the object is a member of the second class. If this is within a few days of obstime and there is velocity information then this is a member of the third class.

• Explicitly assign an attribute that denotes the class of the object, such as catalog, ephemeris or non-ephemeris.

• Create actual subclasses for each type of object with individual propagation methods included.

For the first class of objects, the realization would be apply_space_motion(). For the second class of objects, this would be get_body(). For the third class of objects, this would be something similar to the SolarSystemLinearPropagator() defined in this PR. Although, I'll admit that I don't really understand that as it seems to only add light time delay. I don't see it actually propagating the the position from ref_epoch to obstime using the velocity information, but that may be happening in a way that I don't understand.

More generally, wouldn't it be useful to have this option as part of the transformation network instead of having the custom _get_apparent_body_position() used by get_body()?

I would tend to call _get_apparent_body_position() customary rather than custom. I've been dealing with solar system ephemerides for nearly 40 years. From the time when they were distributed on reels of 9 track tape. The algorithm in _get_apparent_body_position() is exactly how it's always been done. I don't think you should change that in any way for the JPL ephemerides. The lookups are quite fast and you want full accuracy. However, you might consider using a linear propagator for the builtin ephemeris to speed up the light time delay calculation so long as doing so doesn't degrade the accuracy too much. Anyone using the builtin ephemeris isn't all that interested in accuracy anyway.

@mkbrewer - thanks! Agreed with needing some way to distinguish at least between objects with ephemerides and those with only velocities, and, between the latter, between close and distant ones (perhaps simply based on how much the velocity over a given amount of time changes the fractional position). I'm tending towards subclasses, but an attribute may end up being better - actually, the distinguishing attribute itself could be apply_space_motion or get_body (well, some wrapper with an identical signature).

On _get_apparent_body_position(): definitely didn't want to change the implementation - I just wondered if it couldn't be made an explicit part of the transformation network. But thinking more, I am getting more instead of less confused, so perhaps the best thing is simply to try to implement a bit of the above instead.

Hopefully, my posts are less confusing than these instructions on how to carve a pumpkin. 😉

Some targeted comments:

The third type of object that we discussed is an object with no catalog position or ephemeris. Just an instantaneous position and velocity. I'm not sure what these objects are. Comets or interplanetary spacecraft come to mind. For this type of object, the proposed linear transformation would likely be fine. A simple two body Keplerian propagator might be better.

When I also suggested "transient phenomena", I had specifically in mind coronal mass ejections (CMEs). The "object" could be the 3D position/velocity of the peak of the erupting magnetic flux rope. A Keplerian propagator wouldn't be physically accurate in this case, although in truth it may not be appreciably worse than a linear propagator.

I think that what you noted, that an ICRS catalog position makes no sense for solar system bodies

Given the uniqueness constraint, astropy's current ICRS frame makes, like SOFA, the implicit assumption of being a catalog position - specifically, the position as would be measured by an imaginary observer sitting at the SSB on J2000. This makes it useless for solar system objects

Now, if I understood correctly, in this PR (see SolarSystemLinearPropagator), for solar system objects, the ICRS frame class is used for such coordinates at an absolute time. But, having (finally) wrapping my head around it, this seems a mistake, and I wonder if we should not instead define a new frame class for coordinates in the BCRS but in which the time always is without light travel time effects.

I disagree with the notion that Astropy's ICRS has an implicit assumption of being a catalog position at J2000. As the transformation framework stands (before this PR, and not counting SkyCoord.apply_space_motion()), obstime is more accurately thought of as the time that defines the frame's origin/axes, which may or may not correspond to the notion of the reference epoch. The fact that ICRS is missing obstime is a reflection of the fact that the frame's definition is unchanging in time. It is valid (outside of SOFA's distinction of ICRS) to talk about the ICRS position of a solar-system body.

I don't know of any place where the current transformation framework has a notion of J2000 with ICRS. This PR adds that notion by having J2000 as the default reference epoch, but only when using the stellar side of the propagators.

First we need to settle on exactly what we mean by a BCRS frame. I would suggest that this represents the ICRS aligned position of the object referred to the SSB, adjusted for space motion (not parallax as that is within the purview of the astrometric frame) and light time to the observer.

Why do you want an Astropy BCRS to account for light time to the observer? And you want the observer to not necessarily be the SSB, so the frame would have an observer frame attribute? Why not defer it to a more observer-aware astrometric frame (e.g., what I tried to do in #10372)?

I don't see it actually propagating the the position from ref_epoch to obstime using the velocity information, but that may be happening in a way that I don't understand.

No, you're correct: ref_epoch is not used for the solar-system side of that propagator. (In the example above, stars would be propagated from the reference epoch of J2000, but that would be a wild choice of reference epoch for Venus.) The assumption is that the solar-system body's location/velocity is true at the specified obstime, and the propagator deals solely with the light-travel-time propagation. More general propagation from one obstime to a different obstime is a grander notion, and definitely would require the Keplerian propagator you've asked for.

p.s. It does seem one of those cases where a quick call would have eliminated a lot of this confusion, mostly in distinguishing the actual ICRS/BCRS/GCRS from the corresponding astropy classes, and how to make explicit assumptions about epoch, observer location, etc.

I support having some sort of meeting to hash this all out, as I think there are still a lot of details to achieve consensus on.

@ayshih Thanks for chiming in on this. I was beginning to think that only @mhvk was interested.

Why do you want an Astropy BCRS to account for light time to the observer? And you want the observer to not necessarily be the SSB, so the frame would have an observer frame attribute? Why not defer it to a more observer-aware astrometric frame (e.g., what I tried to do in #10372)?

Because, that is how it is currently handled. In this PR, you have modified get_body() to return the solar system object's position to be in the BCRS frame and that includes light time delay to the observer. Your SolarSystemLinearPropagator does the same. SOFA's pmpx() also corrects for the light time delay to the observer here:

/* Proper motion time interval (y) including Roemer effect. */
   dt = pmt + iauPdp(p,pob)*AULTY;

I agree that this means that both BCRS and astrometric need to be aware of the observer's position. The alternative is to skip BCRS and go straight to astrometric.

I disagree with the notion that Astropy's ICRS has an implicit assumption of being a catalog position at J2000.

Is this really true? If one defines a SkyCoord in the FK4 frame and then transforms this frame to ICRS, does that transform not result in a position referenced to the epoch J2000?

As the transformation framework stands (before this PR, and not counting SkyCoord.apply_space_motion()), obstime is more accurately thought of as the time that defines the frame's origin/axes, which may or may not correspond to the notion of the reference epoch.

I'm not really understanding this. If the frame's origin and axes are defined by obstime, then the frame is not ICRS. To be ICRS, the origin must be the SSB and the axes must be aligned with those of the ICRS. They can't have a time dependency defined by obstime as that would imply a rotation of the axes such as precession and nutation.

It is valid (outside of SOFA's distinction of ICRS) to talk about the ICRS position of a solar-system body.

Well, we are talking about using SOFA's distinction here to distinguish between ICRS and BCRS. For solar system bodies, one could imagine backing out the light time delay and calling that the ICRS position, but this would only be valid at obstime and would be similar to removing the correction for light time to the observer from space motion for stars, but not the space motion from the catalog epoch to obstime. I don't think this is what you want for a definition of a position realized in the ICRS frame. This is why I suggest either making the BCRS<->ICRS transformation for solar system bodies a no-op or just not defining any BCRS<->ICRS transformation for this class of objects at all. If you have some other proposal for defining an ICRS position for solar system bodies, please feel free to state it.

Actually, backing out the light time delay for a solar system object and calling that the ICRS position is a perfectly reasonable thing to do. This would imply adding a reference epoch to the ICRS frame as I stated here:

Assign an attribute specifying the reference epoch similar to the equinox attribute used in the FK5 class. If this is J2000, then the object is a member of the first class. If this is obstime and there is no velocity information, then the object is a member of the second class. If this is within a few days of obstime and there is velocity information then this is a member of the third class.

Basically, get_body_barycentric() would return the object's position in the ICRS frame, whereas get_body() would return the object's position in the BCRS frame. This could be extended to get_body_barycentric_posvel(), but that would make it difficult to distinguish between ephemeris and non-ephemeris bodies in this paradigm. Some additional information would be necessary to distinguish the two.

@ayshih, @mkbrewer - We may all be fairly aligned. But perhaps it still is helpful to be even clearer about what is needed. In terms of implementation, I think there are actually two items: 1. how positions are defined, and 2. how they are transformed.

1. Position models. Conceptually, I guess for both stars and solar system objects one starts with a model in ICRS coordinates, where the model differs in that for stars it is necessarily in time measured by an observer at the SSB, while for solar-system objects it is generally in absolute time. So, one question is how to distinguish these options in an implementation. The fact that not only the time frame but also the model are generally very different (just pos + vel * (obstime-J2000) vs ephemeris(t_eph)), suggests a different class.

2. The transformations are also different, since for stars one is just transferring between observers, while for planets we need to correct for the travel time (i.e., move from t_eph to obstime). What I'm less clear about is, though, is which transformations this affects. For full precision, do all transformations have to refer to the ephemeris? And are some paths we currently define impossible? To what extent is should solar system objects really have a separate transformation graph? Or would a hook that does the propagation in (some of the) transformations suffice?

@mhvk I'm not sure where you came up with t_eph there. In 1), I would suggest that the model is ephemeris(obstime). In 2), please state what you mean by "transferring between observers". If you mean moving the observer from the SSB to some other location, I fail to see how these transformations are fundamentally different. In the case of stars, this transformation retards or advances the position by (p dot pob) * AULTY where p is the unit vector from the observer to the source, pob is the position vector of the observer wrt the SSB in AU, and AULTY is the light travel time for 1 AU. In the case of planets, this transformation similarly moves the observer from the source to another location retarding the position from obstime to obstime - t_light where t_light is the light travel time from the source to the observer.

@mkbrewer - t_eph is what it is called in the JPL ephemeris (if I recall correctly) - I meant it as short-hand for the TDB time used as input to plan94 or jplephem.

The transformations I'm not as clear about, though it seemed to me that for a given obstime, a transformation from the SSB to any position on earth, for a planet one has to worry about its ephemeris because of the change in light travel time, while for a star one does not worry about the star's proper motion (though thinking about Barnard's star at GAIA-like precision, one perhaps should...).

But really it was mostly about implementation, for which it seems important to recognize that the planetary ephemerides are given in terms of locations at some absolute time, while the stellar ones are given in terms of positions as observed at the SSB. But when you mentioned ephemeris(obstime), perhaps you are assuming that one wraps the ephemeris so that it can give the position at some location and obstime? That is, wrapped by (some version of) get_body()?

p.s. Just a heads up: feature freeze for astropy 4.2 is close, so expect (even) slower responses than usual!

Since this PR has been been mentioned elsewhere twice today, I thought I should affirm that I haven't abandoned this effort, and I'm continuing to mull over ways to proceed. My minor epiphany recently has been that we shouldn't be trying to treat stars and solar-system bodies differently, but rather we should be treating them the same (cf. @mhvk's mention of Barnard's star). The API should instead distinguish between apparent locations and true locations. I've been tinkering on possible API approaches, which may include some breaking changes to SkyCoord.

Well, they are somewhat different as, in the case of stars, proper motion and radial velocity are applied, while for planets, a new ephemeris lookup is required.

@ayshih - sounds good! I think it will also (continue to) be good to try to break things up in small pieces - e.g., first just doing proper motion and velocity, and worrying about full ephemerides later. E.g., just having BCRS and self-transforms of it would seem an excellent first step.

Well, you can't do a BCRS[t1]->BCRS[t2] transform without looking up the position of a planet at the new time. This gets back to the idea of having a new attribute on the SkyCoord that indicates what type of object is associated with it.

@mkbrewer - absolutely! And it is not that I wouldn't want to deal with planets, but rather that I'd like to break things up in pieces, and deal with the simple ones first!

Here is yet another reason to implement an astrometric frame: #11701 (comment)

And another: #11755 (comment)

Hi humans 👋 - this pull request hasn't had any new commits for approximately 1 year, 3 months. I plan to close this in 30 days if the pull request doesn't have any new commits by then.

In lieu of a stalled pull request, please consider closing this and open an issue instead if a reminder is needed to revisit in the future. Maintainers may also choose to add keep-open label to keep this PR open but it is discouraged unless absolutely necessary.

If this PR still needs to be reviewed, as an author, you can rebase it to reset the clock.

If you believe I commented on this pull request incorrectly, please report this here.

@ayshih @mhvk Is there a plan for going forward with this? Are we going to break it into smaller pieces as Marten has suggested?

Hi humans 👋 - this pull request hasn't had any new commits for approximately 1 year, 3 months. I plan to close this in 30 days if the pull request doesn't have any new commits by then.

In lieu of a stalled pull request, please consider closing this and open an issue instead if a reminder is needed to revisit in the future. Maintainers may also choose to add keep-open label to keep this PR open but it is discouraged unless absolutely necessary.

If this PR still needs to be reviewed, as an author, you can rebase it to reset the clock.

If you believe I commented on this pull request incorrectly, please report this here.

Ops, sorry for the double bot post. Hopefully I fixed the bug upstream.

I'm going to close this pull request as per my previous message. If you think what is being added/fixed here is still important, please remember to open an issue to keep track of it. Thanks!

If this is the first time I am commenting on this issue, or if you believe I closed this issue incorrectly, please report this here.