Frame conversion between ITRF and EME2000

ipark · November 15, 2018, 8:13pm

Hi, Orekit users,

It seems that a frame conversion would be the most popular topic nowadays.
I also have a question about the conversion, especially, between ITRF and EME2000.

Here is what I’m usually using for converting position and velocity vectors from EME2000 to ITRF:

inertial_frame = FramesFactory.getEME2000()
earth_frame = FramesFactory.getITRF(IERSConventions.IERS_2010, False)
t_epoch = AbsoluteDate(int(t0[0]), int(t0[1]), int(t0[2]), int(t0[3]), int(t0[4]), float(t0[5]), TimeScalesFactory.getUTC())
position = Vector3D([2127263.8228368093, 5718032.172429709, 3480933.682526727])
velocity = Vector3D([-154.1568002046, -3894.204237608, 6449.7773756905])
pv = PVCoordinates(position, velocity)
pv_itrf = inertial_frame.getTransformTo(earth_frame, t_epoch).transformPVCoordinates(pv)

And what I’m doing now is if it is possible to use a ‘rotation matrix between ITRF and EME2000.’ Followings are what I’m trying for that.

inertial_frame = FramesFactory.getEME2000()
earth_frame = FramesFactory.getITRF(IERSConventions.IERS_2010, False)
t_epoch = AbsoluteDate(int(t0[0]), int(t0[1]), int(t0[2]), int(t0[3]), int(t0[4]), float(t0[5]), TimeScalesFactory.getUTC())
position = Vector3D([2127263.8228368093, 5718032.172429709, 3480933.682526727])
velocity = Vector3D([-154.1568002046, -3894.204237608, 6449.7773756905])
pv = PVCoordinates(position, velocity)
R_eme2000_to_itrf_JArray = eme2000.getTransformTo(earth_frame, t_epoch).getRotation().getMatrix()
R_eme2000_to_itrf_components = [orekit.pyhelpers.JArray('double').cast_(x) for x in R_eme2000_to_itrf_JArray]
R_eme2000_to_itrf = np.array([
        [R_eme2000_to_itrf_components[0][0], R_eme2000_to_itrf_components[0][1], R_eme2000_to_itrf_components[0][2]],
        [R_eme2000_to_itrf_components[1][0], R_eme2000_to_itrf_components[1][1], R_eme2000_to_itrf_components[1][2]],
        [R_eme2000_to_itrf_components[2][0], R_eme2000_to_itrf_components[2][1], R_eme2000_to_itrf_components[2][2]]
        ])

rotated_pos_itrf = np.matmul(rot_matrix_eme2000_to_itrf, np.array([position.getX(), position.getY(), position.getZ()]))
rotated_vel_itrf = np.matmul(rot_matrix_eme2000_to_itrf, np.array([velocity.getX(), velocity.getY(), velocity.getZ()]))

I expected I could get the same position and velocity in ITRF frame, but it was not true for the velocity.
Here are differences in position and velocity components w.r.t. the first conversion.

dposition = rotated_pos_itrf-(position from 'pv_itrf')
dvelocity = rotated_vel_itrf-(velocity from 'pv_itrf')
('dposition: ', array([ 0.00000000e+00,  9.31322575e-10, -4.65661287e-10]))
('dvelocity: ', array([ 4.24467562e+02, -1.32768188e+02, -5.92213273e-04]))

As shown above, the rotation matrix seems to work for the position, but velocity.

Does anyone have ideas or suggestions or advice?
Hope I can find the way to use the matrix for the conversion.

Thank you in advance for your concern and help.
Have a great day.

Best,

:: inkwan

evan.ward · November 15, 2018, 9:12pm

The issue is that treating a Transform as a rotation matrix ignores the motion between the reference frames, in this case the rotation rate of ITRF w.r.t. EME2000. See how the ω x r term is applied in [1].

What do you need the matrix for? Perhaps then I could make a suggestion.

