Skip to content

Astrometric versions of coordinate frames - #10372

Closed
ayshih wants to merge 2 commits into
astropy:mainfrom
ayshih:astrometric
Closed

ayshih wants to merge 2 commits into
astropy:mainfrom
ayshih:astrometric

Conversation

@ayshih

@ayshih ayshih commented May 19, 2020 •

Copy link
Copy Markdown
Contributor

This PR was inspired by discussions in #10230. Ping @mhvk @mkbrewer @astrojuanlu

This PR introduces AstrometricFrame for creating an astrometric version of a coordinate frame, relative to a specified observer location. (It does so in a similar manner to SkyOffsetFrame.) These astrometric frames integrate astrometric calculations directly into the transformation framework.

  • For solar-system bodies, the coordinate is propagated backwards in time assuming linear motion to account for light travel time (i.e., "planetary aberration"). The documentation includes a table showing the size of errors due to the assumption of linear motion.
  • For cosmic objects, proper motion is added for the difference between the observation time and the specified reference epoch (basically a wrapper to SkyCoord.apply_space_motion()).

There's a tutorial in the PR, and here are some excerpts:

Venus

    >>> earth = SkyCoord(0*u.deg, 0*u.deg, 0*u.AU,
    ...                  frame='gcrs', obstime='2012-06-06 04:07:29')
    >>> ast_gcrs = AstrometricFrame(observer=earth)
    >>> ast_gcrs
    <AstrometricGCRS Frame (observer=<GCRS Coordinate (obstime=2012-06-06 04:07:29.000, obsgeoloc=(0., 0., 0.) m, obsgeovel=(0., 0., 0.) m / s): (ra, dec, distance) in (deg, deg, AU)
        (0., 0., 0.)>, ref_epoch=J2000.000)>

    >>> venus  # instantaneous position/velocity, not corrected for any light travel time
    <SkyCoord (GCRS: obstime=2012-06-06 04:07:29.000, obsgeoloc=(0., 0., 0.) m, obsgeovel=(0., 0., 0.) m / s): (ra, dec, distance) in (deg, deg, km)
        (74.22533319, 22.7724864, 43190032.91493926)
     (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)
        (-7.87423767e+08, -4.03196311e+08, 0.09882208)>

    >>> venus_ast = venus.transform_to(ast_gcrs)
    >>> venus_ast
    <SkyCoord (AstrometricGCRS: observer=<GCRS Coordinate (obstime=2012-06-06 04:07:29.000, obsgeoloc=(0., 0., 0.) m, obsgeovel=(0., 0., 0.) m / s): (ra, dec, distance) in (deg, deg, AU)
        (0., 0., 0.)>, ref_epoch=J2000.000): (ra, dec, distance) in (deg, deg, km)
        (74.23245792, 22.77361214, 43190037.9814361)
     (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)
        (-7.87503393e+08, -4.03424187e+08, 0.09816952)>

A star

    >>> star_ref = ICRS(ra=10*u.deg, dec=45*u.deg, distance=100*u.pc,
    ...                 pm_ra_cosdec=10*u.arcsec/u.yr, pm_dec=-20*u.arcsec/u.yr,
    ...                 radial_velocity=100*u.km/u.s)
    >>> star_ref
    <ICRS Coordinate: (ra, dec, distance) in (deg, deg, pc)
        (10., 45., 100.)
     (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)
        (10000., -20000., 100.)>

    >>> ast_icrs = AstrometricFrame(observer=ICRS(0*u.deg, 0*u.deg, 0*u.AU),
    ...                             obstime=earth.obstime, ref_epoch='J2010')
    >>> ast_icrs
    <AstrometricICRS Frame (observer=<ICRS Coordinate: (ra, dec, distance) in (deg, deg, AU)
        (0., 0., 0.)>, obstime=2012-06-06 04:07:29.000, ref_epoch=J2010.000)>

    >>> star_ref.transform_to(ast_icrs)
    <AstrometricICRS Coordinate (observer=<ICRS Coordinate: (ra, dec, distance) in (deg, deg, AU)
        (0., 0., 0.)>, obstime=2012-06-06 04:07:29.000, ref_epoch=J2010.000): (ra, dec, distance) in (deg, deg, pc)
        (10.00953953, 44.98650549, 99.99978632)
     (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)
        (9997.59524072, -20001.07628946, 102.79019574)>

This PR needs tests and even more documentation, but it's more than ready for feedback.

