Converting coordinates from TEME frame to ECEF frame

[hongzhaozheng]

Hi there,

I am looking for a way to convert the coordinates in TEME frame to ECEF frame. I got the code which utilized both sgp4.api and skyfield.sgp4lib to do the conversion, but the converted coordinates are not accurate as I have the true ECEF coordinates for comparison. The code is provided below. I wonder if this piece of code should be modified since it was published online in 2018.

---

from sgp4.api import Satrec
from sgp4.api import jday
import pandas as pd
import csv
import skyfield.sgp4lib as sgp4lib
from astropy import coordinates as coord, units as u
from astropy.time import Time 
import datetime
import np

TLE_file = "C:/Users/bitry/Desktop/Research/code/u-center/data/2021.12.08/gps_2021.12.08.txt"
TLEs = open(TLE_file, 'r')

L_Name = []
L_1 = []
L_2 = []

i = 1

for line in TLEs:
    j = i
    if i == 1:
        L_Name.append(line)
        j = 2
    elif i == 2:
        L_1.append(line[:69])
        j = 3
    elif i == 3:
        L_2.append(line[:69])
        j = 1
    i = j
    

dataframe = pd.DataFrame(columns = ['Satellite_name', 'Line_1', 'Line_2', 'Position_vector', 'Speed_vector']) 
dataframe.Satellite_name = L_Name
dataframe.Line_1 = L_1
dataframe.Line_2 = L_2
jd, fr = jday(2021, 12, 8, 20, 45, 20.008)  # the starting epoch given in .obs file
time = jd + fr
L_PosVector = []
L_SpeedVector = []

for i in range(len(dataframe)):
    s = dataframe.Line_1[i]
    t = dataframe.Line_2[i]
    satellite = Satrec.twoline2rv(s, t)
    e, r, v = satellite.sgp4(jd, fr)
    r,v = sgp4lib.TEME_to_ITRF(time,np.asarray(r),np.asarray(v)*86400)
    v = v/86400
    L_PosVector.append(r)
    L_SpeedVector.append(v)

dataframe.Position_vector = L_PosVector
dataframe.Speed_vector = L_SpeedVector

---

Many thanks,
Hongzhao

[brandon-rhodes]

The Skyfield conversion from TEME to ECEF is tested against an example calculation from a recent update to the SGP4 algorithm. See the comment above this test for a link; the test itself includes the numbers themselves:

python-skyfield/skyfield/tests/test_earth_satellites.py

Lines 67 to 99 in aa59e2d

	# Note that the following test is based specifically on Revision 2 of
	# "Revisiting Spacetrack Report #3" AIAA 2006-6753 (earlier versions of
	# the PDF use different numbers):
	#
	# http://ww.celestrak.com/publications/AIAA/2006-6753/AIAA-2006-6753-Rev2.pdf
	def test_appendix_c_conversion_from_TEME_to_ITRF():
	rTEME = array([5094.18016210, 6127.64465950, 6380.34453270])
	vTEME = array([-4.746131487, 0.785818041, 5.531931288])
	vTEME = vTEME * 24.0 * 60.0 * 60.0 # km/s to km/day
	jd_utc = julian_date(2004, 4, 6, 7, 51, 28.386)
	d_ut1 = -0.439961
	jd_ut1 = jd_utc + d_ut1 / 86400.0
	xp = -0.140682 * arcsecond
	yp = 0.333309 * arcsecond
	rITRF, vITRF = TEME_to_ITRF(jd_ut1, rTEME, vTEME, xp, yp)
	epsilon = 5e-8 # Why not 1e-8, which would match all of their digits?
	assert abs(-1033.47938300 - rITRF[0]) < epsilon
	assert abs(+7901.29527540 - rITRF[1]) < epsilon
	assert abs(+6380.35659580 - rITRF[2]) < epsilon
	vITRF_per_second = vITRF / seconds_per_day
	epsilon = 7e-8 # Why not 1e-9, which would match all of their digits?
	assert abs(-3.225636520 - vITRF_per_second[0]) < epsilon
	assert abs(-2.872451450 - vITRF_per_second[1]) < epsilon
	assert abs(+5.531924446 - vITRF_per_second[2]) < epsilon