[1] https://www.orekit.org/site-orekit-9.2/xref/org/orekit/utils/AngularCoordinates.html#L594

luc · November 16, 2018, 9:31am

As Evan wrote, the transformPVCoordinates method in the Transform class applies velocity composition (and also acceleration composition for that matter). The reason is that the semantics of this method is to fully convert the state and get something consistent, i.e. be fore conversion the velocity is the time derivative of the position, and after conversion the converted velocity is the time derivative of the converted position. This
property is mandatory for the shiftedBy() methods to work. So the velocity is the velocity relative to the frame.

Here is an example that should explain it more intuitively. Consider a perfect geostationary satellite and ignore all small motions (precession, nutation and pole wandering), keeping only the Earth rotation about its Z axis. In the inertial frame, its position and velocity are some vectors in the inertial XY plane and the velocity vector norm is about 3075.76m/s. If we use transformPVCoordinates to convert it to Earth frame, we will get a position and velocity in the terrestrial XY plane (which with our simplified hypotheses is the same as the inertial XY plane), but the velocity norm will be 0 because the satellite is geostationary and in Earth frame its position is fixed, so the derivative of the position must be 0! The computation will correctly find this 0 vector because when it converts the inertial velocity, it first applies the rotation (hence not changing its norm) and then subtract the ω x r term which will cancel the first term.

These semantics may not suit your needs. If what you want is the inertial velocity vector just projected into the terrestrial frame but not converted to a terrestrial velocity (i.e. in the geostationary case you want a 3075.76m/s vector pointing towards the local East direction), then you should call transformPosition(pv.getPosition()) and transformVector(ps.getVelocity()) and reconstruct the PVCoordinates from these two vectors. The transformVector method only applies the rotation and nothing else (no translation, no derivatives). Beware however that what you get has a very special meaning and has hidden inconsistencies: typically you cannot use shiftedBy() on this object since it would for example move your geostationary satellite at 3075.76m/s to the East when in fact it should not.

ipark · November 19, 2018, 6:52pm

Thank you for the responses Evan and Luc.

I’m trying to apply your advice for my tests to understand those clear.
Hopefully, I will figure them out asap.

There were two reasons why I started to think about this conversion.

I was curious if Orekit provided multiple options for this conversion.
I’m considering if there’s a simple method to project a covariance matrix from ITRF to EME2000, or vice versa.
The rotation matrix was tested for that purpose, and I found Orekit rotated position correctly, but velocity.

Thank you very much again for your advice. And I might ask the following questions (or test results) about this conversion.

Have a great day.

:: Inkwan

MaximeJ · November 20, 2018, 12:25pm

Hi @ipark,

To project a covariance matrix from ITRF to EME2000 I suggest you use the method getJacobian from the Transform class.

If you only use the rotation matrix you will miss a part of the Jacobian.

Suppose you already have:

covItrf: a double[][] covariance matrix (6x6) in ITRF frame;
date: an AbsoluteDate corresponding to your covariance matrix.

Consider the following Java code:

// Frames definition
final Frame j2000 = FramesFactory.getEME2000();
final Frame itrf  = FramesFactory.getITRF(IERSConventions.IERS_2010, false);

// Jacobian from ITRF to J2000 at date
final double[][] jacItrfToJ2000 = new double[6][6];
itrf.getTransformTo(j2000, date).getJacobian(CartesianDerivativesFilter.USE_PV, jacItrfToJ2000);
        
// Covariance transformation, using Hipparchus RealMatrix class to perform the multiplication
final RealMatrix pItrf = MatrixUtils.createRealMatrix(covItrf);
final RealMatrix jItrfToJ2000 = MatrixUtils.createRealMatrix(jacItrfToJ2000);

// pJ2000 contains the covariance matrix in J2000        
final RealMatrix pJ2000 = jItrfToJ2000.multiply(pItrf.multiplyTransposed(jItrfToJ2000));

Have a great day too,
Maxime