Add topocentric CIRS and fix the CIRS to AltAz transformations - #10994

Hello @StuartLittlefair 👋! It looks like you've made some changes in your pull request, so I've checked the code again for style.

There are no PEP8 style issues with this pull request - thanks! 🎉

Comment last updated at 2020-11-19 17:34:34 UTC

I just want to comment that while I think this is great, I also think this should not make it into a release without the faster, direct ICRS to AltAz transformation.

A new release with slower transformation between ICRS and AltAz to gain some arcseconds for close-by objects would be catastrophic for us who are mostly interested in fast, not terribly precise ICRS to AltAz transforms for FAR away objects.

Totally agree @maxnoe

These test failures look unrelated - a timeout on readthedocs, and the coverage test is failing due to stale code - that I am leaving in place since I intend to use it in a later PR. I'm not sure why the windows build on travis failed to tear down properly.

Ideally, we keep the implementation of ICRS<->AltAz decoupled from this PR (which I have not looked at but suspect will not be super-simple!), and it could be either in parallel or after - since we have quite a while before 4.3, either would seem fine.

I think it is best to leave the ICRS transform to another PR. There’s enough to review here already, although it’s not that complicated in the end given how long I’ve spent on it.

The cost of this precision is a rather severe slowdown in the transformation code - approximately a factor of 10.

I am curious about this slowdown. Is this in the ICRS<->CIRS or CIRS<->AltAz transforms? I don't see it being in the former as not much has really changed there.

This is because the calls to erfa.apcs and calculation of the celestial-to-intermediate matrix are more expensive that the calls to erfa.apci we use in the master branch

This isn't true. Calculation of the celestial to intermediate matrix has always been required as it is needed by erfa.apci. Then there's the fact that erfa.apci ends up calling erfa.apcs anyway.

If the slowdown is in the CIRS<->AltAz transformations, this is because of the call to altaz_coo.location.get_gcrs_posvel() necessitated by your choice of GCRS position and velocity to specify observer location for CIRS frames. While I completely understand the convenience and consistency of this approach, that call does require the costly calculation of the celestial to intermediate matrix. An alternate approach would be to use CIRS position and velocity for the observer location for CIRS frames. Then an EarthLocation.get_cirs_posvel() function could be written for use here that would be much less costly. Using CIRS position and velocity wouldn't be all that inconvenient either. In the GCRS<->CIRS transforms, it is necessary to calculate the celestial to intermediate matrix anyway, so this could be applied to the observer's position and velocity also during these transformations. In the ICRS<->CIRS transforms, the observer's GCRS position and velocity required by the call to erfa.apcs could also be easily obtained by using the celestial to intermediate matrix which is also calculated in erfa_astrom.get().apci() anyway.

@mkbrewer, you are wrong about the cause of the slowdown, but then again so am I. I have profiled the code to see where the time is being sent and it’s not in any calls to get_gcrs_posvel

The way it is coded is slightly slower but the bulk of the slowdown comes from the fact that CIRS->AltAz contains an ICRS self transform to get topocentric CIRS.

So what was a fairly direct ICRS->CIRS->AA becomes ICRS->CIRS->ICRS->CIRS’->AA so there is double the number of actual transforms, each of which is slower.

There is something about the speed here that I don't understand - if you make a large number of times and ICRS coordinates and broadcast them to make an altaz frame then the second line of code below is about 18x faster than the first...

coos.transform_to(aa_frame)
coos.transform_to(CIRS()).transform_to(aa_frame)

edit: found the problem, and the cause of some of the worst slowdown. Suppose obstime is a (N, 1) array of times and we have M targets. The first line of code first makes a CIRS frame with shape (N, M) and converts that to topocentric CIRS and thus AA. The CIRS->CIRS' transform takes ages due to the large size of the frame.

The second line of code only does the CIRS->CIRS' transform for an input frame of shape (M,) so it's much quicker.

Not sure how to address this, since the behaviour above is baked into the coordinate transform machinery, which I don't really understand.

The way it is coded is slightly slower but the bulk of the slowdown comes from the fact that CIRS->AltAz contains an ICRS self transform to get topocentric CIRS.

Is this true in the case of a solar system body where one already has a topocentric GCRS position? Or is this just in the case where one is asking for the AltAz coordinates of an ICRS SkyCoord where the time and location are necessarily derived from the AltAz frame itself. In that case, couldn't you just derive the CIRS obstime and location from the AltAz frame before leaving the ICRS in the first place? You must already be doing that forobstime as an ICRS SkyCoord has no obstime attribute.

@mkbrewer - see my comment above yours. The majority of the slowdown is not in the transforms themselves, but in the astropy coordinate transform machinery...

