Why isn't this orbit on the equator?

[troyrock]

I have the following TLE which is a nearly circular orbit with no angle of inclination. I would expect the satellite to circle exactly around the equator:

1 70000U 99999ZZZ 21079.40069444 .00022359 00942-5 37766-3 0 1009 2 70000 0.0000 0.0000 0000100 0.0000 0.0000 14.87482098 12341

I am using the following code:

from skyfield.api import load, EarthSatellite
from datetime import timedelta

sat = EarthSatellite(tle1, tle2)
minutes = 0
ts = load.timescale(builtin=True)
time = ts.from_datetime(sat.epoch.utc_datetime() + minutes * timedelta(minutes = 1))
geocentric = sat.at(time)
print(geocentric.position.km)

Gives me the result of:

[6980.64196792 -33.11916156 -14.1663118 ]

I'm getting a z values of -14.1663118 when I expected it to be on the equator. If I continue to follow the orbit, it oscillates between this value and positive 14. When I use the sgp4 module (see following code), I get exactly 0.0 for the z coordinate. Are the two using different models of the Earth?

from sgp4.api import Satrec
from skyfield.api import load
from datetime import timedelta

sat = Satrec.twoline2rv(tle1, tle2)
minutes = 0
ts = load.timescale(builtin=True)
time = ts.from_datetime(sat.epoch.utc_datetime() + minutes * timedelta(minutes = 1))
jd = time.whole
fr = time.tai_fraction - times[k]._leap_seconds() / 86400.0
e_sat, pos_sat, v_sat = sat.sgp4(jd, fr)
print(pos_sat)

Gives me the result of: (6980.734907411272, 0.000136108758132057, 0.0)

In general, I would expect these numbers to be exactly the same.

[brandon-rhodes]

The two libraries use different references frames. Try looking them up and see what you can find! If either library leaves you guessing about its reference frame, then let me know what documentation you consulted, and we'll get it improved.

[troyrock]

I am not as familiar as I need to be with the different coordinate systems! It appears to me that the sgp4 library gives me a TEME coordinate system answer, just like the one the TLEs use and geocentric gives me a GCRS coordinate system answer. Is there a way to convert between the two and how do these relate to ITRF? I can see that there is a deprecated function to convert TEME to ITRF but I didn't understand how to "use the TEME and ITRS frame objects instead."

I would like to plot the orbits using something similar to:

fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
ax.stock_img()
subsat = wgs84.subpoint(geocentric)
plt.scatter(subsat.longitude.degrees, subsat.latitude.degrees, transform=ccrs.PlateCarree(), color='red')
plt.show()

Is this the correct way to plot GCRS coordinates or do I have to do some transform to keep track of the earth spinning?

When I plot GCRS together with TEME, the GCRS orbits leave a gap equal to the rotation of the Earth but the TEME come back to their starting point (which I think is incorrect since the earth has moved during the orbit so it won't rejoin the same latitude.)

[brandon-rhodes]

Is there a way to convert between the two and how do these relate to ITRF?

(The ITRF merely helps define the ITRS reference frame, so I'll try to use ITRS when I type the acronym. We'll see how successful I am.)

The TEME and GCRS are similar in that both coordinate systems are fixed in space and don't rotate every day; the difference is that the GCRS uses permanent axes pointing in space, whereas the TEME axes gradually shift with precession (the “mean equinox" mentioned in its name).

The ITRS instead takes as its reference the Earth's surface, and the whole coordinate system spins around each day. (But something like, say, Paris, keeps the same x,y,z all day.)

As for converting, yes, there's a way! The Skyfield position docs show, for example, converting a Skyfield position to ITRS:

https://rhodesmill.org/skyfield/positions.html#coordinates-in-other-reference-frames

If you need to convert between TEME and GCRS, you'll want to take a look at the TEME object shown here:

https://rhodesmill.org/skyfield/api.html#earth-satellites

But Skyfield should convert the coordinates out of the TEME frame for you, so most users don't need that object? If you find yourself needing it, let me know why, and maybe I'll need to expand the docs some to cover your use case.

Is this the correct way to plot GCRS coordinates or do I have to do some transform to keep track of the earth spinning?

Those will be ITRS coordinates, not GCRS, because the latitude/longitude system rotates every day in a circle. For GCRS coordinates, you'd just ask the position for its x,y,z as shown at:

https://rhodesmill.org/skyfield/positions.html#position-attributes

Those coordinates should be fixed in space and the only motion will be the satellite's orbit, not the motion of the Earth beneath it.

Let me know where I might update the docs to explain this where users can find it. Thanks!

[troyrock]

I'm very confused at this point. I did the following for all the points in one revolution of an orbit:

geocentric = newsat.at(times)
newlatlon = geocentric.frame_latlon(itrs)
subsat = wgs84.subpoint(geocentric)
plt.plot(subsat.latitude.degrees - newlatlon[0].degrees)

I expected more of a difference, like at a completely different location but they are essentially the same.

[brandon-rhodes]

The ITRS is the frame of reference that rotates with the Earth, providing the x,y,z coordinates that go with our latitude,longitude coordinate system. I think your code is asking for latitude in two slightly different ways, and so would be expected to return very nearly the same angles?

[troyrock]

Thanks for your patience and help. I feel like I understand much better what I'm doing and how to proceed.

[brandon-rhodes]

Let me know what documentation tweaks might have made things easier, and I'll consider whether there are edits or additions I might make.

[troyrock]

Nearly all of the difficulty for me boiled down to understanding the different reference frames involved. I wouldn't be surprised if I still wasn't completely squared on everything. I'm fairly new to this so I fully expected to have a steep learning curve. For me a short tutorial of the different coordinate systems and where they appear would have been helpful. For instance, I found the following somewhere on the web that was helpful to me:

The difference between GCRF and TEME resides in the fact that in the GCRF frame, the X-axis and Z-axis are aligned respectively with the mean equinox and rotation axis of Earth at 12:00 Terrestrial Time on the 1st of January, 2000, while in the TEME frame they are aligned with the mean equinox and rotation axis at the time of the corresponding TLE. Due to the change of the direction of the vernal equinox and the rotation axis over time, coordinates in the two frames differ slightly.

It wasn't clear to me that when you ask for x, y, z coordinates you get GCRF coordinates but when you ask for latitude and longitude, they are aligned to the ITRS coordinate system (so that you can plot them correctly). You are doing exactly the best thing but it wasn't clear to me that there were all of these different coordinate systems involved.

If I had understood the difference between TEME and GCRF, I wouldn't have been surprised that they give different answers for x, y, and z but I still needed to understand that asking for latitude and longitude is yet another reference frame and might have still been surprised that the latitude and longitude given from geocentric.frame_latlon(itrs) is not the same as you would get if you tried to convert x, y, z coordinates to latitude and longitude using the readily available formulas (i.e. https://web.archive.org/web/20080920155754/http://www.ferris.edu/faculty/burtchr/papers/cartesian_to_geodetic.pdf) because of the additional conversion from GCRF to ITRS.

You have been very careful about the code and all of the different reference systems. I have looked through several areas of the code that are both clean and very well documented. I'm very impressed with the tool: it's a pleasure to have such annoying complexities completely taken care of in the correct way.