Coordinate rotation issues

[cnebelecky]

I'm new to using Skyfield. I've been looking for a full python astrodynamics library that can support satellite research and have previously used Orekit and PyOrekit for this. In trying to get a feel for things I am coding examples in both skyfield and Orekit (java) and noticing some discrepancies when representing satellite positions computed from TLEs and I wonder if it is due to the underlying EOP data or the algorithms used. In my example, given a LEO satellite in near polar orbit I'm computing RMS errors of about 2 meters when comparing the skyfield (GCRS) values to a J2000/GCRF value from Orekit - this isnt a big deal. However, when I rotate for itrs using the following

from skyfield.framelib import itrs
state_itrs = state.frame_xyz_and_velocity(itrs)

and compare that against Orekit's ITRF position the error is 15 meters which is concerning. I was worried that perhaps the EOP data was not correctly loaded so I went to IERS and downloaded a new finals2000A.all file and load it in using

from skyfield.api import Loader
load = Loader('~/skyfield-data/')
ts = load.timescale()

I checked using multiple ITRF and IERS conventions in Orekit and eliminated that as the cause. I also checked with a different LEO satellite, and different date and get the same ~2 meter difference in J2000 and ~15 meter difference in ITRF.

I'm not sure if this is a data or algorithm issue. Can anyone with more skyfield experience weigh in? Thanks!

[brandon-rhodes]

Those differences are so small compared to the inherent few-kilometer inaccuracy of SGP4, that one would almost be tempted to blame them on quirks of how the propagator might have been tweaked between different versions. But it’s indeed interesting that the magnitude of the differences are consistent between two different satellites! That does seem to suggest some systematic difference in how the two libraries handle coordinate systems. (Though, of course, two or three more satellites would make it even more likely that the difference isn’t simply noise, but a consistent 2m / 15m offset!)

Does the Orekit library include the small effect of polar motion? That makes about a ~3m difference at the Earth's surface, and a proportionally bigger difference at larger distances from the Earth’s center. You can turn it on in Skyfield with the maneuver described here:

https://rhodesmill.org/skyfield/accuracy-efficiency.html

Have you double-checked that both libraries are using the same timescale, presumably UTC for the moments you are specifying? You might take a look at whether the 2m and 15m differences are always in the same direction: are they always that the satellite is ahead in its orbit? Or behind? Or farther from the Earth's center, or closer? Or are they not systematic at all and flip around between different directions? The answer would be a very strong hint about where we would investigate next.

Or, instead of guessing: does Orekit expose the underlying vectors? If you could have it print the TEME vector generated by SGP4, then the GCRS vector, then the ITRS vector, then we could have Skyfield print the same three vectors and identify at exactly which step a difference is creeping in.

[cnebelecky]

First thanks for pointing me for instructions how to turn on polar motion. I thought it was by default. I am turning it on as:

load = Loader('~/skyfield-data/')
ts = load.timescale()
with open('./skyfield-data/finals2000A.all') as f:
    finals_data = iers.parse_x_y_dut1_from_finals_all(f)
    iers.install_polar_motion_table(ts, finals_data)

That certainly made the differences in ITRF values smaller by about half. I ran some additional tests using current TLE data for today from SpaceTrack. I grabbed various orbit regimes. I obtained the following differences when comparing results from Orekit and Skyfield:

NORAD 25544		28884		40730		48127
Orbit    LEO            GEO                    GPS                LEO
-------------------------------------------------------------
TEME:  0.000	     0.0000		0.000		0.000
J2000: 1.098	     3.575		       4.040		1.083
GCRF:  0.842	     3.337	               3.679		0.866
ITRF:   6.134        38.292	             19.380	       5.902
NASA: 5.526       35.482               15.961               5.233
Angle: 0.186        0.187                0.151                 0.181     

Additionally, I also used a TEME to ITRF conversion available on github from NASA (PosVelConvert -- See case TEME2ECF ). I compared that against the ICRS values from Skyfield.

All of these are RMS in meters. The SGP4 values are right on for all satellites so its not going to be an in-front, behind, etc thing. It seems the discrepancy is with the actual coordinate rotations. I think the GCRS values are good. The ICRS values are concerning though. As you move outward the discrepancy increases so I'd like to track down where its coming about. One thing I did was to try and figure out the approximate angular difference between the Skyfield and Orekit values by dividing the RMS error by the instantaneous orbit radius. I'm finding that pretty consistently there is about a .18 arcsecond difference as shown in the table above.

With Orekit I'm not sure I get get into the internals to see the individual EO parameters and rotation matrices that make up the rotations but I can with the Matlab code. Does skyfield give a way to see what EO parameters are being used and what the individual rotation matrices are that comprise the various transformations? Maybe that will allow us to hunt down where the discrepancy originates.

[brandon-rhodes]

I think the GCRS values are good. The ICRS values are concerning though.

I had thought of “GCRS” as the Earth’s standard ICRS reference frame, so I'm not sure what’s meant in your table by an ICRS reference frame for satellite vectors separate from the GCRS. Could you clarify the difference between the two?

Does skyfield give a way to see what EO parameters are being used and what the individual rotation matrices are that comprise the various transformations?

At the moment its source code has to be consulted. It looks like the TEME coordinates are handed to the TEME frame of reference for translation into the GCRS:

python-skyfield/skyfield/sgp4lib.py

Line 192 in aa59e2d

R = _T(TEME.rotation_at(t))

That reference frame is defined further down in the file. It appears to consist entirely of a z-axis rotation combined with M, the big matrix that converts between the GCRS/ICRS reference frame and the Earth’s dynamical reference system of-date:

python-skyfield/skyfield/sgp4lib.py

Line 292 in aa59e2d

def rotation_at(t):

I think those are the only two steps, the z-axis rotation and M, that are made to perform the transform?

[cnebelecky]

Sorry when I said "ICRS" in the text I mean "ITRF". What is shown in the table is "ITRF" which is the Earth fixed frame that I get using

from skyfield.framelib import itrs
state_itrs = state.frame_xyz_and_velocity(itrs)

I will consult the source code and see if I can tease anything out.