Speed up EarthLocation.get_gcrs - #11073

This speed up is amazing! In that call to get_body(), are you using the JPL ephemeris or the built in one? If the later, I'd like to see the timing for the JPL ephemeris also as it has less overhead.

I am very pleased that you have come around to adding this additional method to the EarthLocation class. This is exactly what I wished for in the first place. I would, however, request that you rename it get_posvel() and make it public. That way the user can easily retrieve the ground based observer's position and velocity in any of CIRS, TETE or TEME by simply calling this method with gcrs_to_refset to an identity matrix as is done in the calls to erfa.c2tcio() in the relevant intermediate rotation transforms. In fact, you could make gcrs_to_ref optional with the identity matrix as the default.

I don't think your calculation of the velocity in the target frame is unclear, in fact, I think it is quite innovative

With this change, we should revisit the intermediate rotation transforms for CIRS and TETE as this new method makes the helper function that you added there superfluous.

There does seem to be a slight issue here though:

Expected:
    <SkyCoord (GCRS: obstime=2014-09-22 23:22:00.000, obsgeoloc=(3949481.68990863, -550931.91188162, 4961151.73733451) m, obsgeovel=(40.15954083, 287.47876693, -0.04597867) m / s): (ra, dec, distance) in (deg, deg, AU)
        (136.91116209, 17.02935409, 5.94386022)>
Got:
    <SkyCoord (GCRS: obstime=2014-09-22 23:22:00.000, obsgeoloc=(3949481.69182405, -550931.91022387, 4961151.73597633) m, obsgeovel=(40.159527, 287.47873161, -0.04597922) m / s): (ra, dec, distance) in (deg, deg, AU)
        (136.91116253, 17.02935396, 5.94386022)>

One odd thing is that in solarsystem.rst there are two instances of retrieving the position of Jupiter. One using the built in ephemeris and the other using the JPL ephemeris. Both of these use the same location and time and yet the values of obsgeopos and obsgeovel are different. This PR's test using the built in ephemeris results in these values being closer to those shown in the JPL test.

In EarthLocation.get_posvel() it might be easier for the user if the second argument was just 'cirs', 'tete' or 'teme'. Then if the third argument is a BPN rather than the default, this could be passed to tete_to_itrs_mat() in the case of 'tete' to yield the GCRS position and velocity required by the TETE intermediate rotation transforms.

@mkbrewer - I made a private method somewhat on purpose, as I think it helps separating getting the speed-up from how to expose a more flexible get_posvel to the user.

In an ideal world, I feel we should not need any of these methods as the transformation machinery itself would be fast enough; from that perspective, keeping everything in intermediate_rotation_transformations might actually make more sense than having a method on EarthLocation (indeed, I started on that yesterday and this morning had just decided to revert to that, but concluded again it was less good).

Anyway, let's do the speed-up first!

The failure is very puzzling: I can understand a difference in obsgeovel, but where does this large difference in obsgeoloc come from? Investigating...

Then in the TETE<->GCRS transforms, we could simply have:

    loc_gcrs = tete_frame.location.get_posvel(tete_frame.obstime, 'tete', gcrs_to_ref=rbpn) 

Hmm, I cannot reproduce that locally. I suspect the output is fairly old, and the position already did not match, but is OK for float comparison, while the velocity now changed a little (probably for the better...) and is triggering the failure.

As you noticed, the JPL based example is much closer, even though of course obsgeoloc and obsgeovel do not depend on the ephemeris. Looking through the history, I see that that example was updated more recently, so something else must have (slightly) changed - maybe the IERS B tables?

How about trying it using the JPL ephemeris then? If you do, please time it for both the master and this PR.

OK, I updated the docs with the new values (which are within tolerance). As for timings with JPL:

solar_system_ephemeris.set('de432s')
# Same el and t as on top
%timeit get_body('mars', t, el)
# master  47.9 ms ± 214 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
# current PR 14.7 ms ± 529 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

So, speed-up is very similar (since time was dominated by get_gcrs_posvel....)

Hmm, I thought the JPL lookups would be quite a bit faster than the built in calculations. Perhaps, it's just the overhead associated with file I/O which would be mitigated for an array of times.

Since this results in a slight change in position, I'd appreciate it if you could run this script and post the output:

skyfield_topo_test.py.txt

I haven't converted apparent position to use the TETE frame there, but you are welcome to do that.

Whoops, wrong script.

astropy_topo_test.py.txt

astropy_topo_test_output.txt

