Fix coordinate transformations near Earth

[eteq]

This PR closes #4221. It turns out the problems there were a bit more complicated because there were two intertwined bugs:

• The way distances were being calculated for the ICRS<->CIRS transformations and the CIRS<->AltAz transformations was wrong in a somewhat subtle way. This lead to errors ranging from a few thousand km to as much as an AU, which is primarily important inside the solar system. Because the angular locations are computed using ERFA, they are not substantially affected, only distances.
• When doing the CIRS<->AltAz transformations, computing the distance change requires working with the AltAz's EarthLocation. It turns out that ITRS coordinate corresponding to the EarthLocation was not correctly being set to the obstime of the AltAz, but rather defaulting to J2000. So distance were all computed relative to the J2000 location of the Earth rather than the current obstime. So this lead to distance errors of order 1-2 AU at any time other than J2000. This PR addresses that problem by creating a way to get the ITRS from an EarthLocation with a particular obstime, and then using that in the transformation.

There's actually a lot more in this PR then just those fixes, but it's essentially all tests (new ones or reorganization of old ones), because figuring out the root cause of these issues involved a ton of testing the various transform steps. In any event, these tests should make solar-system scale transformations a fair bit better tested.

Note that I also tried backporting this to v1.0.x - there are some conflicts... nothing major, but I should probably do the backporting once this gets merged so that we can be sure the right things get backported. Backporting to 1.1.x should be clean, though.

@embray @astrofrog - this fix should actually go into both 1.0.x and 1.1.x, so I'm correct that I should milestone it for 1.0.6, right?

cc @bmorris3 @adrn @mhvk

[bmorris3]

Background for those other than @eteq: this bug was discovered when I was trying to compute the position of the moon using jplephem and transforming cartesian coordinates to AltAz.

My test code computes the alt/az position of the moon using PyEphem and compares it to the alt/az position I get from jplephem after transformations with SkyCoords. The test code runs on v1.0.5 and demonstrates the large difference between the answers (i.e., the bug):

PyEphem (Alt, Az): (55.4920153458 deg, 239.031895666 deg)
jplephem (Alt, Az): (34.177292986 deg, 135.374232837 deg)

When I run the same script with astropy running from this branch (@eteq:fix-gcrs-trans), I get RuntimeError: maximum recursion depth exceeded in __instancecheck__.

[mhvk]

If we define obstime=None here (or whatever is necessary to make it use the "default obstime"), you can just do itrs = property(get_itrs) and remove the duplicate definition above.

[eteq]

good point! will do.

[mhvk]

@eteq - so far only looked at the EarthLocation changes -- see the small inline comments. A larger one (probably for a separate issue if at all relevant): EarthLocation was always a bit of a hack to get sidereal times going. Should it just be a ITRS subclass?

[eteq]

@bmorris3 : Hmm, I got a few errors relating to EarthLocation (unrelated to this specific set of changes) when I first ran it, but after tweaking it a bit I got this:

With astropy v1.1.dev13915
PyEphem (Alt, Az): (55.4920153458 deg, 239.031895666 deg)
jplephem (Alt, Az): (57.0399296347 deg, 238.705846051 deg)

which is a lot closer than before. Is there any chance you're running from an old version of my fix-gcrs-trans branch? I was getting recursion errors a while back for some in-progress work, but those should have been fixed... In the meantime, I'll try to take a closer look at your example code to try to figure out the origin of the remaining few degrees discrepancy.

[eteq]

@mhvk - You may well be right that EarthLocation should be re-implemented as ITRS or something like it... There is one important distinction right now, though: ITRS has a specific obstime attached, while EarthLocation does not. So EarthLocation is sort of a "location on the earth at any time", while ITRS is "location in space-time, but written down in a way that's somewhat more convenient for translating to specific parts of the Earth itself". I think this distinction is useful, but it certainly might be worthwhile to take a look at how we might re-work EarthLocation into something that looks more like a frame class.

I think that's better tackled separately for v1.2, though, as it's definitely not just a "bug fix".

[mhvk]

@eteq - OK, new issue on EarthLocation possibly being an ITRS subclass: #4261

[eteq]

Good point @mhvk - just pushed up a new commit doing just what you suggested there.

[embray]

So is this ready to be merged? I'd like to release v1.0.6 ASAP.

[eteq]

@embray Well, it doesn't look like anyone has reviewed the transformation code yet... So we could merge it given the tests are passing and such, but it would probably be better if someone other than me would review the distance calculations... (@adrn or @astrofrog?)

