Calculated Classical Orbital Elements of the ISS seem to differ from the actual ones

Question

I am trying to calculate the Classical Orbital Elements of the ISS using Python, and then compare them with the "actual" ones provided on the Internet. However, my obtained data differs by a lot from the "actual" one.

I am using the ISS TLE from this website, which at the time I am writing this (2020-12-04 04:19) it says:

1 25544U 98067A   20339.83283773  .00016717  00000-0  10270-3 0  9038
2 25544  51.6418 224.1133 0001925 115.0430 245.0920 15.49154811 18464
Epoch (UTC):    04 December 2020 19:59:17

Now, using a modified version of the code from the answer from this question: How do I obtain geocentric position vectors at three successive times before I use Gibbs' method?,

from skyfield.api import EarthSatellite, Topos, load
import numpy as np
ts = load.timescale()
minutes = np.linspace(0, 60, 3)
t = ts.utc(2020, 12, 4, 20, minutes, 0)
l1 = '1 25544U 98067A   20339.83283773  .00016717  00000-0  10270-3 0  9038'
l2 = '2 25544  51.6418 224.1133 0001925 115.0430 245.0920 15.49154811 18464'
satellite = EarthSatellite(l1, l2, name='ISS (ZARYA)')
geocentric = satellite.at(t)
r1, r2, r3 = geocentric.position.km
print(f"x/y/z of r1: {r1}")
print(f"x/y/z of r2: {r2}")
print(f"x/y/z of r2: {r3}")

I was able to obtain these positional vectors:

x/y/z of r1: [-4755.01172716  4947.37677595   377.21687117] # at '2020-12-04 20:00:00 UTC'
x/y/z of r2: [-4851.24226498  -367.40786871  5177.39228056] # at '2020-12-04 20:30:00 UTC'
x/y/z of r2: [  268.17962108  4636.21164042 -4395.73178942] # at '2020-12-04 21:00:00 UTC'

After getting these, I used Gibbs' method to obtain the following v2

v2 = [3.0406; -6.1187; 2.3591]

At this point, I use poliastro's rv2coe function to find the COEs:

from poliastro.constants import GM_earth
from poliastro.core.elements import rv2coe
from astropy import units as u
import numpy as np
k = GM_earth.to(u.km ** 3 / u.s ** 2).value
r2 from the obtained positional vectors
r = np.array([-4851., -367., 5177.])
v2 calculated from Gibbs' method
v = np.array([3.040, -6.118, 2.359])
p, ecc, inc, raan, argp, nu = rv2coe(k, r, v)
print("p:", p, "[km]")
print("ecc:", ecc)
print("inc:", np.rad2deg(inc), "[deg]")
print("raan:", np.rad2deg(raan), "[deg]")
print("argp:", np.rad2deg(argp), "[deg]")
print("nu:", np.rad2deg(nu), "[deg]")

Which resulted in:

p: 6613.627129899017 [km]
ecc: 0.06923948841179417
inc: 53.14716567039852 [deg]
raan: 131.4224364765848 [deg]
argp: 241.26105012632252 [deg]
nu: 184.34321702974734 [deg]

Now, these result differ a lot from the results shown on the website. At the time I am writing this, the results from the website mentioned above are:

Epoch (UTC):    04 December 2020 19:59:17
Eccentricity:   0.0001925
inclination:    51.6418°
perigee height:     418 km
apogee height:  420 km
right ascension of ascending node:  224.1133°
argument of perigee:    115.0430°
revolutions per day:    15.49154811
mean anomaly at epoch:  245.0920°
orbit number at epoch:  1846

My question being: why do my results differ so much? Did I miss something?

Answer 1

(Ok, this is my attempt of answering my own question based on the suggestions in the comments, full credits to @uhoh and @David-Hammen)

The main issue

As pointed out by @uhoh, the main issue here is parsing of the positional vectors. This part:

r1, r2, r3 = geocentric.position.km

Parses r1 = [r1.x, r2.x, r3.x], r2 = [r1.y, r2.y, r3.y] and r3 = [r1.z, r2.z, r3.z], but you want to have r1 = [r1.x, r1.y, r1.z], etc..

So, the correct parsing would be:

x, y, z = geocentric.position.km
r1 = [x[0], y[0], z[0]]
r2 = [x[1], y[1], z[1]]
r3 = [x[2], y[2], z[2]]

How could this be spotted earlier earlier:

• By checking the altitudes of the three geocentric positions.

What else could have gone wrong:

• The epoch time is 04 December 2020 19:59:17. SGP4 can be very time sensitive.
• Too much truncation
• Calculating Keplerian osculating elements 30+ minutes after the epoch (this is way too large spacing!)
• Gibb's method doesn't quite apply here. A sanity check shows that you should expect two places of accuracy.
• TLEs are not Keplerian osculating elements, Keplerian orbits (which is what rv2coe assumes) do not account for effects such as Earth's oblateness and atmospheric effects.

