GCRS to ITRS conversion results differ from Astropy

[shfxia]

I compare the transform results from GCRS to ITRS between astropy and skyfield. There is a considerable difference in the results.

import astropy.units as u
from astropy.coordinates import GCRS, ITRS, CartesianRepresentation
from astropy.time import Time

from skyfield.api import load
from skyfield.positionlib import Geocentric
from skyfield.constants import AU_M
from skyfield.framelib import itrs

import numpy as np

# GCRS
rc = np.array([2395723.894320751, 5701333.894320751, 2337419.2469580634])

# astropy GCRS to ITRS
obstime = Time("J2000.0")
rc_ecef = (
    GCRS(CartesianRepresentation(x=rc, unit=u.m), obstime=obstime)
    .transform_to(ITRS(obstime=obstime))
    .cartesian.xyz.to_value(u.m)
)

# skyfield GCRS to ITRS
t = load.timescale().J2000
p_gcrs = Geocentric(rc / AU_M, t=t)
p_itrs = p_gcrs.frame_xyz(itrs)

print("astropy", rc_ecef)
print("itrs   ", p_itrs.m)

Results,

astropy [-5187390.74533245  3367000.53859377  2337202.2840885 ]
itrs    [-5187391.23614591  3367004.82139493  2337195.02486631]

I'd like to know if I'm using the feature incorrectly, or if there are any other differences in details that need to be considered.

[brandon-rhodes]

Those results look very close — try two things:

1. First, do a bit of trigonometry to figure out what angular difference those few meters of difference make at that range. Are we looking at a big difference (multiple arcseconds or even arcminutes), or are we below the arcsecond level?
2. If the latter, maybe Skyfield and Astropy have loaded up different estimates for ∆T. See:

https://rhodesmill.org/skyfield/time.html#ut1-and-downloading-iers-data

[shfxia]

@brandon-rhodes Thanks for your advice. Finally, I found the few meters of difference is caused by polar motion. Because I didn't load polar motion data in skyfield. After I made the following changes,

with load.open("tests\\finals2000A.all") as f:
    finals_data = iers.parse_x_y_dut1_from_finals_all(f)
ts = load.timescale()
iers.install_polar_motion_table(ts, finals_data)

the results difference is quite small.

astropy : [-5187390.74533245  3367000.53859377  2337202.2840885 ]
skyfield: [-5187390.74414268  3367000.5410562   2337202.28318177]
astropy - skyfield: [9.04944e-05 "  3.02421e-05 "  3.72529e-09 m]

By, the way without load polar motion data, the spherical difference is :

astropy - skyfield: [0.110905 "  0.242116 "  3.72529e-09 m]

[brandon-rhodes]

Excellent, I'm glad you found the source of the difference you were seeing.

[shfxia]

@brandon-rhodes I found ITRF_to_GCRS2 in positionlib.py and ITRF_to_GCRS2_fast in sgp4_parallel.py didn't apply polar motion too. While in ITRF_to_GCRS in position.lib the polar motion is applied.

    def rotation_at(t):
        R = mxm(rot_z(-t.gast * tau / 24.0), t.M)
        if t.ts.polar_motion_table is not None:
            R = mxm(t.polar_motion_matrix(), R)
        return R

[brandon-rhodes]

Correct; several old deprecated routines indeed did not account for polar motion, which is why they are no longer used, and are not mentioned in Skyfield's documentation.