For a single point in space and a single obstime, you are correct that the slowdown is due to the call to get_gcrs_posvel

Yeah, I seem to be consistently running one comment behind you. That (N, M) business is way out of my league, but it looks like you're on it.

Well I know what the problem is, but not how to fix it! You are right about the cause of the slowdown for a single point though

@StuartLittlefair - your diagnosis for (N, M) makes sense. Apparently, we broadcast before transforming. That should not be necessary and can probably be solved, but needs some thought... (broadcasting is super powerful but tricky to get exactly right)

Could you raise a separate issue? It does seem quite important generally since presumably this means that single-ICRS-position to AltAz for many times -- surely one of the more common transforms! -- is also much slower than it needs to be.

Thanks @mhvk - I've raised the issue at #10997 and will not try to address it as part of this PR.

I now have a working branch of the direct ICRS->AltAz transform, but it requires many of the same changes to the tests included in this PR, and many of the tests are failing since ICRS->AltAz now gives a different result to ICRS->CIRS->AltAz (because of the innacuracies this PR should fix).

Comments on including that transform as part of this PR? I'd prefer to wait until this was merged and re-base that branch upon this one.

@StuartLittlefair - I'm currently reviewing this PR. It's looking good so my tendency would be to get it in and then do the other one right after - it doesn't matter too much if master is slow for a week or so!

Is the method of specifying the observer's location still up for discussion here? If so, after giving this a lot of thought, I truly believe that the cleanest thing to do would be to use an EarthLocation. Not only is this the natural attribute to use in CIRS<->AltAz transformations, but it makes things a bit easier for the user also as it doesn't require a call to EarthLocation.get_gcrs_posvel() with its unnecessary attendant overhead. It also follows the way SOFA does things more closely. In this scenario, the ErfaAstrom method apci() would be replaced with apco() like this:

    @staticmethod
    def apco(frame_or_coord):
        '''
        Prepare astrometry context used in conversions CIRS <-> ICRS

        Arguments
        ---------
        frame_or_coord: ``astropy.coordinates.BaseCoordinateFrame`` or ``astropy.coordinates.SkyCoord``
            Frame or coordinate instance in the corresponding frame
            for which to calculate the calculate the astrom values.
            For this function, a CIRS frame is expected.
        '''
        jd1_tt, jd2_tt = get_jd12(altaz_frame.obstime, 'tt')
        xp, yp = get_polar_motion(altaz_frame.obstime)
        sp = erfa.sp00(jd1_tt, jd2_tt)
        x, y, s = get_cip(jd1_tt, jd2_tt)
        era = erfa.era00(*get_jd12(altaz_frame.obstime, 'ut1'))
        earth_pv, earth_heliocentric = prepare_earth_position_vel(altaz_frame.obstime)

        astrom = erfa.apco(
            jd1_tt, jd2_tt, earth_pv, earth_heliocentric, x, y, s, era,
            frame_or_coord.location.lon.to_value(u.radian),
            frame_or_coord.location.lat.to_value(u.radian),
            frame_or_coord.location.height.to_value(u.m),
            xp, yp, sp, 0.0, 0.0
        )
        return astrom

This makes the intermediate rotation transform gcrs_to_cirs() very slightly more complicated in that it would be necessary to transform the ITRS coordinates of the CIRS EarthLocation into the GCRS by calling erfa.pvtob() followed by multiplication by the transpose of the celestial to intermediate matrix which is required for this transform anyway. Frankly, I don't understand the intermediate rotation transform cirs_to_gcrs() as it seems to violate the principle of transforms always taking their attributes from the target frame, but the same could be done there as well.

The same could actually apply to TETE as well. Since this is derived from GCRS, the intermediate rotation transforms supply the celestial to intermediate matrix necessary to convert an EarthLocation observer's position attribute there also.

There are a few refinements that could be made to this. The calls to pvtob() which include a geodetic to geocentric transformation could be eliminated (including the one in apco()) by inclusion of get_cirs_posvel() (and get_tete_posvel()) methods in the EarthLocation class. ErfaAstrom.apco() could be used in a direct ICRS to AltAz transform with the inclusion of arguments for the refraction constants. Of course, the refraction constants could be added in to the astrom data structure after the call to ErfaAstrom.apco() also.

It's certainly up for discussion @mkbrewer, but this is probably not the best place, since it would have to be a separate PR to make these changes. The discussion at #10994 is probably the best place to talk about it, since the discussion is heading that way.

I've linked to your comment on that issue.

p.s. Sounds like you did the same thing as I did, but it is very easy to just get one thing out of order and then everything breaks...

Indeed - thanks to all. Especially @mkbrewer without which this wouldn’t have happened.

For lack of a better place: 4.2 is not yet out, and it is the release in which TETE is introduced - should be have a "bug fix" to make it use location instead of obsgeoloc/vel? (We won't be loved for it...)

p.s. I'm trying to create a quick PR so we have something to discuss

I think we should ask politely. There’s a strong case for it, and I can turn it around tomorrow but if it causes too much pain maybe not?

Who handles the 4.2 release? Is it @pllim ?

Who handles the 4.2 release?

Maybe @astrofrog

See #11061 for a PR that does it, if only to show that it does not have to cause delay...

@bsipocz is handling 4.2

Delay is caused by other factors. I can get to the release over the weekend but not much before, so if the bugfix is merged by then it's all good from my side.

Does this mean that the CIRS fix may make it into the 4.2 release also?

Thanks for this!! This saved us from a recent rabbit-hole chasing 100 mas errors on Solar System bodies 🎉 😅

... transformations into the AltAz frames did not properly handle aberration ...

I would like to understand the fix better. Is it really (diurnal) aberration? I would have thought that aberration would be independent of the object distance. It seems more like parallax, with larger effects for closer objects, but it's also smaller than the full geocentric vs topocentric parallax.

I wonder if it has something to do with light-travel time. Take the case of the Moon. It takes up to 20 ms for light to travel the distance from the surface of the Earth to its centre. In that time the Moon moves up to 600 metres in BCRS coordinates (assuming a speed of 30 km / s around the barycentre), which translates to a maximum direction offset of 300 milli-arcseconds. This matches the order of magnitude errors we see.

Oops, the geocentric Moon velocity is more like 1 km / s, so this is more like a 10 mas effect. Still a bit stumped :-)

@ludwigschwardt — The angle produced by aberration is indeed completely independent of object distance. Light does not know how far it’s traveled, and has the same relative velocity to the observer as it arrives whether it came from nearby or from very far away.

But applying the same angle to two vectors of different length will require two different offsets: if v₂ has twice the length of v₁, then the correction applied will need to have twice the length to produce the same angle.

Does that perhaps make better sense of what the code is doing? Or does it not quite match the line you’re looking at?

@ludwigschwardt The problem was that the component of diurnal aberration produced by diurnal parallax was missing. Aberration is dependent on the dot product of the position vector and the velocity vector and the position vector has a diurnal component for nearby objects. Aberration was being computed at the geocenter and, while a correction for diurnal velocity was added, no correction was made for diurnal parallax.

For more on this, please see my comment here: #10356 (comment)

@brandon-rhodes @mkbrewer Thanks a lot for the feedback, it's much appreciated!

I carefully reread your posts several times until something dawned. Let's see if I got it this time...

Let aberration be a function a(v, d) of the velocity vector v of the observer and the direction vector d towards the object (just before applying aberration), which typically involves the dot product of v and d. Then the total aberration for the geocenter (= annual aberration) is

a( v_E, d_E )

and the total aberration for the topocentric observer (= annual + diurnal aberration) is

a( v_T, d_T ) = a( v_E + Δv_T, d_T ).

If I understand you correctly, Astropy instead used to calculate

a( v_E + Δv_T, d_E ).

On top of this we also have the diurnal parallax p = d_T - d_E which is hopefully already taken care of (up to 1 degree for the Moon!).

If everything is linear and you wanted to split up the total topocentric aberration into components, one could then write it as

a( v_T, d_T ) = a( v_E, d_E ) + a( Δv_T, d_E ) + a( v_E, p ) + a( Δv_T, p )
= annual + diurnal + diurnal due to parallax + tiny

So you're saying they left out the last two terms.

For the Moon the maximum error in aberration is on the order of (max parallax) * (annual aberration) = 1/60 * 20" ~ 0.3", which is what I got before, but along a way different route... 😂 At least I now realise that the Moon's own velocity does not affect aberration.

As an aside, it's also interesting that the Moon is at the breakpoint where this matters, because the second and third terms are coincidentally nearly equal in maximum size. That's because |Δv_T| ~ 1/60 |v_E| (0.4 vs 30 km/s), and |p| ~ 1/60 |d_E| for 1 degree parallax. So both terms are <0.3" 😄

That's correct. Aberration was being calculated at the geocenter (your first term), then the position was moved to the EarthLocaton adding diurnal parallax, then the aberration correction for diurnal velocity (your second term) was applied, but your third and fourth terms were, unfortunately, neglected.

In the current master, aberration is calculated at the EarthLocation as it should have been done originally.