How to get earth-centered, earth-fixed coordinates from Skyfield?

Question

I recently started using Skyfield (python) to find (x, y, z) coordinates of things in space. Now I'd like to transform between Earth-Centered, Earth-Fixed (rotating) and Earth-Centered Inertial (non-rotating) coordinates, since "rocket speeds" may sometimes be webcasted in earth-fixed (rotating) frames . See this answer and this answer

This is what I have so far. The barycentric positions are straight from Skyfield, but I'm just subtracting the earth position from the moon's position to get earth-centered inertial position - I think this is OK for non-precision results, I don't know if it is really that easy or if there is something subtle I'm missing.

However, I'm wondering if there is an "official" Earth-fixed (rotating) frame, and if Skyfield has a method to convert to it, or to any rotating frame in general?

Python:

from skyfield.api import load, JulianDate

de421 = load('de421.bsp')

earth = de421['earth']
moon  = de421['moon']
slc40 = earth.topos('28.562 N', '80.577 W')

jd = JulianDate(utc=(2015,12, 22, 1, 48, 0))

epos     = earth.at(jd).position.km # barycentric (ICRS) position of earth
mpos     = moon.at(jd).position.km  # barycentric (ICRS) position of moon
slc40pos = slc40.at(jd).position.km # barycentric (ICRS) position of SLC-40

mpos_ec     = mpos     - epos # earth centered (inertial) position of moon
slc40pos_ec = slc40pos - epos # earth centered (inertial) position of SLC-40

Answer 1

Skyfield ≥1.34

The most recent versions of Skyfield introduce an explicit itrs reference frame that rotates with the Earth. You can learn about it starting here:

• https://rhodesmill.org/skyfield/coordinates.html#eci-versus-ecef-coordinates
• https://rhodesmill.org/skyfield/coordinates.html#geographic-itrs-latitude-and-longitude

Skyfield <1.34

First: the main problem with your technique of subtracting the two positions is that it does not account for light-time travel. In the case of the Moon observation that you set up in your code, the error is only about 38 km, which might be well within your tolerances. The way to ask for a relative position that is properly light-time backdated is through the observe() method:

mpos_ec = earth.at(jd).observe(moon).position.km

Second: You are right that Skyfield currently has no built-in first-class support for a reference frame that rotates with the Earth. I should probably add one. A workaround that should get you up and running immediately, though, is to tell Skyfield that you want to observe:

• From an Earth location, so the Skyfield does the rotation calculation that you want.
• But from an Earth location whose latitude and longitude are zero, so that Skyfield applies no further rotations beyond moving into the Earth's rotating reference frame.
• And from a location whose elevation is negative enough to place it exactly at the center of the Earth, so that you are not applying an offset to the position of the Moon or whatever.

The one trick is that we have to do the rotation ourselves because, as you note, the way Skyfield is currently written it's only willing to apply the rotation if it's about to reduce the coordinates to polar alt/az for you. So try out the following:

# Create an observer at the center of the Earth,
# rotated zero degrees from the Earth's reference frame:

gcrs = earth.topos(
    latitude_degrees=0,
    longitude_degrees=0,
    elevation_m=-6378136.6,
    )

# Ask where the observer was at a given time:

g = gcrs.at(jd)

# Then ask for a relative position as usual:

p = g.observe(moon).position.km
print p

# Now the tricky part: each observer.at() has
# an attribute `altaz_rotation` which is the rotation
# matrix you want to apply:

print einsum('ij...,j...->i...', g.altaz_rotation, p)

Let me know if the numbers you get back seem to work! And once I have the time-system rewrite finished up and can release Skyfield 1.0, I'll turn my attention to make this computation a bit more convenient.