Could you explain which numbers you would like Skyfield to print in your own case, and which numbers it is printing instead for comparison?

[hongzhaozheng]

Hi brandon-rhodes,

I used downloaded TLE file to compute the satellite positions as I did not have eph data directly outputed from my GNSS unit. Forget about the date I used in my previous post. I got reference satellite positions in ECEF from my supervisor for comparison. The epoch I want to use is week number: 1485 and TOW: 152143. And the SVs of interest are [22, 31, 14, 32, 12, 20, 5, 30]. The TLE file I used is given below. Could you verify what I did is right or wrong?

combined.txt

Kindly,
Hongzhao Zheng

[brandon-rhodes]

1. Could you define "TOW" and "SV"? Alas, I am not familiar yet with those acronyms.
2. Could you maybe choose one satellite, one exact time, and share the positions in ECEF you got from your supervisor? Without them, we cannot see what result your supervisor wanted.

[hongzhaozheng]

Of course!

TOW: time of week.
SV: satellite vehicle

The image below contains the satellite positions in ECEF at '2008-06-23 18:15:43'. This time is obtained by converting week number: 1485 and TOW: 152143 using the GPS time calculator provided by LabSat. The positions are organized orderly with the array [22, 31, 14, 32, 12, 20, 5, 30] for the satellites of interest. For example the first column from the left corresponds to the position of SV 22.

[brandon-rhodes]

Thank you! That gives me enough information to move forward.

A first problem appears to be that your script expects 3-line groups "Name / Line 1 / Line 2", but your combined.txt file does not have names for the satellites, and has blank lines and text comments sometimes too, which your script tries to load as TLE lines. If you will print(L_1), for example, you will see lines from the file that do not start with 1, including comment lines. This is probably why the first satellite position I am getting back is [nan nan nan]?

[hongzhaozheng]

The TLE file I attached is assembled by myself. Those TLE data are old, I requested those data via email. Try the first four set of TLE data instead. I have given sat PRN number to each set of TLE data. Note that some of them contains two versions of TLE data for a specific satelltie, you can try any one of them if you wish, just don't use both.

combined_v2.txt

[brandon-rhodes]

After hand-editing the combined_v2.txt file and editing your script to use a date of jd, fr = jday(2008, 6, 23, 18, 15, 43), since that is the date your advisor's coordinates are for, I do see a result that is different from that of your advisor:

[-26356.92540716   2168.507144    -3310.48192975]

(That's the result of print(r) for the first satellite.) As I don't see anything obviously wrong with your script — you are calling TEME_to_ITRF() as I would expect it to be called — I don't know any explanation for why the coordinates don't match those of your advisor. My suggestion would be to try another astronomy library as a double-check. For example, what ITRF coordinates would Skyfield give for one of those satellites? That would provide a third "vote" that might help resolve whether the problem is in your script, or in the coordinates from your advisor (maybe they are in a coordinate frame we don't expect?).

[hongzhaozheng]

I am glad to hear that what I did is right. And the skyfield library is trustworthy. I noticed that many other astronomy libraries require additional inputs which I don't have. May I ask what the other astronomy library would you recommend me to use? Also I want to confirm that the input for jday() should be the epoch of my interest in UTC format, correct?

[brandon-rhodes]

I'm not sure what to recommend beyond Skyfield, as that's of course the one I'm most familiar with. Yes, jday() should be given a UTC date and time, as the TLE format also (so far as I know) uses UTC.

[dcajacob]

What do you mean the TLE is assembled yourself? That may be the issue. Also, If you are trying to match the output of a GPS receiver by using SGP4 to propagate the GPS satellites, you will probably fail. RXs use a simpler propagator based on fitted coefficients that are transmitted by the satellites in the almanac.