@ayshih

ayshih commented May 19, 2020

Copy link
Copy Markdown
Contributor Author

The name is certainly up for discussion too, particularly since the name AstrometricFrame was used in past versions of Astropy.

@mkbrewer

Copy link
Copy Markdown
Contributor

The proper way to do this from ICRS is exactly like icrs_to_gcrs() except replacing the call to atciqz() with a call to Pmpx(). For coordinates with no distance, this will take care of space motion and annual parallax. For a call that includes a distance, you can zero out px in the call to Pmpx() as you are accounting for annual parallax by subtracting the observer's position astrom_eb. This call to Pmpx() should actually be included in icrs_to_gcrs() also and would be if you used atciq() in that function.

For solar system bodies, you already have these in the GCRS frame, so the fastest way to get the astrometric position is to use gcrs_to_icrs() without adding on astrom_eb at the end.

As a side note, the expression (x ** 2 - y**2) should always be evaluated as (x - y) * (x + y) to avoid round off error. I noticed that in your expression for c2_V2.

@mkbrewer

Copy link
Copy Markdown
Contributor

Oh, never mind about zeroing px. I see that is used to incorporate radial velocity into space motion also. In order to only disable annual parallax in the call to Pmpx() you actually need to zero out pob.

@mkbrewer

mkbrewer commented May 20, 2020 •

Copy link
Copy Markdown
Contributor

I looked into this more carefully and I see that the function apply_space_motion() does what Pmpx() does. You just derive annual parallax from distance is all, so you are good for stars here. Sorry for the confusion on my part about that.

On the other hand, for solar system bodies you are assuming that the given coordinates do not already include light time delay. This is incorrect for the very important case of planets, the Sun and the Moon. The GCRS positions of these returned by get_body() already include light time delay. Converting these positions to the ICRS frame removes aberration and gravitational light deflection from the Sun, then translates the origin to the SSB, but doesn't remove light time delay. The only thing that you need to do in this case is translate the origin back to the observer.

For solar system objects other than these, the situation becomes unclear. If given an ICRS frame containing position, distance within the solar system and velocity, you don't a priori know whether or not these already include light time delay to the observer.

@ayshih

ayshih commented May 20, 2020

Copy link
Copy Markdown
Contributor Author

On the other hand, for solar system bodies you are assuming that the given coordinates do not already include light time delay.

Indeed, I should probably add a giant warning. As you note, people should not combine both get_body() and AstrometricFrame, because that would over-correct for light travel time. Note that in my tutorial, I instead create the coordinate for Venus using get_body_barycentric_posvel(), which does not factor in any light-travel-time delay.

@ayshih

ayshih commented May 20, 2020

Copy link
Copy Markdown
Contributor Author

I'll point out that the analogous function to get_body() in SunPy – get_body_heliographic_stonyhurst() – takes what I feel is a better approach. It returns the instantaneous location by default, and only if an observer coordinate is specified does it adjust for light travel time to that observer.

Of course, from the coordinate itself, it's impossible to divine whether light-travel-time delay has been included or not, as you also note. The user presumably has to bear some responsibility here.

@mhvk

mhvk commented May 20, 2020

Copy link
Copy Markdown
Contributor

For ICRS, is there a standard on whether light time delay has been included? By analogy with stars, I would think it would be the position as observed at the SSB at a given time, but then I'm known to get things like these wrong...

@ayshih

ayshih commented May 20, 2020 •

Copy link
Copy Markdown
Contributor Author

For ICRS, is there a standard on whether light time delay has been included? By analogy with stars, I would think it would be the position as observed at the SSB at a given time, but then I'm known to get things like these wrong...

I hope the "standard", if there is one, is to have no light-travel-time delay included. The SSB is not a sane choice for the observer location (for solar-system bodies).

For what it's worth, JPL Horizons returns the instantaneous positions of bodies when asking for vectors, so it's straightforward to avoid having light-travel-time delay baked in.

@ayshih

ayshih commented May 20, 2020

Copy link
Copy Markdown
Contributor Author

Incidentally, it occurred to me to check how inaccurate the linear-motion approximation would be for Parker Solar Probe during a perihelion and an Earth observer. It looks like it's off by ~36 km and ~43 mas. For most purposes, that level of discrepancy may be acceptable.

For an Earth observer, I suspect the worst case for something orbiting the Sun may be Sun-grazing comets. The largest discrepancies would presumably be for satellites orbiting other planets.