So let's re-write this based on these suggestions:

1. Obtaining three positional vectors over 10 minutes, where the satellite's epoch is the time of the second positional vector.

from skyfield.api import EarthSatellite, Topos, load
import numpy as np
ISS
l1 = '1 25544U 98067A   20339.83283773  .00016717  00000-0  10270-3 0  9038'
l2 = '2 25544  51.6418 224.1133 0001925 115.0430 245.0920 15.49154811 18464'
ts = load.timescale()
t1 = ts.utc(2020, 12, 4, 19, 54, 17)
t2 = ts.utc(2020, 12, 4, 19, 59, 17) # Epoch
t3 = ts.utc(2020, 12, 4, 20, 4, 17)
time = [t1, t2, t3]
satellite = EarthSatellite(l1, l2, name='ISS (ZARYA)')
for index, t in enumerate(time):
    index = index + 1
geocentric = satellite.at(t)
x, y, z = geocentric.position.km
mag = np.linalg.norm(np.array([x, y, z]))

print(f"r{index} = [{x}; {y}; {z}];")
print(f"r{index} magnitude: {mag}\n")

We get that:

r1 = [-5599.602701921242; -3436.3402635476973; -1756.6275176484512]; 
r1 magnitude: 6800.715040494088
r2 = [-4903.0245824528865; -4709.731428621669; 10.226801416237763]; # at '2020-12-04 19:59:17 UTC' <- Epoch
r2 magnitude: 6798.626682893481
r3 = [-3651.287058898108; -5449.812055367599; 1775.887224806899]; 
r3 magnitude: 6796.03737927766

2. Using a GNU Octave implementation of the Gibbs' method, we obtain the velocity of the second positional vector with:

mu = 398600;
R1 = [-5599.602701921242; -3436.3402635476973; -1756.6275176484512];
R2 = [-4903.0245824528865; -4709.731428621669; 10.226801416237763];
R3 = [-3651.287058898108; -5449.812055367599; 1775.887224806899];
r1 = norm(R1);
r2 = norm(R2);
r3 = norm(R3);
c12 = cross(R1, R2);
c23 = cross(R2, R3);
c31 = cross(R3, R1);
N = r1c23 + r2c31 + r3*c12;
D = c12 + c23 + c31;
S = R1(r2 - r3) + R2(r3 - r1) + R3*(r1 - r2);
v = sqrt(mu./norm(N)./norm(D)).*(cross(D,R2)./r2 + S);

Which will result in:

v = [3.3074; -3.4186; 5.9972]

3. Plugging the results in Poliastro's rv2coe

from poliastro.constants import GM_earth
from poliastro.core.elements import rv2coe
from astropy import units as u
import numpy as np
k = GM_earth.to(u.km ** 3 / u.s ** 2).value
v = np.array([3.3073, -3.4186, 5.9972])
r = np.array([-4903.025, -4709.731, 10.227])
p, ecc, inc, raan, argp, nu = rv2coe(k, r, v)
print("Semi-latus Rectum of Parameter:", p, "[km]")
print("Eccentriity:", ecc)
print("Inclination", np.rad2deg(inc), "[deg]")
print("Right Ascension of Ascending Node:", np.rad2deg(raan), "[deg]")
print("Argument of Perigee:", np.rad2deg(argp), "[deg]")
print("True Anomaly:", np.rad2deg(nu), "[deg]")

Will results in:

Semi-latus Rectum of Parameter: 6794.203064821717 [km]
Eccentriity: 0.001220230040456721
Inclination 51.58120030298049 [deg]
Right Ascension of Ascending Node: 223.77968914802696 [deg]
Argument of Perigee: 122.33400836496219 [deg]
True Anomaly: 237.77599772309975 [deg]

When compared to the results obtained from the website mentioned in the question, i.e.:

Eccentricity:   0.0001925
inclination:    51.6418°
right ascension of ascending node:  224.1133°
argument of perigee:    115.0430°

Might seem "close", but it is not perfect. However, plotting the obtained results:

from poliastro.bodies import Earth
from poliastro.twobody import Orbit
from poliastro.plotting.core import OrbitPlotter2D as op2d
from astropy.time import Time
epoch = Time("2020-12-4 19:59:17", scale="utc")
r = np.array([-4903.025, -4709.731, 10.227])
v = np.array([3.3073, -3.4186, 5.9972])
r = r * u.km
v = v * u.km / u.s
orbit = Orbit.from_vectors(Earth, r, v, epoch)
fig = op2d().plot(orbit)
fig.show()

Will give us:

Which is pretty close to the orbit plot shown at website from the question.