Speed up GCRS<->(TETE, CIRS) by smarter location conversions by mhvk · Pull Request #11069 · astropy/astropy

mhvk · 2020-11-20T20:06:15Z

Following the suggestion by @mkbrewer in #11067, this speeds up TETE<->GCRS considerably by using the rbpn matrix also for getting the sidereal time for the conversion of the observer location from ITRS to TETE.

fixes #11067

~~Note: need to re-unify code so for now just a draft so we have something to look at.~~ (done)

EDIT: quick timing:

from astropy.coordinates import SkyCoord, EarthLocation
from astropy import units as u

tete = SkyCoord(10*u.deg, 11*u.deg, frame='tete', location=EarthLocation(11*u.deg, 12*u.deg))
%timeit tete.transform_to('gcrs')
# current master
# 16.6 ms ± 1.13 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
# this PR
# 6.51 ms ± 609 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)

mhvk · 2020-11-20T21:10:00Z

Sorry, this is definitely not right: obsgeoloc and obsgeovel become all 0; I'm doing something stupid.

mhvk · 2020-11-20T21:12:18Z

More generally, a sequence of matrix multiplications cannot make the velocity non-zero, so obsgeovel is wrong even if I get obsgeoloc to be right.

mkbrewer · 2020-11-20T21:14:57Z

That check against Horizons is off by 288 mas. This suggests that there is something wrong at the polar motion level.

mkbrewer · 2020-11-20T21:17:10Z

Oh, never mind. I see you have already discovered the error.

This helps performance of the TETE/CIRS<->GCRS transforms.

mhvk · 2020-11-20T21:53:43Z

OK, I get tests to pass again by reverting to what used to be done in get_gcrs_posvel() before we had velocities in our coordinates -- just calculating the velocity from the earth's rotation rate. So, this avoids the finite difference calls to get the velocity (which is why the speed-up is as large as it is: it not just halves the number of calls to pnm06, but quarters it!)

mhvk · 2020-11-20T21:57:01Z


