Possible Inertial coordinate system discrepancy with astropy

[curtispuetz]

I wanted to verify that the position vector I get from using skyfield's EarthSatellite was the same as what I would get from taking skyfields geodetic coordinates and then using astropy to convert to gcrs. There is a slight discrepancy in the vectors, which I think is non negligible:

from astropy import coordinates as coords
from astropy.time import Time
from skyfield.api import load, EarthSatellite, utc
from datetime import datetime

line1 = '1 44031U 98067PX  19083.14584174  .00005852  00000-0  94382-4 0  9997'
line2 = '2 44031  51.6393  63.5548 0003193 165.0023 195.1063 15.54481029  8074'
satellite = EarthSatellite(line1, line2)
ts = load.timescale()
time_track = datetime(2019, 3, 24, 18, 35, 1, tzinfo=utc)
t = ts.utc(time_track)
geo = satellite.at(t)
subpoint = geo.subpoint()

now = Time(time_track)
ecef = coords.EarthLocation.from_geodetic(subpoint.longitude.degrees, subpoint.latitude.degrees, subpoint.elevation.m)
gcrs = ecef.get_gcrs(obstime=now).cartesian.xyz.value
print(geo.position.m)
print(gcrs)

result:

[ 4153133.07088156 -1095531.83670525 -5255198.0848858 ]
[ 4153155.43621461 -1095479.82910933 -5255191.81122256]

[brandon-rhodes]

To narrow the question, we can observe that where the coordinates come from doesn't affect the result. If we take the same x,y,z coordinates as raw AU numbers and ask Skyfield "what's the subpoint" and then ask it "what's the x,y,z of that subpoint", Skyfield can convert back and forth with an error that's a bit under a millimeter:

from astropy import coordinates as coords
from astropy.time import Time
from skyfield.api import load, utc
from skyfield.positionlib import Geocentric
from datetime import datetime

ts = load.timescale()
time_track = datetime(2019, 3, 24, 18, 35, 1, tzinfo=utc)
t = ts.utc(time_track)
geo = Geocentric([2.776As19798e-05, -7.32317801e-06, -3.51288294e-05], t=t, center=399)
subpoint = geo.subpoint()
print(geo.position.m)
print(subpoint.at(t).position.m)
print(subpoint.longitude.degrees, subpoint.latitude.degrees, subpoint.elevation.m)

#time_track = time_track.replace(second=0, microsecond=720060)

now = Time(time_track)
ecef = coords.EarthLocation.from_geodetic(subpoint.longitude.degrees, subpoint.latitude.degrees, subpoint.elevation.m)
gcrs = ecef.get_gcrs(obstime=now).cartesian.xyz.value
print(gcrs)

399
[ 4153133.06449641 -1095531.83705306 -5255198.07842356]
[ 4153133.06474804 -1095531.83711916 -5255198.07821222]
-115.19912019925971 -50.81491889072899 421850.1753570454
[ 4153155.4298366  -1095479.82942796 -5255191.80476125]

The question is why the AstroPy code, given the same latitude and longitude and elevation, doesn’t produce these same x,y,z coordinates. Playing around with the numbers a bit, it looks like AstroPy is using the old WGS84 values for the Earth's radius and flattening, but plugging those into Skyfield constants.py, I only see a difference in the output values of about 4cm — much smaller than this difference of >20 meters.

As you can see in the commented-out line, I also tried adjusting the time to see if it was a calendar difference of some sort, but finding a moment that brought the x-coordinates quite close didn't make the y and z match exactly.

Is there a third tool we could use to turn the latitude, longitude, and elevation into x,y,z coordinates to decide between the two libraries? I looked at HORIZONS but didn't immediately see a tool to do it.

[curtispuetz]

We could use pysofa I think. One of my colleagues wrote a function that outputs a rotation matrix from an inertial coordinate system to the ECEF coordinates with pysofa. I can use this function combined with pysofa's gd2gc function to turn latitude, longitude, elevation to this inertial system (which may be a little different than gcrs? I'm not well versed in these inertial systems). Continuing from the code I posted before:

rot_matrix = icrf_to_fixed(time_track)

import pysofa
import numpy as np
gc = pysofa.gd2gc(1, subpoint.longitude.radians, subpoint.latitude.radians, subpoint.elevation.m)
inertial_pysofa = rot_matrix.T @ np.array(gc).reshape(3)
print(geo.position.m)
print(inertial_pysofa)
print(gcrs)

[ 4153133.07088156 -1095531.83670525 -5255198.0848858 ]
[ 4153155.2741622  -1095448.38028367 -5255198.49497931]
[ 4153155.41207798 -1095479.78970774 -5255191.83854414]

The icrf_to_fixed function:

import numpy as np
import pysofa
import datetime

def icrf_to_fixed(epoch: datetime.datetime, seconds: np.ndarray=0.0):
    """
    Calculate rotation matrices for transforming vectors from the ICRF (International Celestial Reference Frame) to
    a Fixed frame that follows the attitude of the Earth. These matrices should work well for STK ephemeris (.e) and
    attitude (.a) files in the ICRF and Fixed frames.

    The UT1-UTC offset, the CIP offsets, and the polar motion corrections are all ignored to match default STK
    settings. These corrections are only necessary for arcsecond precision.

    The input date is in UTC.

    The returned array has shape (3, 3) if seconds is a scalar, and (len(seconds), 3, 3) if seconds is a vector.
    """
    sec = np.atleast_1d(seconds)

    djmjd0, date = pysofa.cal2jd(epoch.year, epoch.month, epoch.day)
    time0 = (60. * (60. * epoch.hour + epoch.minute) + epoch.second) / 86400.
    dat = pysofa.dat(epoch.year, epoch.month, epoch.day, time0)

    rc2ti = np.zeros((len(sec), 3, 3))
    for i, s in enumerate(sec):
        time = (60. * (60. * epoch.hour + epoch.minute) + epoch.second + s) / 86400.
        utc = date + time
        tai = utc + dat / 86400.
        tt = tai + 32.184 / 86400.

        # cip and cio, IAU 2006/2006A
        x, y, s = pysofa.xys06a(djmjd0, tt)

        # gcrs to cirs matrix
        rc2i = pysofa.c2ixys(x, y, s)

        # earth rotation angle
        era = pysofa.era00(djmjd0 + date, time)

        # form celestial-terrestrial matrix (no polar motion yet)
        rc2ti[i, :, :] = pysofa.rz(era, rc2i)

    if isinstance(seconds, float):
        return rc2ti[0, :, :]
    else:
        return rc2ti

So, all the results are different? Could be that this pysofa inertial system is different, or that my colleague did some approximations along the way.

[brandon-rhodes]

I think that the "UT1-UTC offset" difference it mentions could be a contributing factor. It says "These corrections are only necessary for arcsecond precision" but what does that mean for an x,y,z coordinate? How big is an arcsecond in meters on the ground? If I'm doing my math right:

6378000 * 2 * pi / 360 / 60 / 60

30.921417

— then an arcsecond measured from the geocenter is 30 meters at the Earth's surface. That's around the precision of these differences, I think?

[brandon-rhodes]

Pending a response, I am going to close this issue for now. Please re-open it if in the future you have further information, or a reply to my error bar computation above. Thanks!