@mhvk

mhvk commented May 20, 2020

Copy link
Copy Markdown
Contributor

@ayshih - I think it is down to what you noted before, that we are missing infrastructure. Right now, everything is based around obstime, so I would argue light travel delay is included implicitly (and arguably that is the correct way of thinking about it - there is no absolute time after all).

The solution for having an "instantaneous" position would seem to have a way to indicate that the observer is at the object itself (possibly just having obsgeoloc and obsgeovel to be of the planet itself? Though that would make the "geo" in those a bit of a misnomer).

Anyway, far from being able to properly wrap my head around this, but the need for "giant warnings" is worrying.

@mhvk

mhvk commented May 20, 2020

Copy link
Copy Markdown
Contributor

p.s. I'll try to have a more proper look instead of reacting to the discussion!

@mkbrewer

mkbrewer commented May 20, 2020 •

Copy link
Copy Markdown
Contributor

One further question regarding astrometric star positions. I see that apply_space_motion() doesn't actually apply parallax. Is this done in the final transformation to the reference coordinates? For instance if one specifies an EarthLocation as the astrometric frame origin, will the final transformation subtract the observer's location wrt the SSB from the star position thereby handling parallax implicitly? Not being an astropy expert, this isn't obvious to me from looking at the code.

Indeed, I should probably add a giant warning. As you note, people should not combine both get_body() and AstrometricFrame.

Well, one major point of my request for an astrometric frame was to be able to compare the positions returned by get_body() with those returned by JPL Horizons and Skyfield. There is currently considerable confusion in that regard. People call get_body() and just assume that the returned positions are, in fact, astrometric. I've seen this happen many times. If they even bother to check against Horizons, they wonder why the positions don't match. In the long run, it may be best to just have get_body() return planet positions in the astrometric frame rather than GCRS. After all, the main purpose of get_body() is to interact with a planetary ephemeris to calculate an accurate light time delayed position wrt the given location which could then be the origin of the returned astrometric frame. Since this astrometric position is only one step of applying aberration and gravitational light deflection from the GCRS frame, it would be quite natural to introduce astrometric to GCRS as a separate transformation (see Figure 1 in the SOFA Astrometry Tools documentation). In fact, this might be a way to solve the problem with solar system bodies for which no ephemeris is available also. You could create a function apply_light_time() similar to apply_space_motion() that takes an origin as an argument and returns an astrometric frame with the position adjusted linearly with velocity for objects less than 1 light year distant that contain both position and velocity information. This wouldn't apply to planets as they would already be in the astrometric frame and don't contain any velocity information in any case.

By the way, considering the rather huge error in position and distance that you report for Venus, I took your tutorial as a lesson on what not to do. Astropy really needs to aim for less than 1 mas accuracy for planets.

For what it's worth, JPL Horizons returns the instantaneous positions of bodies when asking for vectors, so it's straightforward to avoid having light-travel-time delay baked in.

Are you sure of this? Horizons vector output includes a light time delay and I've always assumed that this was applied to the position and velocity vectors.

@mkbrewer

Copy link
Copy Markdown
Contributor

Oh, I do see that there are three choices in Horizons under Optional vector-table settings:

Geometric states (no aberration; instantaneous ephemeris states)
Astrometric states (includes target to obs. light time aberration)
Apparent states (astrometric states with stellar aberration)

@mkbrewer

Copy link
Copy Markdown
Contributor

get_body() is currently earthcentric, but I could envision that constraint being relaxed at some point. If the origin was required to be wrt the SSB, then get_body() could return astrometric positions of the planets wrt to an observer in space, on the Moon, Mars, etc.

@mkbrewer

Copy link
Copy Markdown
Contributor

Actually, have you tried this with the output of get_body()? I should think that it would work fine as currently implemented. Since there is no velocity information, there is no danger of adding an additional light time delay. If it works, the result should compare directly to the output from Horizons.

@ayshih

ayshih commented May 20, 2020

Copy link
Copy Markdown
Contributor Author

Actually, have you tried this with the output of get_body()? I should think that it would work fine as currently implemented. Since there is no velocity information, there is no danger of adding an additional light time delay.

True. Because get_body() doesn't return the velocity – although that functionality would be straightforward to add – there wouldn't be an over correction for light travel time. However, using any coordinate with no velocity information with AstrometricFrame would not be particularly fruitful because, for example, the GCRS<->AstrometricGCRS transformations would do nothing.