[mhvk]

@eteq - from #4254 (comment), it seemed there was still a problem at the ~degree level. Has that been resolved?

[embray]

It can be delayed for 1.0.7 if need be. Or merged provisionally since it at least superficially fixes the problem...??

[adrn]

Sorry -- won't have time to review in detail before, say, Nov. 1 :(

[eteq]

Alright, I made some more adjustments to @bmorris3's moon test code and also tested this on a branch that includes both this (#4254) and #4255 - see https://gist.github.com/eteq/08dba433082e57b7f2af for what actually yields this:

With astropy v1.1.dev13939
PyEphem (Alt, Az, Dist): (55.4920153458 deg, 239.031895666 deg, 356424848.107 m)
Skyfield (Alt, Az, Dist): (55.32651159 deg, 239.298837015 deg, 0.00238248020632 AU)
Skyfield(de430) (Alt, Az, Dist): (55.3265117716 deg, 239.298837058 deg, 0.00238248020483 AU)
jplephem/Astropy (Alt, Az, Dist): (56.0578669005 deg, 239.043478848 deg, 356423493.806 m)
jplephem/Astropy(ICRS) (Alt, Az, Dist): (56.0596259457 deg, 239.03677597 deg, 356421587.215 m)

Whereas a branch with just this gives the last two as:

With astropy v1.1.dev13936
PyEphem (Alt, Az, Dist): (55.4920153458 deg, 239.031895666 deg, 356424848.107 m)
Skyfield (Alt, Az, Dist): (55.32651159 deg, 239.298837015 deg, 0.00238248020632 AU)
Skyfield(de430) (Alt, Az, Dist): (55.3265117716 deg, 239.298837058 deg, 0.00238248020483 AU)
jplephem/Astropy (Alt, Az, Dist): (57.0399296347 deg, 238.705846051 deg, 357213053.992 m)
jplephem/Astropy(ICRS) (Alt, Az, Dist): (56.0596259457 deg, 239.03677597 deg, 356421587.215 m)

So with that, the difference @mhvk is mentioning between PyEphem and astropy are at the <~half a degree level, which is also the disagreement between Skyfield and PyEphem (which are both by the same author and conceptually similar, but with somewhat different ephemerides). At that level I'm willing to chalk it up to subtle differences in how the various frames are to be interpreted (e.g., the ICRS/GCRS difference in the last two lines are because ICRS->AltAz includes aberration, light deflection by the sun, etc., which of course is not relevant for the moon but might be for objects elsewhere in the solar system).

[eteq]

So unless someone thinks they can merge it within the next, say, day (is that your timescale for 1.0.6, @embray?), we should probably just merge and do a "live" review, as this is a big improvement over master/1.0.5 ...

Of course there's nothing to stop us from merging and then getting a more thorough review when someone has time, with a follow-on PR fixing any major issues.

[embray]

I'm thinking more like the next hour, but yeah :)

[eteq]

This block is in particular could use review/somebody else considering it. The geometry is rather confusing, but most importantly, I'm not certain how numerically stable this solution to the SSA trig problem is. So if anyone else knows a trick for this particular geometric calculation I'd love to hear it.

[eteq]

Hah, ok, that works too @embray - if you want to pull the trigger to make sure this makes 1.0.6, feel free to do so. I'd still appreciate whatever review this can get by someone willing to take the time to wrap their brain around the distance calculations, though!

[embray]

@eteq Maybe open an issue to ensure that this gets further review?

[mhvk]

Agreed that it is better to have something we at least think is correct! I'm happy to have this merged as long as we have a new issue with a very specific test case (the moon seems fine...). Ideally, for that, we also involve @brandon-rhodes.

[embray]

Just to note, when I backported this to v1.0.6 I had to remove some stuff related to galactic coordinates, I think it was?, that wasn't applicable to 1.0.x. But I confirmed that the other tests passed.

[eteq]

@embray - am I right that 60edc62 is the commit that backported this to v1.0.x ? That seems to be the right diff, but the commit message doesn't seem right... That's not really a big deal, except that I wonder if it will confuse some of the backporting tools/log parsers?

But the backporting all looks find in terms of the actual content! Will create a new issue shortly to prompt someone to do a review.

[embray]

Weird about the commit message. I think I might know how that happened. Now I need to make sure the other commit messages on that branch are correct.... I will fix it.