-def get_location_gcrs(reference, gcrs_to_reference_matrix):
-    """Create a GCRS frame at the location and obstime of a given reference.
+OMEGA_EARTH = u.Quantity(7.292115855306589e-5, 1./u.s)


Actually, still defined in earth.py, although it is not used there. I wonder if we should revert to the old code there too...

If we do revert to the old and much faster way, we can also revert out the "include_velocity" argument....

And if we do that, we perhaps can instead allow passing in a matrix in get_gcrs_posvel - though let's leave that for another PR

Yes - let’s leave earth.py alone for now.

mhvk · 2020-11-20T22:27:57Z

Ah, this is neat: with this change, the GCRS->CIRS->ICRS and GCRS->TETE->ICRS paths again give consistent results within 1 mm! I added a commit with the tightened test.

mhvk · 2020-11-20T22:28:46Z

-                obsgeovel=gcrs_cart.differentials['s'].to_cartesian())
+    itrs_cart = reference.location.get_itrs(obstime).cartesian
+    obsgeoloc = itrs_cart.transform(matrix_transpose(gcrs_to_itrs_matrix))
+    vel_x = -OMEGA_EARTH * obsgeoloc.y


@mkbrewer - is there an ERFA routine that does this (maybe even in a better way)?

Hmm, looking through why we changed it in earth.py, I came back to #6699, which planned to do:

# CIRS velocity rotation_vector = CartesianRepresentation(0/u.s, 0/u.s, OMEGA_EARTH) cirs_cart = itrs.transform_to(cirs_frame).cartesian cirs_vel = rotation_vector.cross(cirs_cart) # Transform to GCRS (rotation) obsgeovel = cirs_vel.transform(matrix_transpose(gcrs_to_cirs_mat(obstime))) return obsgeoloc, obsgeovel

But then in #6861 it was changed to its present state. @StuartLittlefair - both were your PRs, so you probably remember better what is right/wrong with it.

This was removed when velocities were added to coordinates purely because it was simpler.

mkbrewer · 2020-11-20T22:48:59Z

The ERFA routine is pvtob(), but it does the geodetic to geocentric conversion also which is not necessary with an EarthLocation.

mkbrewer · 2020-11-20T22:56:32Z

Hmm, I wonder why my comment was automatically deleted there. Anyway, yes, we both thought of that at basically the same time.

mhvk · 2020-11-20T22:57:27Z

Ah, from that routine, it looks like I should be calculating the velocity in CIRS, as was done in #6699, so this needs an update.

mhvk · 2020-11-20T22:58:08Z

It is because I deleted it, following your example. Sorry!!

mhvk · 2020-11-20T23:00:22Z

If we need to go to CIRS to add the velocity, maybe it makes more sense after all to go for your second option, i.e., transform the ITRS location to GCRS via CIRS also for TETE positions.

mkbrewer · 2020-11-20T23:02:13Z

I think you are fine there. You are just calculating the velocity after the final rotation. Otherwise, you'd need to rotate the velocity also. See apco() for an example of this. In fact, we could optimize ErfaAstrom.apco() a bit by using your transform there also.

mkbrewer · 2020-11-20T23:05:38Z

btw, in fact, pvtob() works in either CIRS or TETE, see note 3) there.

mhvk · 2020-11-20T23:10:12Z

OK, that's good, so the plan would be to get the location it in ITRS, transform it to either CIRS or TETE using one matrix, then add the velocity as I do now, and finally transform to GCRS with the second matrix (which applies directly to the velocity as well). That should be a straightforward change!

mkbrewer · 2020-11-20T23:14:21Z

I don't see the point of making that change. Unless you actually want the velocity in CIRS or TETE, you might as well calculate it in GCRS. That saves one extra rotation.

mhvk · 2020-11-20T23:27:07Z

But I now assume the velocity is perpendicular to the GCRS Z axis, while I think it should be perpendicular to the real Earth's polar axis (and the matrix looked like it had rotations around x and y, but perhaps I misinterpreted c2i06a)

mkbrewer · 2020-11-20T23:32:22Z

Well, it is perpendicular to the CIP. The transformation to GCRS just aligns the CIP with the ICRF. That's why apco() rotates both the position and velocity into the GCRS. Remember, in that transformation we are going from of date to the mean equator and equinox of J2000.

mhvk · 2020-11-21T02:14:30Z

It's not quite perpendicular, though, and for some applications that small difference may matter (which was the reason for #6699 and the subsequent #6861)

from astropy import units as u
from astropy.time import Time
from astropy.coordinates import CIRS, GCRS
obstime = Time('J2010')
cirs = CIRS(x=10000*u.km, y=0.*u.km, z=0.*u.km, representation_type='cartesian', obstime=obstime)
cirs.transform_to(GCRS(obstime=obstime)).cartesian
# <CartesianRepresentation (x, y, z) in km
#    (9999.99496686, -0.00014031, -10.03308023)>

mkbrewer · 2020-11-21T02:39:53Z

That's actually a quite interesting demonstration of the stability of the CIO. In ten years it only moved by 14 cm while the equator moved by 10 km.

mkbrewer · 2020-11-21T02:45:20Z

Since it seems that I'm not going be able to convince you, here's an experiment that you could do. First calculate the velocity in the CIRS frame and then rotate the position and velocity into the GCRS frame. Then recalculate the velocity from the GCRS position and compare the two.

mhvk · 2020-11-21T03:10:16Z

My lack of being convinced is just because of @StuartLittlefair's old #6699... But strengthened by the fact that with the current PR, loc.get_gcrs_posvel() and my quick routine give inconsistent velocities. If I apply the velocity in the CIRS or TETE frame, then they become consistent. I've just pushed a new commit that includes tests for this - do have a look, since if this is wrong, it affects people trying to get accurate radial velocities to hunt for planets...

mhvk · 2020-11-21T03:17:00Z

If we go this route, we should make sure OMEGA_EARTH is exactly consistent with the value used in ERFA's pvtob (unfortunately, not exported in their .h file). It might also be a good idea to rewrite EarthLocation.get_gcrs_posvel() to use the same method - I had not really realized how by just using the standard transformations, the velocity requires multiple evaluations. Though perhaps it is better to try to rewrite the transformation mechanism itself so that it uses knowledge of how velocities change.

mkbrewer · 2020-11-21T06:50:18Z

OK, you're right. The frame rotates, but the vectors don't. Therefore the spin picks up a z component in the new frame. I was erroneously thinking that the vectors stayed aligned with the CIP. Sorry about that.

Omega is defined in pvtob.c:

/* Earth rotation rate in radians per UT1 second */
   const double OM = 1.00273781191135448 * D2PI / DAYSEC;

This is just a sidereal second converted to radians.

StuartLittlefair · 2020-11-21T09:33:25Z

A lot happens when I’m asleep! It looks like you’ve sorted out all the issues; as you say velocity should be added in CIRS or TETE and then rotated from there.

I’m very pleased this also fixes the slight inaccuracies in the ICRS round trip as well as providing a massive speed up.

mkbrewer · 2020-11-21T15:01:36Z

@mhvk You might add that this fixes issue #11067

mkbrewer · 2020-11-21T15:59:48Z

There is a slight discrepancy between SOFA and Astropy regarding OMEGA. Using SOFA's definition of D2PI, I'm seeing:

OM = 7.292115146706980e-05

mhvk · 2020-11-21T17:11:29Z

@mkbrewer - thanks, I think we might as well be consistent with what SOFA has.

I'll make another commit that updates the value and adds a changelog entry about the speed and accuracy improvements.

mhvk · 2020-11-21T19:40:49Z

OK, I changed OMEGA_EARTH to match ERFA (and found where they got the value from...), and also added changelog and what's-new entries. @mkbrewer - for the latter, I mention you specifically, as I wanted to ensure your really helpful advice is recognized! (Normally, people are recognized by commits, but that process is not the only one...)

mhvk · 2020-11-21T19:53:17Z

I put some follow-up items in #11072

mhvk added coordinates Performance no-changelog-entry-needed labels Nov 20, 2020

mhvk mentioned this pull request Nov 20, 2020

More efficient algorithm for deriving observer's GCRS position and velocity in TETE frame. #11067

Closed

mhvk requested a review from StuartLittlefair November 20, 2020 20:08

Speed up location to obsgeoloc/vel transforms.

ce7e1c0

This helps performance of the TETE/CIRS<->GCRS transforms.

mhvk force-pushed the coord-tete-performance-improvement branch from db502ad to ce7e1c0 Compare November 20, 2020 21:50

mhvk commented Nov 20, 2020

View reviewed changes

Tighten comparison between GCRS-ITRS via CIRS and TETE

de54d05

mhvk commented Nov 20, 2020

View reviewed changes

astropy deleted a comment from mkbrewer Nov 20, 2020

Add velocity in CIRS/TETE, not GCRS

49592bf

StuartLittlefair approved these changes Nov 21, 2020

View reviewed changes

StuartLittlefair marked this pull request as ready for review November 21, 2020 09:33

mhvk changed the title ~~Rough check to see if we can speed up TETE~~ Speed up GCRS<->(TETE, CIRS) by smarter location conversions Nov 21, 2020

mhvk removed the no-changelog-entry-needed label Nov 21, 2020

mhvk added this to the v4.3 milestone Nov 21, 2020

mhvk added the precision label Nov 21, 2020

mhvk added 2 commits November 21, 2020 14:42

Ensure our OMEGA_EARTH matches that used in SOFA.

6a00bbf

Add changelog and whats-new entries.

177e7dd

mhvk added the zzz 💤 merge-when-ci-passes Do not use: We have auto-merge option now. label Nov 21, 2020

mhvk mentioned this pull request Nov 21, 2020

Further performance improvements for earth-centred transformations #11072

Open

1 of 3 tasks

mhvk mentioned this pull request Nov 21, 2020

Cannot directly transform output from GCRS to AltAz / need ITRS location? #10996

Open

mhvk merged commit 2c6fd69 into astropy:master Nov 21, 2020

mhvk deleted the coord-tete-performance-improvement branch November 21, 2020 21:08

Conversation

mhvk commented Nov 20, 2020 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

mhvk commented Nov 20, 2020

Uh oh!

mhvk commented Nov 20, 2020

Uh oh!

mkbrewer commented Nov 20, 2020

Uh oh!

mkbrewer commented Nov 20, 2020

Uh oh!

mhvk commented Nov 20, 2020

Uh oh!

mhvk Nov 20, 2020 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

StuartLittlefair Nov 21, 2020

Choose a reason for hiding this comment

Uh oh!

mhvk commented Nov 20, 2020

Uh oh!

mhvk Nov 20, 2020

Choose a reason for hiding this comment

Uh oh!

mhvk Nov 20, 2020

Choose a reason for hiding this comment

Uh oh!

StuartLittlefair Nov 21, 2020

Choose a reason for hiding this comment

Uh oh!

mkbrewer commented Nov 20, 2020 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

mkbrewer commented Nov 20, 2020

Uh oh!

mhvk commented Nov 20, 2020

Uh oh!

mhvk commented Nov 20, 2020

Uh oh!

mhvk commented Nov 20, 2020

Uh oh!

mkbrewer commented Nov 20, 2020

Uh oh!

mkbrewer commented Nov 20, 2020

Uh oh!

mhvk commented Nov 20, 2020

Uh oh!

mkbrewer commented Nov 20, 2020 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

mhvk commented Nov 20, 2020

Uh oh!

mkbrewer commented Nov 20, 2020 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

mhvk commented Nov 21, 2020

Uh oh!

mkbrewer commented Nov 21, 2020

Uh oh!

mkbrewer commented Nov 21, 2020

Uh oh!

mhvk commented Nov 21, 2020

Uh oh!

mhvk commented Nov 21, 2020

Uh oh!

mkbrewer commented Nov 21, 2020

Uh oh!

StuartLittlefair commented Nov 21, 2020

Uh oh!

mkbrewer commented Nov 21, 2020 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

mkbrewer commented Nov 21, 2020

Uh oh!

mhvk commented Nov 21, 2020

Uh oh!

mhvk commented Nov 21, 2020

Uh oh!

mhvk commented Nov 20, 2020 •

edited

Loading

mhvk Nov 20, 2020 •

edited

Loading

mkbrewer commented Nov 20, 2020 •

edited

Loading

mkbrewer commented Nov 20, 2020 •

edited

Loading

mkbrewer commented Nov 20, 2020 •

edited

Loading

mkbrewer commented Nov 21, 2020 •

edited

Loading