Responses to earlier comments, working backwards:

Oh, I do see that there are three choices in Horizons under Optional vector-table settings

I'd forgotten that options other than "geometric" were even available...

Astropy really needs to aim for less than 1 mas accuracy for planets.

AstrometricFrame shouldn't ever be the only way for a user to get the astrometric location of a planet. I'd pitch AstrometricFrame as allowing a user to have motion propagation integrated directly into the transformation framework, with the price being a small (sometimes negligible) loss of accuracy. For example, when it comes to the positions of planets in solar imagery, we can accept far worse than 1 mas, and AstrometricFrame would be fine. Of course, we ought to document the level of inaccuracies and point users to how to get more accurate answers if they need them.

I see that apply_space_motion() doesn't actually apply parallax. Is this done in the final transformation to the reference coordinates? For instance if one specifies an EarthLocation as the astrometric frame origin, will the final transformation subtract the observer's location wrt the SSB from the star position thereby handling parallax implicitly?

For AstrometricFrame, maybe the keyword should be renamed from origin to observer. If the observer location is also the origin of the coordinate frame (e.g., GCRS with an appropriately set obsgeoloc), then the parallax will be handled implicitly as you say. Things may get a little non-intuitive if the observer location is not at the origin of the base coordinate frame.

Right now, everything is based around obstime, so I would argue light travel delay is included implicitly (and arguably that is the correct way of thinking about it - there is no absolute time after all).

I think this perspective would majorly conflict with the current design of the transformation framework. For the most part, despite its name, I think of obstime as the coordinate-frame time, not as an observation time. When you convert from a frame with one obstime to a different frame with a different obstime, you're referring to the same point in inertial space using differently defined coordinate axes. You're not opening the can of worms of light travel time.

For coordinates with no distance, this will take care of space motion and annual parallax.

I'll point out that AstrometricFrame as currently implemented does not accept coordinates with no distance (i.e., 2D coordinates in the sky). I use the distance to determine whether the coordinate is in the solar system or not. I opted not to assume that any 2D coordinate referred to a cosmic object.

@mkbrewer

mkbrewer commented May 20, 2020 •

Copy link
Copy Markdown
Contributor

for example, the GCRS<->AstrometricGCRS transformations would do nothing.

Why do you say that? The only velocity used in icrs_to_gcrs() is that of the observer wrt the SSB. This is used in beta - 1 for aberration and gravitational light deflection. All of this is handled by the call to apcs() where the observer's position and velocity wrt the geocenter (if any) are taken from the GCRS frame. No target body velocity information is needed. An astrometric to GCRS transformation could do the same except that the observer's position and velocity could be derived from the astrometric frame's origin directly.

Things may get a little non-intuitive if the observer location is not at the origin of the base coordinate frame.

I'm afraid that I don't understand this. Aren't there just two choices, topocentric or geocentric observer? I suppose other choices would be an orbiting observatory such as Hubble which is really just an extension of topocentric or an observatory in orbit around the Sun which would also experience annual parallax, but that could easily be dealt with as well.

I'll point out that AstrometricFrame as currently implemented does not accept coordinates with no distance.

I understand that, so that statement is moot assuming that any catalog star with parallax also has a distance included. The reason for including parallax in star catalogs is so that annual parallax can be computed without distance information as is done by Pmpx(), but astropy doesn't use that method, so distance is required in order to compute annual parallax.

Regarding accuracy, I am entirely focused on getting planetary astrometric positions at full SOFA accuracy, whereas your focus seems to more on comets, asteroids and spacecraft, so we may need to agree to disagree on that.

@mhvk

mhvk commented May 20, 2020

Copy link
Copy Markdown
Contributor

Regarding accuracy, I am entirely focused on getting planetary astrometric positions at full SOFA accuracy.

The goal in astropy is to do exactly that!

But now having looked a bit at the PR, I got a bit confused. In #10230, the Astrometric frame that was discussed was, at least as I understood it (#10230 (comment)), a frame that sits between BCRS and GCRS in the chain of transformations, and really is just was something that would only need calls to SOFA routines to calculate.

