ITRS/GCRF transform validation against Spice & SOFA

Astropy
1

Hello Astropy community.

I’m struggling to get a ITRF/GCRF conversion from Astropy, Spice and Sofa that are in agreement with each other.

Basically, Spice and Sofa are within 4 arcsecs, using the following implementation for Sofa

However, Astropy and Sofa differ by nearly 12 arcsec.

The transform evaluation time is defined as a TAI string

time_string_TAI = "2020-12-31T00:00:00" # TAI timestring

The ITRF_GCRF DCM is computed like so:

def gcrs_to_itrs(t):
    return SkyCoord(x=np.array([1,0,0]),
                    y=np.array([0,1,0]),
                    z=np.array([0,0,1]),
                    representation_type='cartesian'
                   ).transform_to(ITRS(obstime=t)).cartesian.get_xyz().value

...

INPUT_TIME = Time(time_string_TAI,scale = "tai")
ITRF_GCRF_astropy = gcrs_to_itrs(INPUT_TIME)

The UT1-UTC offset provided to SOFA is the same as computed by Astropy at the evaluation time.

The EOP parameters are left to zero when calling SOFA, which contributes to the overall error (I’m guessing that Astropy retrieves the EOPs corresponding to the evaluation tome internally) but to a much smaller degree than the few arcsecs of agreement I’m expected to find.

The two DCMs returned from Astropy and SOFA respectively are

ITRF_GCRF_astropy : [[-1.64352104e-01  9.86420571e-01  3.30480317e-04]
 [-9.86399708e-01 -1.64239026e-01  2.08095097e-03]
 [ 2.00068122e-03 -6.22881226e-06  9.99997780e-01]]

ITRF_GCRF_sofa : [[-1.64336961e-01  9.86404205e-01  3.28918388e-04]
 [-9.86402214e-01 -1.64337290e-01  1.98161946e-03]
 [ 2.00873133e-03  1.20749412e-06  9.99997982e-01]]

The error between the two, expressed as a principal rotation vector is

prv_error : [-4.50870748e-05  3.77193352e-06  4.21492558e-05]

The Z component is along the angular velocity vector. Because this component is not exceedingly larger than the other two, a time offset can be ruled out (otherwise the orientation error would predominantly be around the angular velocity vector).

I would be grateful is someone could point me in the right direction regarding this inconsistency. I don’t think Astropy expects additional transform parameters to correctly instantiate the ITRF_GCRF DCM ?

Thanks !

2

Found the culprit: by default, SkyCoord assumes that it is tied to the ICRS. Adding frame = gcrs to the SkyCoords constructor did the trick

def gcrs_to_itrs(t):
    return SkyCoord(x=np.array([1,0,0]),
                    y=np.array([0,1,0]),
                    z=np.array([0,0,1]),
                    representation_type='cartesian',
                    frame = 'gcrs'
                   ).transform_to(ITRS(obstime=t)).cartesian.get_xyz().value
3

i’m not quite there yet.

i’ve displayed the error in the ITRF/GCRF transform over the course of several months : it features sinusoidal patterns repeating every year

image
image1839×975 126 KB

The transform is still extracted from Astropy using

def gcrs_to_itrs(t):
    return SkyCoord(x=np.array([1,0,0]),
                    y=np.array([0,1,0]),
                    z=np.array([0,0,1]),
                    representation_type='cartesian',
                    frame = 'gcrs'
                   ).transform_to(ITRS(obstime=t)).cartesian.get_xyz().value

The sofa transform is obtained from the following :

if (iauTaiutc ( tai1, tai2 ,&utc1, &utc2 ) !=0){
        return 1;
    };

    /* From UTC get TT and UT1. */
    if (iauTaitt ( tai1, tai2, &tt1, &tt2 ) != 0){
        return 1;
    }

    if (iauUtcut1 ( utc1, utc2, dut1, &ut11, &ut12 ) != 0){
        return 1;
    };

    /* ========================= */
    /* IAU 2006/2000A, CIO based */
    /* ========================= */

    /* CIP and CIO, IAU 2006/2000A. */
    iauXys06a ( tt1, tt2, &x, &y, &s );

    /* Add CIP corrections. */
    x += dx06;
    y += dy06;

    /* GCRS to CIRS matrix. */
    iauC2ixys ( x, y, s, rc2i );

    /* Earth rotation angle. */
    era = iauEra00 ( ut11, ut12 );

    /* Form celestial-terrestrial matrix (no polar motion yet). */
    iauCr ( rc2i, rc2ti );
    iauRz ( era, rc2ti );

    /* Polar motion matrix (TIRS->ITRS, IERS 2003). */
    sp = iauSp00 ( tt1, tt2 );
    iauPom00 ( xp, yp, sp, rpom );

    /* Form celestial-terrestrial matrix (including polar motion). */
    iauRxr ( rpom, rc2ti, rc2it );

    for (size_t i = 0; i < 3; ++i){
        for (size_t j = 0; j < 3; ++j){
            ITRF_GCRF[i * 3 + j] = rc2it[i][j];
        }
    }

I’m obviously missing something …

4

SOFA’s official documentation provides a validation case on page 25/26 of this document. The time at which the iTRF/GCRF transform is evaluated is

2007 April 05, 12h 00m 00s.0 UTC

The document reports the following ITRF/GCRF transform at the prescribed time

9.73104317697536e-1,   2.30363826239128e-1, -7.03163481769e-4
−2.30363800456036e-1 , 9.73104570632801e-1,  1.18545368117e-4
7.11560162594e-4,      4.6626402444e-5 ,     9.99999745754024e-1

I’m getting the following from SOFA

9.73104318e-01,  2.30363824e-01, -7.03165729e-04
-2.30363798e-01, 9.73104571e-01,  1.18543985e-04
7.11562031e-04,  4.66282649e-05,  9.99999746e-01

whereas Astropy is returning

9.73087832e-01,   2.30484277e-01, -5.99206719e-04
-2.30433522e-01,  9.73076049e-01,  2.04479220e-05
6.80499271e-04,   1.57197846e-05 , 9.99999820e-01

Finally, Spice is giving the following transform

9.73104319e-01   2.30363819e-01 -7.03177223e-04
-2.30363793e-01  9.73104572e-01  1.18553063e-04
 7.11575308e-04  4.66220747e-05  9.99999746e-01

Astropy appears to deviate from both SOFA, Spice and SOFA’s reference document for a reason I can’t nail down, can anyone shed some light on this ?