Thanks. I'll compare this output to my program.

Woo-Hoo! Everything is now dead on to one digit in the last decimal place in RA/Dec and Alt (3.6 uas) and the maximum offset in Az is 3 decimal places (10.8 uas).

For the most part the distance is good to 1 mm. There are just a few samples that are off by 1.5 - 2 mm.

I'd need to add another decimal place to my output to see any real difference here.

@mkbrewer - for a proper comparison, I wondered if one should transform to TETE rather than use _apparent_position_in_true_coordinates - as the latter infers location from obsgeoloc and obsgeovel rather than use the location directly, so may introduce small rounding errors. But I ran a quick check with a suitably modified version of your code and there are no differences at all in the numerical output.

Yeah, I said that you were welcome to do that. I just wasn't real sure of myself on the proper syntax is all.

I profiled the conversion from GCRS to AltAz:

        1    0.000    0.000    0.364    0.364 intermediate_rotation_transforms.py:171(gcrs_to_cirs)
        1    0.000    0.000    0.096    0.096 intermediate_rotation_transforms.py:45(gcrs_to_cirs_mat)
        1    0.000    0.000    0.002    0.002 intermediate_rotation_transforms.py:50(cirs_to_itrs_mat)
        1    0.000    0.000    0.001    0.001 intermediate_rotation_transforms.py:96(get_location_gcrs)
        1    0.000    0.000    0.262    0.262 icrs_cirs_transforms.py:180(gcrs_to_gcrs)

Looks like most of the time was spent in gcrs_to_gcrs(). Since I don't see any GCRS<->ICRS calls, I guess this is just how long it took to decide that the location and times were equal.

No problem! Here's the output (just in case: code is missing the import of TETE).

astropy_topo_test_output3.txt

Here are three sets of plots.

Astropy_10994-11073_Comparison.pdf

The first in each set is from the last data that I got from you. This was from PR #10994 before the CIRS attributes were changed from GCRS position and velocity to EarthLocation. The second in each set is from this PR.

There is still a very small residual in RA/Dec. This is entirely in the astrometric position as the apparent position tracks it exactly.

To me, the second set of plots is the most interesting. This is the residual in the distance to the Moon. The ~+/-2 cm oscillation and the jumps seen previously are gone. Since all of the samples include the same number of ephemeris lookups, the variation in distance with time of observation from the ephemeris should be smooth. This implies that the oscillation and jumps were somehow associated with obsgeopos which is subtracted off to form the astrometric position. There is still a slight diurnal variation in the distance residual which I attribute to drift within your relatively wide (3 meter) window for acceptance of light time delay.

The way that I calculate light time delay is quite different. I do the ephemeris lookups sequentially in time using the light time delay for each time as the starting guess for the next time. In this way, the light time delay stays very well converged throughout the set of observation times so that there is no drift within an acceptance window.

This is beautiful. The speed up is one thing, but the fact that we have narrowed down the remaining discrepancies with @mkbrewer's code to the FiniteDifferenceTransform's poor precision in velocity is particularly pleasing.