Now I think that I understand that in principle, "astrometric" can make sense for other systems, but just so I get things right: in the context of this PR, the one in Fig. 2 of the SOFA tools documentation would be AstrometricGCRS, is that correct? [I'll add that I find the SOFA documentation a bit confusing in its nomenclature: simply "astrometric" in figure 1, but in the example in section 6 (page 15), it is called "astrometric place" and in the code on page 17 it is called "ICRS (astrometric)"]

In any case, are there fewer issues with that specific frame? E.g., can one do without the assumptions about distance for solar system vs stars? (Which I'd really rather not have anywhere!) If so, would it make sense to just only define that one for now, since it can be constructed directly from SOFA routines?

p.s. Happy to be told that I'm misunderstanding things, but if that happens to a poor astronomer like me who actually has written papers on astrometry (neutron star parallaxes with HST), it is clear that this better all be documented excruciatingly well!

@ayshih

ayshih commented May 21, 2020 •

Copy link
Copy Markdown
Contributor Author

(I've renamed the primary keyword argument to observer to avoid the incorrect impression that the specified observer location becomes the origin of the coordinate frame, at least as currently implemented.)

However, using any coordinate with no velocity information with AstrometricFrame would not be particularly fruitful because, for example, the GCRS<->AstrometricGCRS transformations would do nothing.

...
An astrometric to GCRS transformation could do the same except that the observer's position and velocity could be derived from the astrometric frame's origin directly.

I was referring specifically to the implementation in this PR, where AstrometricGCRS is "simply" adding proper motion and light-travel-time delay effects to GCRS coordinates. Thus, if the coordinate has no velocity, the transformation does nothing. To be clear, AstrometricGCRS is not the same as SOFA's "astrometric" frame (see below).

Things may get a little non-intuitive if the observer location is not at the origin of the base coordinate frame.

I'm afraid that I don't understand this. Aren't there just two choices, topocentric or geocentric observer?

This PR currently chooses not to shift the origin of the astrometric version of the base coordinate frame to the observer location. It won't matter for GCRS if it is defined with the correct obsgeoloc, but it will matter for ICRS (see below).

Regarding accuracy, I am entirely focused on getting planetary astrometric positions at full SOFA accuracy, whereas your focus seems to more on comets, asteroids and spacecraft, so we may need to agree to disagree on that.

Yup, understood. I'm primarily interested in supporting seeing solar-system bodies in solar imagery, which means that accuracy typically doesn't need to be better than ~100 mas. It's critical that there be a way to get the highest accuracy, but I feel that it's useful to provide additional tools for use cases where a compromise in accuracy may be acceptable to achieve other benefits. In this case, it's the increase in usability due to integration into the transformation framework. See #6068 for a proposed compromise in accuracy to significantly increase performance.

in the context of this PR, the one in Fig. 2 of the SOFA tools documentation would be AstrometricGCRS, is that correct?

No. If I understand the SOFA documentation correctly:

  • "BCRS" in Figure 1 is equivalent to AstrometricICRS, with the treatment of space motion
  • "astrometric" in Figure 1 is shifted relative to AstrometricICRS by the observer location (i.e., the parallax correction)
  • "GCRS" in Figure 1 is closest to AstrometricGCRS (not GCRS, because of the lack of space motion in GCRS), except unfortunately the light-deflection calculation will be done before the space motion, which isn't the correct sequence

In any case, are there fewer issues with that specific frame? E.g., can one do without the assumptions about distance for solar system vs stars? (Which I'd really rather not have anywhere!) If so, would it make sense to just only define that one for now, since it can be constructed directly from SOFA routines?

The problem is handling the motion of solar-system bodies. The SOFA transformations are really intended for stars, so I don't know if you can "hack in" solar-system bodies by thinking in terms of parallax and proper motion.

Edit: As pointed out later (#10372 (comment)), I managed to confuse myself when I wrote parts of the above comment. This PR's AstrometricGCRS does in fact do the corrections in the correct order, so it is consistent with SOFA's GCRS.

@mhvk

mhvk commented May 21, 2020

Copy link
Copy Markdown
Contributor

@ayshih - thanks for explaining! I do feel rather at a loss here, and it helps to understand the goals. I should add that I'm rather weary of adding things which are not precise - I think the goal is that we always provide a best-of-class precision - and am super-happy that for Time at least we are at the level that it can be used for pulsar timing. And while I can see an eventual argument for allowing one to trade precision for performance, I think we should first ensure the code is reasonably close to maximum performance.

Anyway, for the present PR, I guess it does not cover the case I thought it did. The main question would seem to be whether one can do better, and, if not, whether it is better to have a good-enough approximation or whether one should just refer to other packages which do it right.

@mkbrewer

Copy link
Copy Markdown
Contributor

But now having looked a bit at the PR, I got a bit confused. In #10230, the Astrometric frame that was discussed was, at least as I understood it (#10230 (comment)), a frame that sits between BCRS and GCRS in the chain of transformations, and really is just was something that would only need calls to SOFA routines to calculate.

The main issue here is that SOFA is only concerned with extra-solar system bodies (i.e. stars). You can see that in Figure 1 of the SOFA Astrometry Tools documentation as there is no mention of light time delay (planetary aberration, if you will). The same difficulty exists when it comes to gravitational light deflection where SOFA assumes that there is no difference in the direction from the target body to an observer on Earth than from the target body to the Sun. (i.e. the target body is at a distance such that the Earth and the Sun are coincident from its point of view, see issue #10291). In order to provide planetary ephemerides, it is necessary for astropy to meld an ephemeris package such as the one provided by JPL onto the SOFA package. An essential part of this melding is to provide the light time delayed BCRS position of the solar system target body as seen from the observer at any given time. The astrometric position of the target body is then calculated by subtracting the BCRS position of the observer from the light time delayed BCRS position of the target body.

In an analogous manner, the proper way to provide the astrometric position of a non-planetary solar system body (comet, asteroid, spacecraft) is to provide an ephemeris for this body along with a method for calculating the light time delayed position of the body from the point of view of the observer. This ephemeris could be a simple two body propagator that works from a given set of orbital elements or, as in the case of this PR, simply a position and velocity at a given time that is then time delayed by linearly offsetting the given position by using the velocity.

It is important to note here that the act of providing the time delayed position (adding planetary aberration) isn't actually a function of the astrometric frame transformation at all. This is, in fact, analogous to adding space motion in the above referenced Figure 1 and so is, strictly speaking, the step before this where the ICRS catalog position of a star or the ICRS instantaneous ephemeris position of a solar system body is transformed to the BCRS position of the target body at the time of observation from the point of view of the observer. The astrometric frame transformation is then performed by the subtraction of the observer's position in the BCRS from the BCRS position of the target body. This is why the BCRS to astrometric transformation shown in Figure 1 only includes parallax.

That said, the ICRS to BCRS step is an important step in the ICRS->BCRS->astrometric frame transformation, so that functionality needs to be provided. I'm just saying that the paradigm for this should be some sort of ephemeris routine like get_body() rather than anything inherent in the astrometric frame transformation which is why I suggested some sort of apply_light_time() function analogous to apply_space_motion() for instantaneous ICRS coordinates of non-planetary solar system bodies be implemented. A get_non-planetary_body() function might be a better analogy here, but then you'd need a way to define such a non-planetary body which would include its instantaneous ICRS ephemeris/position and a method for providing the light time delayed BCRS position.

Now I think that I understand that in principle, "astrometric" can make sense for other systems, but just so I get things right: in the context of this PR, the one in Fig. 2 of the SOFA tools documentation would be AstrometricGCRS, is that correct?

I assume that you meant Figure 1 here. There is no such thing as "AstrometricGCRS". For a topocentric observer, the position of the observer wrt the Earth's geocenter is expressed in the GCRS frame, but this frame is aligned with the BCRS and, of course, aberration and gravitational light deflection aren't an issue for the observer, so this is actually the BCRS position of the observer since the position of the Earth's geocenter is also expressed in the BCRS.

I was referring specifically to the implementation in this PR, where AstrometricGCRS is "simply" adding proper motion and light-travel-time delay effects to GCRS coordinates.

There seems to be some misunderstanding here. The GCRS coordinates of the target body include stellar and dirunal aberration and gravitational light deflection, so any frame transformation from astrometric to GCRS will include these effects regardless of any velocity present in the target body's BCRS position.

"GCRS" in Figure 1 is closest to AstrometricGCRS (not GCRS, because of the lack of space motion in GCRS), except unfortunately the light-deflection calculation will be done before the space motion, which isn't the correct sequence.

I'm afraid that you lost me here. Space motion is applied during the transformation of the star's ICRS catalog position to the BCRS position at the time of the observation. It is this BCRS position that is subject to gravitational light deflection during the transformation from BCRS to GCRS, so that order is correct.

As a side note, it appears to me that the numerical instability in near parabolic elliptical orbits that @astrojuanlu referred to in the discussion of issue #10230 only occurs when the true anomaly is near 180 degrees as seen in this plot from poliastro/poliastro#475

comet_instability

It makes sense that there would be difficulty in closing the orbit at aphelion in very near parabolic situations, but for the most part comets in these types of orbits are well beyond the range of observation near aphelion so there is likely no practical situation that astropy would be presented with that would preclude the use of a simple two body propagator to provide an ephemeris constructed from orbital elements.

Sorry for the "wall of text", but some things need a lot of exposition to be as absolutely clear as possible.

@mkbrewer

Copy link
Copy Markdown
Contributor

One simple way to handle this might be to add some sort of flag attribute to SkyCoord() objects to indicate what type of object (star, planet, non-planetary solar system object) they are so that you can tell what method to use when transforming from ICRS to BCRS. (i.e. apply_space_motion, do nothing, apply linear time delay).

@ayshih

ayshih commented May 21, 2020 •

Copy link
Copy Markdown
Contributor Author

Sorry, I think you got confused by some of the terminology since backticks are easy to miss and get removed by default when quoting:

  • When I wrote AstrometricGCRS (using backticks), I was referring to the frame that gets created in this PR by calling AstrometricFrame(observer=GCRS(...)), not necessarily any standard astrophyiscal frame.
  • When I wrote GCRS (using backticks), I was referring specifically to Astropy's GCRS class, which doesn't function the same as standard GCRS (in terms of transformations).

I'll be explicit in the following responses, including adding qualifiers to quoted excerpts:

That said, the ICRS to BCRS step is an important step in the ICRS->BCRS->astrometric frame transformation, so that functionality needs to be provided. I'm just saying that the paradigm for this should be some sort of ephemeris routine like get_body()

Yes, having a routine for a ICRS->BCRS conversion that works on both stars and solar-system bodies is a solution. That's equal to the transformation in this PR from Astropy's ICRS to this PR's AstrometricICRS. For reasons, I'm pushing for it to be part of the transformation framework, but it doesn't have to be.

I was referring specifically to the implementation in this PR, where [this PR's] AstrometricGCRS is "simply" adding proper motion and light-travel-time delay effects to [Astropy's] GCRS coordinates.

There seems to be some misunderstanding here. The GCRS coordinates of the target body include stellar and dirunal aberration and gravitational light deflection, so any frame transformation from astrometric to GCRS will include these effects regardless of any velocity present in the target body's BCRS position.

Astropy's GCRS currently includes stellar aberration and gravitational light deflection, but does not include space motion when converting a coordinate from Astropy's ICRS to Astropy's GCRS.
(get_body() does not count.) It does not function the same as standard GCRS. This PR's AstrometricGCRS adds space motion to Astropy's GCRS and gets close to standard GCRS (but not quite, see below). That's why I said the transformation from Astropy's GCRS to this PR's AstrometricGCRS would do nothing for zero body velocity (and hence no space motion). Neither end of that transformation is the same as SOFA's "astrometric".

"GCRS" in Figure 1 is closest to [this PR's] AstrometricGCRS (not [Astropy's] GCRS, because of the lack of space motion in [Astropy's] GCRS), except unfortunately the light-deflection calculation will be done before the space motion, which isn't the correct sequence.

I'm afraid that you lost me here. Space motion is applied during the transformation of the star's ICRS catalog position to the BCRS position at the time of the observation. It is this BCRS position that is subject to gravitational light deflection during the transformation from BCRS to GCRS, so that order is correct.

Again, Astropy's GCRS is not the same as SOFA's GCRS when transforming from ICRS because of the lack of space motion. It essentially assumes no proper motion in the implicit ICRS->BCRS step. This PR's AstrometricGCRS adds space motion to Astropy's GCRS, but because of how it is implemented, it performs the linear-motion propagation to coordinates after gravitational light deflection has already been applied. I know that will make you sad.

Honestly, this last detail might make me scrap this PR. It's one thing to make an assumption of linear motion and accept the resulting inaccuracy. It's a very different thing to do the calculations in the wrong order...

Edit: As pointed out later (#10372 (comment)), I managed to confuse myself when I wrote parts of the above comment. This PR's AstrometricGCRS does in fact do the corrections in the correct order, so it is consistent with SOFA's GCRS.

@mhvk

mhvk commented May 21, 2020

Copy link
Copy Markdown
Contributor

@mkbrewer, @ayshih - thanks both for the explanations!

@ayshih - astropy's GCRS not being consistent with that of SOFA would seem to be a bug, and quite a bad one at that given that I think we are meant to be consistent with SOFA. If so, it would be good to separate it from this PR (though it sounds like the PR might have to rely on solving it to be viable...). Could you do that? I think copying over the relevant parts of your comment above would do.

p.s. My guess is that this may be a bit of a historical accident, since the whole coordinate package was initially written without proper motions, and this may have been simply missed when we added support for those.

@mkbrewer

Copy link
Copy Markdown
Contributor

Astropy's GCRS currently includes stellar aberration and gravitational light deflection, but does not include space motion when converting a coordinate from Astropy's ICRS to Astropy's GCRS.

Yes, I noticed that icrs_to_gcrs() in icrs_cirs_transforms.py does not include space motion. I just assumed that the standard way to transform a star from ICRS to GCRS was to use apply_space_motion() before transforming to GCRS.

This PR's AstrometricGCRS adds space motion to Astropy's GCRS, but because of how it is implemented, it performs the linear-motion propagation to coordinates after gravitational light deflection has already been applied.

I still don't understand this. The first thing that you do in this PR is to transform the target body's position to ICRS, so it seems to me that, if given a GCRS position, that would remove aberration and gravitational light deflection before proceeding to apply space motion or light time delay. But, not being an astropy expert, maybe I'm missing something.

@ayshih

ayshih commented May 22, 2020

Copy link
Copy Markdown
Contributor Author

astropy's GCRS not being consistent with that of SOFA would seem to be a bug

Yes, I noticed that icrs_to_gcrs() in icrs_cirs_transforms.py does not include space motion. I just assumed that the standard way to transform a star from ICRS to GCRS was to use apply_space_motion() before transforming to GCRS.

I guess the question is: If a user transforms in Astropy from ICRS to GCRS, should they expect the transformation to automatically include space motion – which one might expect from how the frames are defined in SOFA – or should they be expected to separately call apply_space_motion() with the appropriate argument? If the former, then it's a bug in the transformation, although how to fix it is non-trivial. If the latter, then it's really just the responsibility of documentation.

This PR's AstrometricGCRS adds space motion to Astropy's GCRS, but because of how it is implemented, it performs the linear-motion propagation to coordinates after gravitational light deflection has already been applied.

I still don't understand this. The first thing that you do in this PR is to transform the target body's position to ICRS, so it seems to me that, if given a GCRS position, that would remove aberration and gravitational light deflection before proceeding to apply space motion or light time delay. But, not being an astropy expert, maybe I'm missing something.

Ha! You're absolutely right. I had a total brain misfire earlier. Too many meetings today. I walk all the transformations through ICRS and do the propagation there, so that should unwind all the corrections and reapply them in the correct order. Whew.

@ayshih

ayshih commented May 22, 2020 •

Copy link
Copy Markdown
Contributor Author

If the latter [that users's should be expected to separately call apply_space_motion()], then it's really just the responsibility of documentation.

Well, documentation and maybe emitting a warning when trying to transform a coordinate with velocity out of ICRS if it doesn't look like it's been propagated (e.g., if it's not in a SkyCoord with a non-J2000.0 obstime)...

@mkbrewer

Copy link
Copy Markdown
Contributor

I guess the question is: If a user transforms in Astropy from ICRS to GCRS, should they expect the transformation to automatically include space motion

An obvious way to automate this would be for transforms from ICRS to GCRS, CIRS or observed to pass through the astrometric frame. Of course, this means that the astrometric frame would need to be very robust in knowing what sort of object it was dealing with in order to take the appropriate action.

@ayshih

ayshih commented May 26, 2020

Copy link
Copy Markdown
Contributor Author

See #10404 for a suggested approach that implements BCRS (in SOFA's Figure 1) and allows a user to define how the ICRS<->BCRS transformations should behave.

@nstarman nstarman mentioned this pull request May 17, 2021
2 tasks
@github-actions

github-actions Bot commented Sep 3, 2021

Copy link
Copy Markdown
Contributor

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.

@github-actions github-actions Bot added the Close? Tell stale bot that this issue/PR is stale label Sep 3, 2021
@github-actions

github-actions Bot commented Sep 4, 2021

Copy link
Copy Markdown
Contributor

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.

@pllim

pllim commented Sep 4, 2021

Copy link
Copy Markdown
Member

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

@github-actions

github-actions Bot commented Oct 5, 2021

Copy link
Copy Markdown
Contributor

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

Close? Tell stale bot that this issue/PR is stale closed-by-bot Closed by stale bot coordinates

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants