Do topocentric frames self-transform properly through ICRS?

[StuartLittlefair]

When adding the TETE frame I mentioned that there is an unusual issue with topocentric frames that perhaps reveals a bug in those frames (TETE, GCRS and the soon-to-be implemented topocentric CIRS).

The issue arises with transforms from ICRS to a topocentric coordinate frame, and I'm not 100% sure yet that it is wrong, but I am raising this issue so I don't forget the detective work I've done to date.

You can most easily see the behaviour if you create an ICRS position for a location just above the Earth's surface:

from astropy.coordinates import *
from astropy.time import Time
from astropy import units as u

t = Time('J2010')
obj = EarthLocation(-1*u.deg, 52*u.deg, height=10.*u.km) 
home = EarthLocation(-1*u.deg, 52*u.deg, height=0.*u.km)
obsloc_gcrs, obsvel_gcrs = home.get_gcrs_posvel(t)

# create geocentric and topocentric GCRS frames
gcrs_geo_frame = GCRS(obstime=t)
gcrs_topo_frame = GCRS(obstime=t, obsgeoloc=obsloc_gcrs, obsgeovel=obsvel_gcrs)

icrs_coo = obj.get_itrs(t).transform_to(ICRS())

Now if you take this ICRS coordinate and transform it back in to the topocentric and geocentric frames, you can see that they are not separated by the observatory position vector as one might expect:

In [333]: crepr1 = icrs_coo.transform_to(gcrs_geo_frame).cartesian - obsloc_gcrs 
     ...: crepr2 = icrs_coo.transform_to(gcrs_topo_frame).cartesian                                                                      

In [334]: (crepr1-crepr2).norm()                                                                                                         
Out[334]: <Quantity 0.63376095 km>

In [335]: np.arccos(crepr1.dot(crepr2)/crepr1.norm()/crepr2.norm()).to(u.deg)                                                            
Out[335]: <Quantity 3.62610794 deg>

I noticed this when I implemented topocentric CIRS, as the same issue caused some of our AltAz tests to break - see this comment.

After a lot of digging I have tied down the root of the issue to the aberration calculation in atciqz. If I comment out the aberration lines then the issue goes away. What I haven't yet decided if the above code is the correct answer. My naiive understanding of aberration is that it should never cause a discrepancy as large as over 3 degrees, and the discrepancy is larger for points closer to the Earth's surface, which doesn't seem to make sense either...

Tagging @mhvk and @mkbrewer in case they have any bright ideas

[mhvk]

Just to get my thinking right: if I simplify your example from the comment a bit (hey, equator is easier), then in current master,

from astropy.time import Time
from astropy.coordinates import EarthLocation, AltAz
from astropy import units as u, constants as const
t = Time('J2010')
obj = EarthLocation(0*u.deg, 0*u.deg, height=1*u.km)
home = EarthLocation(0*u.deg, 0*u.deg, height=0*u.km)
aa = obj.get_itrs(t).transform_to(AltAz(obstime=t, location=home))
aa
# <AltAz Coordinate (obstime=J2010.000, location=(6378.137, 0., 0.) km, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0.0, obswl=1.0 micron): (az, alt, distance) in (deg, deg, km)
#     (90., 89.99991111, 1.)>
(aa.alt-90*u.deg).to(u.arcsec)
# <Angle -0.32000124 arcsec>
((aa.alt-90*u.deg)*const.c).to(u.km/u.s, u.dimensionless_angles())
# <Quantity -0.46510096 km / s>

Hence, the deviation here reflects the rotation of the Earth. This seems correct.

With your new implementation, does the same larger factor occur for AltAz? Do you have a difference with current master somewhere to look at and ponder?

[StuartLittlefair]

Hi @mhvk - thanks for playing. Yep, the current master is ok in this regard. This is because it handles topocentric CIRS by adding in the observatory location manually before transforming to AltAz. This leads to the issues at high precision that @mkbrewer pointed out in #10356.

In my branch with a topocentric CIRS, you get the following:

In [1]: from astropy.time import Time 
   ...: from astropy.coordinates import EarthLocation, AltAz 
   ...: from astropy import units as u, constants as const 
   ...: t = Time('J2010') 
   ...: obj = EarthLocation(0*u.deg, 0*u.deg, height=1*u.km) 
   ...: home = EarthLocation(0*u.deg, 0*u.deg, height=0*u.km) 
   ...: aa = obj.get_itrs(t).transform_to(AltAz(obstime=t, location=home))      

In [2]: aa                                                                      
Out[2]: 
<AltAz Coordinate (obstime=J2010.000, location=(6378.137, 0., 0.) km, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0.0, obswl=1.0 micron): (az, alt, distance) in (deg, deg, km)
    (274.16090083, 57.21235149, 1.18953628)>

In [3]: (aa.alt-90*u.deg).to(u.arcsec)                                          
Out[3]: <Angle -118035.53465035 arcsec>

In [4]: ((aa.alt-90*u.deg)*const.c).to(u.km/u.s, u.dimensionless_angles())      
Out[4]: <Quantity -171556.95975485 km / s>

However, after a lot of digging I realised the reason for the above is a root cause that affects all our topocentric frames, including GCRS *in master. I actually think the snippet I posted in the first post is the clearest way of seeing the issue. Especially if you ignore how I got to icrs_coo and just think of it as any ICRS position that is close to the surface of the Earth.

If you transform that ICRS position into two GCRS frames at the same obstime/orientation, but one geocentric/topocentric respectively, then the difference between them should be the GCRS observatory position, to high accuracy. Diurnal aberration can induce a difference of up to 0.32" as you've seen. However, that's not what we see in astropy master. The "culprit" is the call to erfa.ab in atciqz...

[StuartLittlefair]

Well I now understand what's happening here I believe.

Both the geocentric and topocentric positions have an angular displacement from the natural direction of around 20.5", as expected for aberration. But the parallax between them means that the vector from the observatory to the geocentric coordinate is rather large:

I'm still thinking about what this means the correct value for aa should be in the following snippet...

from astropy.time import Time
from astropy.coordinates import EarthLocation, AltAz
from astropy import units as u, constants as const
t = Time('J2010')
obj = EarthLocation(0*u.deg, 0*u.deg, height=1*u.km)
home = EarthLocation(0*u.deg, 0*u.deg, height=0*u.km)
aa = obj.get_itrs(t).transform_to(AltAz(obstime=t, location=home))

[mkbrewer]

Well, this is because you aren't handling the geocentric to topocentric transformation correctly here. To do that properly, you need to subtract the sum of the position of the Earth wrt the solar system barycenter and the GCRS position of the observer from the ICRS position of the object before transforming to GCRS. You also need the velocity of the observer wrt to the geocenter added to the velocity of the geocenter for the aberration calculation to yield the correct result.

[mkbrewer]

Another very minor issue here is that in transforming the ITRS position of the object to ICRS, polar motion is added in. This is because any ITRS position is assumed to be attached to the Earth. This is the reason that the AltAz transformation yields 89.999 degrees instead of 90.000 for Alt.

[StuartLittlefair]

Well, this is because you aren't handling the geocentric to topocentric transformation correctly here. To do that properly, you need to subtract the sum of the position of the Earth wrt the solar system barycenter and the GCRS position of the observer from the ICRS position of the object before transforming to GCRS. You also need the velocity of the observer wrt to the geocenter added to the velocity of the geocenter for the aberration calculation to yield the correct result.

@mkbrewer Thanks for taking time out on Election Day to reply ;) I don’t understand your answer though, because that’s precisely what we do, isn’t it? In icrs_to_gcrs:

        astrom_eb = CartesianRepresentation(astrom['eb'], unit=u.au,
                                            xyz_axis=-1, copy=False)
        newcart = icrs_coo.cartesian - astrom_eb

        srepr = newcart.represent_as(SphericalRepresentation)
        gcrs_ra, gcrs_dec = atciqz(srepr.without_differentials(), astrom)

        newrep = SphericalRepresentation(lat=u.Quantity(gcrs_dec, u.radian, copy=False),
                                         lon=u.Quantity(gcrs_ra, u.radian, copy=False),
                                         distance=srepr.distance, copy=False)

Since astrom is from erfa.apcs then I think astrom[‘eb’] is vector from the observer to the barycentre.

[mkbrewer]

That is correct, but this is not the example that you presented when you opened this particular issue. There crepr2 does exactly that and therefore calculates the correct topocentric position of an object in the GCRS given its position in the ICRS. However, this is not what crepr1 is doing. As you noticed, this is calculating the difference between the geocentric GCRS position of the object and that of the observer. I'm not sure what relevance this has to calculating the AltAz position of an object directly overhead in the ITRS anyway. The simplest way to do that, is to assume that the given position of the object is in the TIRS, transform that into the CIRS by adding the earth rotation angle, transform the observer's ITRS location into the CIRS by adding in polar motion and the earth rotation angle, take the difference to calculate the topocentric intermediate position and then call atioq() after setting up the proper context as I outlined here: #10356 (comment). If you really want to transform the object's location into the ICRS, you need to be careful how you do that. Aberration is a tricky business since it depends on the observer's point of view (as you noted). To get the proper ICRS position of the object as viewed by the observer here, you need to transform the object's topocentric intermediate position into the BCRS by setting up the proper context using apco() or something equivalent, calling aticq() and finally adding the observer's GCRS position in. In the context of astropy as currently written, this means calculating the BPN matrix, transforming both the object's topocentric (CIRS) intermediate position and the observer's CIRS position into the GCRS and then calling gcrs_to_icrs().

[StuartLittlefair]

@mkbrewer thank you for the comments. I will take time to digest them. What I think you are missing is that astropy does not have, and I hope does not need, bespoke transforms from one frame into every other frame. Instead, we have a chain of transforms that can be linked to get from one frame to any other.

So the way this is done now is to transform the ITRS position into geocentric CIRS, transform that into topocentric CIRS and transform that into AltAz.

It’s all done as you suggest except for the geocentric to topocentric transform. As I’ve mentioned elsewhere this “self transform” must go through ICRS because it is designed to be general - ie to transform between CIRS frames at different times as well as topo/geocentric.

And it is that self transform that I am fighting to understand. I illustrated with GCRS because it shows the same behaviour, ie the GCRS-ICRS-GCRS’ chain is the same as for CIRS and TETE apart from the BPN matrix.

[mkbrewer]

Oh, and you also need to calculate the observer's obsgeovel.

[mkbrewer]

Well, this situation seems quite artificial. Is there any real world application for specifying the position of a source in the ITRS (or more properly in the TIRS)? A far more likely case is given an AltAz, calculate the CIRS position of the source. This naturally yields a topocentric position, so I was taking that as a starting point and, for this test only, suggesting the proper way to get from an observer's location and an object's location in ITRS/TIRS coordinates to the starting topocentric CIRS position. From there, you can use standard transforms to go from one observer's location to another or one obstime to another by transforming the topocentric CIRS position back though topocentric GCRS to ICRS/BCRS. Then changing the observer's position or obstime is just a matter of undoing space motion in the case of stars or doing a new look up in the case of solar system bodies as in the (hopefully ongoing) discussion at issue #10404. The main thing to think about here is reversibility. You want to take the exact same steps going backwards as you do going forwards so that a round trip without changing anything will yield the same set of topocentric CIRS coordinates and hence the same AltAz.

[mkbrewer]

One main point that I'm attempting to convey is that the only proper way to go from topocentric CIRS to geocentric CIRS is through topocentric BCRS (astrometric). One can't separate a topocentric CIRS position into geocentic positions of the source and the observer and then transform the geocentric CIRS position of the source to ICRS/BCRS independent of the observer's position and expect that transform to yield the proper ICRS/BCRS position of the source.

[StuartLittlefair]

Thanks @mkbrewer for your comments. Having slept on them I think I can clarify what my problem is, and the relationship to the AltAz problem. Apologies, this is a long post. It serves to clarify my thinking as much as anything else!

You wrote that the simplest way to get a topocentric CIRS is:

assume that the given position of the object is in the TIRS, transform that into the CIRS by adding the earth rotation angle, transform the observer's ITRS location into the CIRS by adding in polar motion and the earth rotation angle, take the difference to calculate the topocentric intermediate position

I completely agree this is the simplest way. So below I do that within astropy:

from astropy.coordinates import *
from astropy.time import Time
from astropy import units as u
from astropy.coordinates.builtin_frames.intermediate_rotation_transforms import cirs_to_itrs_mat, gcrs_to_cirs_mat
from astropy.coordinates.matrix_utilities import matrix_transpose

t = Time('J2010')
obj = EarthLocation(-1*u.deg, 52*u.deg, height=10.*u.km) 
home = EarthLocation(-1*u.deg, 52*u.deg, height=0.*u.km)
obsloc_gcrs, obsvel_gcrs = home.get_gcrs_posvel(t)

itrs_coo = obj.get_itrs(t)
pmat = cirs_to_itrs_mat(t)  # earth rotation angle and polar motion

# Geocentric CIRS position of object
cartrepr = itrs_coo.cartesian.transform(matrix_transpose(pmat))

# CIRS position of observatory
obsrepr = home.get_itrs(t).cartesian.transform(matrix_transpose(pmat))

# Topocentric CIRS position of object (MKB's method)
obj_cirs_repr1 = cartrepr-obsrepr

The other way to do this is harder, but we do it this way in astropy so that CIRS->CIRS' will work for changes in obstime as well as topocentric->geocentric changes. It is to create a geocentric CIRS coordinate for this object and transform, through ICRS, into topocentric CIRS. Below is that method:

cirs_geo = CIRS(cartrepr, obstime=t)
topocentric_cirs_frame = CIRS(obstime=t, obsgeoloc=obsloc_gcrs, obsgeovel=obsvel_gcrs)
cirs_topo = cirs_geo.transform_to(topocentric_cirs_frame)

# to clarify that this is the same as CIRS->ICRS->CIRS'
cirs_topo_alt = cirs_geo.transform_to(ICRS()).transform_to(topocentric_cirs_frame)

# cartesian representation for comparison with method above
obj_cirs_repr2 = cirs_topo.cartesian
obj_cirs_repr2_alt = cirs_topo_alt.cartesian
print(obj_cirs_repr2_alt - obj_cirs_repr2)
# (0., 0., 0.) km

But now I have demonstrated my problem, because look at the difference between the two methods!

In [21]: obj_cirs_repr2 - obj_cirs_repr1                                                           
Out[21]: 
<CartesianRepresentation (x, y, z) in km
    (0.6361604, 0.08214658, 0.01567519)>

In [22]: (obj_cirs_repr2 - obj_cirs_repr1).norm()                                                  
Out[22]: <Quantity 0.64163372 km>

For some reason, getting from geocentric to topocentric CIRS through ICRS is not working as expected. The difference between the two is in some way caused by the handling of aberration, since if I just comment out aberration from the transforms the difference vanishes. That's what my crappy diagram is supposed to illustrate too.

This ties everything together; because my CIRS->AltAz fix is "failing" astropy tests that expect an altitude of 90 degrees for obj observed from home. The root cause is that ITRS->AA goes through ITRS->CIRS(geo)->ICRS->CIRS(topo)->AltAz. But this gives a different topocentric CIRS to the one we expect, and so a different AltAz:

altaz_frame = AltAz(obstime=t, location=home)
topocentric_cirs_frame.realize_frame(obj_cirs_repr1).transform_to(altaz_frame)
# <AltAz Coordinate (obstime=J2010.000, location=(3934.36115314, -68.67452938, 5002.80334548) km, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0.0, obswl=1.0 micron): (az, alt, distance) in (deg, deg, km)
#    (303.8534313, 90., 10.)>
cirs_topo.transform_to(altaz_frame)
# <AltAz Coordinate (obstime=J2010.000, location=(3934.36115314, -68.67452938, 5002.80334548) km, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0.0, obswl=1.0 micron): (az, alt, distance) in (deg, deg, km)
#    (272.26090326, 86.32878441, 10.02067844)>

Sure - I can special case the ITRS->AltAz transform, but it won't fix the issue that I've now lost confidence in the CIRS "self-transform". That issue still needs addressing and might indicate a bug either in our CIRS->ICRS or our ICRS->CIRS transform. And yes, this is an artificial example, but it is one of our test cases designed to stop bugs from recurring, so it needs to pass.

The reason why I demonstrated this issue in GCRS originally is because it affects all our topocentric frames (GCRS, CIRS, TETE) and seeing it in GCRS re-assures me it's not something I've got wrong in CIRS specifically. Here's the exact same issue, illustrated using GCRS (nb: this is just a re-ordering of my initial example):

# geocentric GCRS
pmat = gcrs_to_cirs_mat(t)
obj_gcrs_geo = cartrepr.transform(matrix_transpose(pmat))  # cartrepr is CIRS
# check against astropy normal way
print(itrs_coo.transform_to(GCRS(obstime=t)).cartesian - obj_gcrs_geo)
# (0., 0., 0.) km

# topocentric GCRS (method 1)
obj_gcrs_repr1 = obj_gcrs_geo - obsloc_gcrs
# method 2
gcrs_topo_frame = GCRS(obstime=t, obsgeoloc=obsloc_gcrs, obsgeovel=obsvel_gcrs)
obj_gcrs_repr2 = GCRS(obj_gcrs_geo, obstime=t).transform_to(gcrs_topo_frame).cartesian
print(obj_gcrs_repr2 - obj_gcrs_repr1)  
# (0.63617581, 0.08214677, 0.01503589) km
print((obj_gcrs_repr2 - obj_gcrs_repr1).norm()) 
# 0.641633716491988 km

which is exactly the same discrepancy seen in CIRS.

So finally perhaps I have managed to explain the concerns I have about the topocentric frames in astropy, which is holding up my fix for #10356? I have two ways of getting a topocentric position which should agree, but don't.

edit: if you see the discussion below, I believe the reason they don't agree is that the two methods disagree over whether obj and home's ITRS coordinates include aberration of a light source or not.

[StuartLittlefair]

My main question is, which calculation above is correct? What I think is that the coordinates calculated "the astropy way" are correct. That is, for a setup like the one below, the altitude of an object obj from home should not be 90 degrees:

obj = EarthLocation(-1*u.deg, 52*u.deg, height=10.*u.km).get_itrs(t)
home = EarthLocation(-1*u.deg, 52*u.deg, height=0.*u.km).get_itrs(t)

My thinking is thus. Since these coordinates are specified in ITRS, which is geocentric and "downstream" of GCRS they include the effects of aberration. What this code is actually saying is that - for a geocentric observer - the source appears to be overhead the observatory at home†. Due to aberration, the proper direction of the source is different by about 20.5".

However, an observer at home, observing a source at the proper position of the source, will disagree strongly about the natural direction of the source, as shown below:

In this case, trig implies a difference in natural direction of about 3.6 degrees, suggesting the observer at home sees the source at an altitude of 86.3 degrees, which is what my CIRS branch predicts.

In short, I believe the tests in astropy, which expect an altitude of 90 degrees, are wrong. My CIRS branch is right. If I get some agreement, I will change the tests and push my CIRS branch for comments.

†the confusion arises because the code encourages the user to think that we are saying obj appears to be directly above home, but in fact that is not the case!

[mkbrewer]

The main issue here is that your topocentric GCRS method 1 is wrong. One can't form a topocentric GCRS position after aberration is applied. Method 2 is correct as it forms the topocentric position first and then applies aberration. The aberration correction depends on both the relative velocities and relative positions of the observer and the source through the expression P dot V. This is the heart of the problem with astropy's current implementation of the CIRS frame. It completely ignores the contribution to the aberration correction from diurnal parallax as I discussed at length here: #10356. Since the topocentric CIRS position is derived from the topocentric GCRS position by the application of the BPN matrix, in order to go backwards to ICRS, one must undo the topocentric aberration before one can transform from topocentric to geocentic.

Regarding what is meant by the ITRS position of the source, both topocentric and geocentric points of view are correct, but mutually exclusive. One must decide which of these is relevant to the case at hand. For your AltAz test, the ITRS position of the source is clearly meant to be from the point of view of the observer (i.e. topocentric) The source should be seen directly overhead after the aberration correction is made. In order for this to happen, the aberration correction must include the topocentric position of the source.

[mkbrewer]

If you transform obj_gcrs_repr2 to CIRS by applying the BPN matrix and then to AltAz, I think that you will find that the source is, in fact, at Alt = 90.

[mkbrewer]

Oh darn. I got that backwards. Method 1, is, in fact correct. It's method 2 that is wrong. Since we are coming from CIRS, method 1 is equivalent to transforming obj_cirs_repr1 to GCRS, whereas method 2 will undo geocentric aberration and then add in topocentric aberration in transform_to(gcrs_topo_frame) leading to the wrong result. Sorry about that.

[StuartLittlefair]

Actually, I believe you were right the first time. If I specify any position in ITRS, CIRS or GCRS it must be the natural direction of the source, including aberration. To get the vector to that source from a topocentric position you must undo the aberration, find the proper direction, correct for parallax and then re-do the aberration, which is what method 2 does.

This is what the figure show..

[mkbrewer]

The following transforms should be correct:

obj_icrs = GCRS(obj_gcrs_repr1, obstime=t, obsgeoloc=obsloc_gcrs, obsgeovel=obsvel_gcrs).transform_to(ICRS).cartesian
obj_gcrs_repr2 = obj_icrs.transform_to(gcrs_topo_frame).cartesian

The first transform will undo topocentric aberration to yield the proper ICRS position and then the second one will add it back in for a consistent round trip. Please pardon any syntax errors that I may have made there.

[mkbrewer]

I think the way to look at this is that in forming obj_gcrs_repr1, you haven't done anything in regards to aberration. You've just taken a topocentic CIRS position of an object directly overhead from the point of view of the observer in the ITRS and transformed that to GCRS. Then the transform of that topocentric GCRS position to ICRS will correctly undo the topocentric aberration.

[mkbrewer]

Another way to look at it is that if you transform obj_cirs_repr1 = cartrepr-obsrepr to GCRS, that is equivalent to transforming cartrepr and obsrepr to GCRS separately and then taking the difference which is essentially what you have done in forming obj_gcrs_repr1 = obj_gcrs_geo - obsloc_gcrs.

[mkbrewer]

My initial confusion came from thinking of obj_gcrs_geo as the position of the object as viewed from the geocenter, which it is not. In fact, it's the geocentric position of the the object as viewed by the topocentric observer.

[mkbrewer]

If you try it, I think that the transformation above #10983 (comment) will yield Alt=90 whereas this:

gcrs_geo_frame = GCRS(obstime=t)
obj_gcrs_repr2 = obj_icrs.transform_to(gcrs_geo_frame).cartesian

will yield the offset shown in your figures.

[StuartLittlefair]

Yes, I agree and the transformation in #10983 (comment) does give Alt=90 because it's the right way to construct a topocentric GCRS position for an object that appears overhead to a topocentric observer (which would not appear overhead for a geocentric observer)

[mhvk]

I've slightly lost track although I think I get the gist - is the existing test that fails basically set up improperly?

In any case, do try to add something to the documentation about these pitfalls as part of the PR! My contribution to reviewing can then be whether that helps me understand it...