Given that, I have 1 minor request before merging, which is to tighten the tolerances in test_aa_high_precision (https://github.com/astropy/astropy/blob/master/astropy/coordinates/tests/test_intermediate_transformations.py#L683).

They're currently set to 0.5mas, but the plots @mkbrewer added clearly show we are doing better than that. I'd suggest reducing them to 0.002 mas, which should prevent anyone making changes that degrade our precision down the road.

Perhaps we can also add a target distance of 376252883.247 meters as well? This value is taken from MKB's code and I this PR shows agreement to 11 mm, so perhaps a tolerance of 15mm is appropriate.

Here's the latest from my program with one extra digit of precision.

TARGET_AZ, TARGET_EL, TARGET_DIST = 15.0326735105*u.deg, 50.3031101339*u.deg, 376252883.2473*u.m

Yes, good point; I tightened the tests, using the new numbers as well. Also slightly rewrote the test to compare the numbers directly (but then added rtol=0 since otherwise the tests are not as accurate as implied; let of a problem if comparing with 0...)

But, if you want that tolerance in azimuth, you are going to need to multiply the difference by cosine elevation.

This is because the Moon crosses the lower meridian at -80.5 degrees, so the azimuth coordinate offset is quite a bit larger than the actual offset on the sky.

(codecov failure is unrelated)

Somehow, it squeaked by. I see a residual of 2.16 uas in azimuth for that point.

Thanks @mhvk. Happy to merge this now unless you think we should wait for others to check it over.

Oh, and that 15 mm tolerance for distance is a bit wide. The actual actual maximum residual is 1.3 mm.
Perhaps you misread the Y axis scale there.

Checking numbers:

(moon_aa.az-TARGET_AZ).to(u.uas)
<Angle -1.10632783 uas>
((moon_aa.az-TARGET_AZ)*np.cos(TARGET_EL)).to(u.uas)
<Angle -0.70664041 uas>
(moon_aa.alt-TARGET_EL).to(u.uas)
<Angle 0.07351559 uas>
(moon_aa.distance-TARGET_DISTANCE).to(u.mm)
<Distance -0.01775334 mm>

So maybe 2 uas is OK and for distance we can use 2 mm (which is about 1uas * distance). I'll still change the latter.

In it went - thanks, @mkbrewer and @StuartLittlefair!

That's odd. Here's what I have from the file you posted.

Astropy     15.0326735099
MKB         15.0326735105
Delta (mas) -0.0021600002

Perhaps we are down to the point where round off introduces a bit of random noise.

One difference is that in the test case, the observer's location gets passed in already in get_body, so its location is dealt with a bit earlier. This may make the rounding errors differ.

@mhvk @StuartLittlefair This is a followup regarding the improved accuracy of this PR. While working with Brandon Rhodes on polar motion, we discovered that the remaining issue that was causing the residual in ICRS RA/Dec to be +/-0.2 mas and the residual in moon distance to be +/-40 cm was that he is using the IERS 2010 ellipsoid rather than the WGS-84 ellipsoid that my site coordinates are referred to. After using a bit of a hack to redefine his ellipsoid as WGS-84, these remaining residuals went away. I then noticed that the noise in the moon distance was +/-1 mm again just as it was for this PR. This led me to realize that the noise was originating on my side.

I was using the method to build two part dates recommended in apco() which is to set the first part to J2000 and the second part to days since. It turns out that this method has the drawback of causing additional round off noise in the JPL ephemeris. When I switched my method to date and day fraction (as is also recommended in era00(), which I hadn't noticed before), I got a large improvement. So here's the revised comparison of my program with this PR:

MKB-Astropy_11073_Comparison.pdf

Most of the remaining residual is due to the need for yet another decimal place in the Astropy output!

Here also is the new comparison of my program with the latest version (still a work in progress) of Skyfield after hacking it to use the WGS-84 ellipsoid:

MKB-Skyfield_WGS84_Comparison.pdf

Most of the remaining residual there is due the fact that neither my program nor Astropy currently supports deflection of light due to the gravitational field of the Earth itself (see issue #10843).

Fantastic!!! (And am happy to run the tests again with different precision...)

OK. If you are willing, this is one decimal above machine precision.

astropy_topo_test.py.txt

Here's the output:
astropy_topo_test_output4.txt

@mhvk Thanks for posting that. It led me to yet another subtle bug in my program that led to ~40 ns of round off jitter in my times. After correcting that, here are the new plots:

MKB-Astropy_11073_Comparison.pdf

The average residuals are now 61 nas (yes, nano-arcseconds) in ICRS RA, 29 nas in ICRS Dec, 66 nas in Apparent RA, 31 nas in Apparent Dec, 5 nas in Azimuth, 0.2 nas in Elevation and 28 microns in distance.
The average LST is off by 0.54 usec which is very close to what I'd expect from the missing TIO locator.

Here are new target values for 2017-04-06 00:00:00 UTC from my program:

TARGET_AZ, TARGET_EL, TARGET_DIST = 15.032673509886*u.deg, 50.303110133928*u.deg, 376252883.247239*u.m

Although, you might as well use your own.

Wow!! Of course, we're now both so similar in part because (I think) we use the same SOFA/ERFA routines. I'm starting on comparisons with CALC11 - ignoring ocean loading for now (which I didn't know was a thing until a week ago...).

I adjusted the test - realizing that at first it failed because I had not included the location in the time object (and this matters slightly for transformation to other time scales...). See #11103.

Hi all, I've been eagerly following these developments, as I'm in the process of switching to Astropy for our radio interferometer delays within the katpoint package. I'm also continuously comparing against CALC11 and trying to squeeze out the maximum reasonable precision.

@mhvk, you would probably be interested in https://gitlab.com/ska-telescope/almacalc/. It's not quite ready for wider release (for one, the license and copyright are probably wrong) but you are welcome to give it a spin 😄 And I would very much like to compare notes.