to_skycoord method does not set correct frame

[lpsinger]

The to_skycoord method does not set the Astropy frame, so it subsequent transformations yield wildly incorrect results. For example:

>>> from skyfield.api import load
>>> ts = load.timescale()
>>> eph = load('de421.bsp')
>>> moon = eph['moon']
>>> coord = moon.at(ts.utc(2020, 3, 1))
>>> coord.to_skycoord().itrs.earth_location
<EarthLocation (0.43176788, -0.84260876, -0.24923164) AU>

This is way off: the moon is obviously not a large fraction of an AU from the Earth.

[brandon-rhodes]

Currently the object is created with the call:

SkyCoord(representation_type='cartesian', x=x, y=y, z=z, unit=au)

How can Skyfield set not only the three Cartesian coordinates, but also the “center” (as SPICE would call it) so that the object can make explicit to AstroPy that it is measuring from the Solar System Barycenter and not from some other position like the center of the Earth?

[lpsinger]

How can Skyfield set not only the three Cartesian coordinates, but also the “center” (as SPICE would call it) so that the object can make explicit to AstroPy that it is measuring from the Solar System Barycenter and not from some other position like the center of the Earth?

That would be the frame. For example...

from astropy.coordinates import GCRS

SkyCoord(x, y, z, unit=au, representation_type='cartesian', frame=GCRS(obstime=t))

Here, t should be an Astropy time object. Replace GCRS with ITRS etc. as needed.

[brandon-rhodes]

So an AstroPy frame provides not only axes, but a coordinate origin? Could you point me at a table that might map to their complete list — it would help me with two issues:

1. Skyfield can easily create coordinates measured from an observer on the Moon, or on Jupiter. This means that I not only need to know how to generate ICRS and GCRS coordinates, but I'll also need to know how to select frames based on the Moon or Jupiter.
2. Do you know what the correct approach is for locations on the Earth's surface? Because of their velocity around the Earth's center, they are not exactly GCRS positions — first, because, obviously, their origin is offset from the Earth's center by roughly the Earth's radius; and, second, because things like aberration affect them differently because their velocity vector does not point exactly where the Earth's does.

Or should Skyfield simply raise an error for all positions that are not either SSB centered or geocentered, and say “AstroPy can’t represent this position”?

[lpsinger]

So an AstroPy frame provides not only axes, but a coordinate origin?

Yes, this is my understanding.

Could you point me at a table that might map to their complete list

https://docs.astropy.org/en/stable/coordinates/index.html#classes

— it would help me with two issues:

1. Skyfield can easily create coordinates measured from an observer on the Moon, or on Jupiter. This means that I not only need to know how to generate ICRS and GCRS coordinates, but I'll also need to know how to select frames based on the Moon or Jupiter.

Astropy currently lacks frames for planets other than the Earth. It would be safer to return a NotImplemented error rather than returning sky coordinates in an incorrect frame.

2. Do you know what the correct approach is for locations on the Earth's surface? Because of their velocity around the Earth's center, they are not exactly GCRS positions — first, because, obviously, their origin is offset from the Earth's center by roughly the Earth's radius; and, second, because things like aberration affect them differently because their velocity vector does not point exactly where the Earth's does.

I don't. I'm a little fuzzy myself about the distinction between GCRS, CIRS, and ITRS. However, Astropy coordinates generally can have both positions and velocities attached to them. The Astropy example based on your sgp4 library has an example of this.

Tangentially related: I was just noticing (and discussing on the Astropy slack) that this example lacks the correct treatment of leap seconds, and I was thinking about writing an Astropy time scale to fix that.

Or should Skyfield simply raise an error for all positions that are not either SSB centered or geocentered, and say “AstroPy can’t represent this position”?

Yup.

[brandon-rhodes]

Thanks for the links, I think that provides everything I need — I hadn't realized that AstroPy was so tightly limited in the kinds of coordinates it can represent. I should be able to get this implemented within the next version or two of Skyfield. In the meantime, let me know if there's anything you need to work around the limitation — it looks like you probably can get your code working for the moment?

[lpsinger]

I hadn't realized that AstroPy was so tightly limited in the kinds of coordinates it can represent.

Indeed. You can return a bare Representation without any frame information (see, for example, astropy.coordinates.SphericalRepresentation), but that would probably be of limited use to other users.

Interestingly, Astropy does have the JPL planetary ephemerides...

https://docs.astropy.org/en/stable/api/astropy.coordinates.get_body_barycentric_posvel.html#astropy.coordinates.get_body_barycentric_posvel

it looks like you probably can get your code working for the moment?

Actually, there is one other thing, although it is only tangentially related... for my application, I actually prefer calling your SGP4 package directly rather than Skyfield, because I want more control over the interaction with Astropy and Numpy. I am using Astropy for all of the time and coordinate frame conversions, but I am only able to reproduce the ITRS position to an accuracy of about 20m. Would you please take a look at nasa/dorado-scheduling#10?

[brandon-rhodes]

Would you please take a look at nasa/dorado-scheduling#10?

An accuracy of 20m is actually very good for an SGP4 position, which are often exact to only a few kilometers. The difference might arise from their implementations of TEME, as it's a rather ill-specified frame?

[brandon-rhodes]

I've finally had time to sit down and put together a simple solution. Hopefully the code's new caution doesn't break Skyfield for anyone who was counting on getting positions back, even if meaningless/inaccurate. We'll see!