ECI <--> ECF coordinate transformation needs to account for Earth's precession

[Ian Culverwell]

Stig Syndergaard (DMI) has highlighted a problem with EUMETSAT's planned upgrade to the PPF, as follows.

Hi Dave, Ian,

I think Ian maybe got some info from Chris Marquardt on Friday 
about these things because Chris Burrows was not in office on 
Friday. I'll send some additional info after this email.

Related to ROPP, the bottom line is this:

There are different versions of ECI systems. EUMETSAT uses J2000, 
CDAAC uses something called ITOE, which is a true-of-date system. 
The rotation from J2000 -> ECF is different (and more complex) than 
the rotation from ITOE -> ECF (which is basically just a rotation 
about the z-axis). This page explains nicely about different 
systems and what true-of-date means:

​http://physics.stackexchange.com/questions/130580/how-to-determine-satellite-position-in-j2000-from-latitude-longitude-and-distan

ROPP's eci2ecf.f90 and ecf2eci.f90 only does the rotation about the 
z-axis. To process EUMETSAT data from level1a with the 
ropp_pp_occ_tool we need code to do the full J2000 -> ECF 
transformation. The Earth's z-axis in J2000 is tilted about 0.1 deg 
relative to the J2000 frame, and that difference matters. If we 
don't account for it in the ellipsoidal correction, I believe we 
will see regional biases of about 0.5% - what postponed the PPF 
day-2 upgrade.

So I'm thinking that what could be coded are extensions to the 
existing subroutines with an additional argument specifying for 
which ECI reference frame the rotation should be done. These 
subroutines can thus be used for both CDAAC and EUMETSAT data (and 
perhaps with additional implementation of other coordinate systems 
later on - we don't know what GFZ is using, WegC or China). It 
should be sufficient to code what is called the RNP matrix (see 
link above), i.e., not taking into account the so-called "Earth 
Orientation Parameters" which are minor corrections, and depends on 
regular Bulletins (see link above) that would be hard to maintain 
for ROPP - well it could be done, I suppose, if the users provides 
these themselves.

More emails to follow.

-Stig

He further says

Hi all,

I learned a bit more about coordinate systems and 
precession/nutation over the weekend. It turns out that the Earth's 
precession has the right size to be able to explain the issue. I 
did not know the difference between nutation and precession until a 
few days ago, nor did I know their sizes.

But as I understand it now, the precession has the largest 
contribution and results in that the Earth's spin axis trace out a 
cone in ECI of about 23.5 deg (the well-known tilt relative to the 
normal of the ecliptic plane) over 26000 years. After only 16 years 
since January 1, 2000, we are not far along the trace of the cone, 
but long enough so that it has resulted in a small tilt of the 
Earth's spin axis relative to the J2000 ECI z-axis of about 0.1 deg.

I found the attached chapter from a book on these things. I assume 
it is describing the J2000 ECI, but it is not totally clear. 
Anyways, plugging in 16.3 years (mid April 2016) I got the tilt to 
be 0.0907 deg (theta in eq. 2.16). I tried to confirm this by 
backing out the tilt from the NetCDf4 files from April 17 
using vectors that are given in both coordinate systems. That gave 
me about 0.0903 deg, but with a little daily variation of about 
0.0001 deg that I don't understand. Oh, well!
Maybe it is not exactly J2000 that is described in the attached 
doc, but close enough.

I made the attached drawing trying to explain the issue 
geometrically. This is also based on what Riccardo told me in 
Barcelona when we made a similar drawing on the table at the 
University restaurant.

Hope it helps understanding all this.

Cheers,
-Stig

The picture and book chapter to which Stig refers are attached.

Clearly this is important; equally clearly, we don't have the resources to do anything about it in the middle of an upgrade of ROPP. So this ticket records the problem but is not assigned to anyone.

[Ian Culverwell]

Chapter 2 "Coordinate and Time Systems"

[Ian Culverwell]

ellipsoid.png

[Ian Culverwell]

Kent Baekgaard Lauritsen has asked the RO data providers to describe the ECI systems they use for their level 1A data.

Here is the information from EUMETSAT.

Hi Kent (et al.),

with respect to the inertial reference system, we are consistent with IGS, which
means that we are using IGS08 (or I think more accurately: IGb08) reference
frame (https://igscb.jpl.nasa.gov/network/refframe.html); the corresponding
paper in GPS Solutions (2012):
ftp://igs.org/pub/resource/pubs/IGS08_The_IGS_Realization_of_ITRF2008.pdf. 

In order to transform between the above celestial reference frame we currently
follow the IERS conventions from 2003:
https://www.iers.org/IERS/EN/Publications/TechnicalNotes/tn32.html. We actually
use the reference software that was provided with that convention; it can be
downloaded from, e.g., http://62.161.69.131/iers/conv2003/conv2003_c5.html.

If you need a EUMETSAT reference, I will update the Product Format Specification
hopefully next week, and it will contain the above references.

For the future, e.g. EPS-SG, we might at some point switch to the more recent
IERS 2010 conventions
(https://www.iers.org/IERS/EN/Publications/TechnicalNotes/tn36.html). @Riccardo,
do we mention those in the ATBD/PGS?

Hope this helps,

   Chris

[Ian Culverwell]

Here is the information from Da Kuang at JPL:

All ground stations are defined in Earth Centered Earth Fixed (ECEF)  system
with a set of x, y, and z coordinates. The Earth Centered  Inertial (ECI) system
is where the equation of motion is described. The  transformation between ECEF
and ECI are 5 rotations:

1. polar motion Xp: rotation around y-axis by a small angle Yp; 
2. polar motion Yp: rotation around x-axis by a small angle Xp; These two
rotations bring the z-axiz to align with Earth's spin axis.
3. Earth rotation U: rotation around z-axis by Greenwich sidereal time  (GMT,
westward from Greenwich meridian to vernal equinox). This bring  the x-axis to
the equinox of date.
4. nutation: 
5. Precession:

These two bring the x-axis to the mean equinox, the x-axis of ECI. In JPL's GPS
products (the Flinn product), all stations are in IGb08,  aligned to ITRF2008
frame. This is the also the ECEF realization for  all orbit products
representation. Polar motion and Earth rotation  parameters are also given in
the product package. Nutation and  Precession models are given by those in IERS
Conventions (2010).

In the current FlinnR product, orbits are given in pos_goa data files  in ECEF
system. If LEO orbit and ground stations are also given in  ECEF, then
transformation may not be needed.

Kent thinks this means that ROPP would treat their Level 1A data correctly.

[Ian Culverwell]

And here is some information from Stefan Heise at GFZ:

1) Our Inertial frame is TDS (True of Date System). Here the reference  frame is
defined by the the mean equator and equinox of the given day.  It is not
referred to the epoch J2000.

2) We do the ellipsoidal correction in the Inertial frame (but  obviously TDS
avoids bias characteristics that may occur if the  Inertial frame refers to a
fixed epoch).

Stig agrees that using the TDS avoids the problem.

[Ian Culverwell]

(Stig also says that he thought Axel had already confirmed that what he sees in the data is something else.)

[Ian Culverwell]

At the ROPP9.0 DRI on 13/2/2017, and some follow up discussions, we agreed to move the deadline of this to ROPP10.0 (scheduled for release Jun 2018). It might even make its way into an intermediate release, ROPP9.1. (This issue is blocking for DMI reprocessing. They will provide the source code.)

[Ian Culverwell]

Changes to support this functionality have been introduced into the branch https://trac.romsaf.org/ropp/browser/ropp_src/branches/dev/Share/ic_sofa at r5569, r5570, r5571, r5572, r5573, r5574, r5575 and r5576.

[Ian Culverwell]

This implementation differs from DMI's in that

• SOFA is now optional. Users need not build it to get ROPP to work. This is done by saying

  COORD = ../coordinates/coordinates.f90       \
          ../coordinates/earth.f90             \
          ../coordinates/curvature.f90         \
          ../coordinates/cart2geod.f90         \
          ../coordinates/geod2cart.f90         \
          ../coordinates/impact_parameter.f90  \
          ../coordinates/occ_point.f90         \
          ../coordinates/tangent_point.f90     \
          ../coordinates/plane_coordinates.f90 \
          ../coordinates/vectors.f90
  
  if HAVE_SOFA
    COORD += ../coordinates/ecf2eci_sofa.f90   \
             ../coordinates/eci2ecf_sofa.f90   \
             ../coordinates/eci2eci_sofa.f90
  else
    COORD += ../coordinates/ecf2eci_nosofa.f90 \
             ../coordinates/eci2ecf_nosofa.f90 \
             ../coordinates/eci2eci_nosofa.f90
  endif
  

in ropp_utils/build/Makefile.am. Here, the *_sofa.f90 files are the same as DMI's, and the *_nosofa.f90 files have the same interface, but don't call the SOFA routines. If the SOFA routines are requested but unavailable, an error message to that effect is issued, saying that the code is going to use the Hoffman-Wellenhof formulae instead. (Currently the only use for the eci2eci routine is to rotate EUMETSAT POD data from its J2000 frame to the true-of-date frame, and this is done by means of the H-W formulae anyway.)

• 'configure' scripts have been generated (and tested) for ifort12, ifort15, ifort16, nagfor52, gfortran and pgf15 on linux. ('configure' is in quotes because SOFA has already been configured, for gfortran on linux. The 'configure' scripts just edit the entries in the given makefile, using sed. The organisation of this rewriting has been changed a bit.)

• Changes to effect the building of the SOFA libraries have been made to the 'ROPP' version of buildpack, not the DMI one, which has some differences.

[Ian Culverwell]

It all seems to build OK. How big are the effects? Test this out by looking at the changes to the POD which follow from applying each of the three available eci2eci transformations to the EUMETSAT POD (in J2000 frame) in the ROPP test file ropp_io/data/eum_test.n4 (which was originally derived from GRAS_1B_M02_20120909000057Z_20120909000118Z_N_T_20141111112929Z_G15_NN.nc.

We find, at the second point in the closed loop timeseries,

BEFORE eci2eci: data_cl%Lev1a%r_gns(2,:) =    -3.85387E+06   -2.17145E+07   -1.46448E+07
AFTER  eci2eci: data_cl%Lev1a%r_gns(2,:) =    -3.77418E+06   -2.17253E+07   -1.46495E+07

BEFORE eci2eci: data_cl%Lev1a%v_gns(2,:) =     2.43448E+03   -1.96950E+03    2.30229E+03
AFTER  eci2eci: data_cl%Lev1a%v_gns(2,:) =     2.43722E+03   -1.96259E+03    2.30529E+03

BEFORE eci2eci: data_cl%Lev1a%r_leo(2,:) =    -3.64800E+06    5.05684E+06   -3.60806E+06
AFTER  eci2eci: data_cl%Lev1a%r_leo(2,:) =    -3.65789E+06    5.04647E+06   -3.61256E+06

BEFORE eci2eci: data_cl%Lev1a%v_leo(2,:) =     3.18962E+03   -2.21595E+03   -6.34154E+03
AFTER  eci2eci: data_cl%Lev1a%v_leo(2,:) =     3.20371E+03   -2.20688E+03   -6.33760E+03

So the changes to the positions, especially, are of a reasonable size.

If we try with SOFA=.TRUE. and CIO=.FALSE.:

BEFORE eci2eci: data_cl%Lev1a%r_gns(2,:) =    -3.85387E+06   -2.17145E+07   -1.46448E+07
AFTER  eci2eci: data_cl%Lev1a%r_gns(2,:) =    -3.77226E+06   -2.17258E+07   -1.46491E+07

BEFORE eci2eci: data_cl%Lev1a%v_gns(2,:) =     2.43448E+03   -1.96950E+03    2.30229E+03
AFTER  eci2eci: data_cl%Lev1a%v_gns(2,:) =     2.43729E+03   -1.96238E+03    2.30540E+03

BEFORE eci2eci: data_cl%Lev1a%r_leo(2,:) =    -3.64800E+06    5.05684E+06   -3.60806E+06
AFTER  eci2eci: data_cl%Lev1a%r_leo(2,:) =    -3.65813E+06    5.04615E+06   -3.61277E+06

BEFORE eci2eci: data_cl%Lev1a%v_leo(2,:) =     3.18962E+03   -2.21595E+03   -6.34154E+03
AFTER  eci2eci: data_cl%Lev1a%v_leo(2,:) =     3.20405E+03   -2.20679E+03   -6.33746E+03

So similar but different differences.

If we try with SOFA=.TRUE. and CIO=.TRUE.:

BEFORE eci2eci: data_cl%Lev1a%r_gns(2,:) =    -3.85387E+06   -2.17145E+07   -1.46448E+07
AFTER  eci2eci: data_cl%Lev1a%r_gns(2,:) =    -3.83537E+06   -2.17148E+07   -1.46491E+07

BEFORE eci2eci: data_cl%Lev1a%v_gns(2,:) =     2.43448E+03   -1.96950E+03    2.30229E+03
AFTER  eci2eci: data_cl%Lev1a%v_gns(2,:) =     2.43157E+03   -1.96945E+03    2.30540E+03

BEFORE eci2eci: data_cl%Lev1a%r_leo(2,:) =    -3.64800E+06    5.05684E+06   -3.60806E+06
AFTER  eci2eci: data_cl%Lev1a%r_leo(2,:) =    -3.64345E+06    5.05676E+06   -3.61277E+06

BEFORE eci2eci: data_cl%Lev1a%v_leo(2,:) =     3.18962E+03   -2.21595E+03   -6.34154E+03
AFTER  eci2eci: data_cl%Lev1a%v_leo(2,:) =     3.19762E+03   -2.21609E+03   -6.33746E+03

Different again.

Note that if we build without the SOFA libraries, we find, in the last two cases:

BEFORE eci2eci: data_cl%Lev1a%r_gns(2,:) =    -3.85387E+06   -2.17145E+07   -1.46448E+07
WARNING (from ropp_io_read_eumdata):  eci2eci called with SOFA = .TRUE., but SOFA library is unavailable. Will use Hoffman-Wellenhof formulae instead.
AFTER  eci2eci: data_cl%Lev1a%r_gns(2,:) =    -3.77418E+06   -2.17253E+07   -1.46495E+07

BEFORE eci2eci: data_cl%Lev1a%v_gns(2,:) =     2.43448E+03   -1.96950E+03    2.30229E+03
WARNING (from ropp_io_read_eumdata):  eci2eci called with SOFA = .TRUE., but SOFA library is unavailable. Will use Hoffman-Wellenhof formulae instead.
AFTER  eci2eci: data_cl%Lev1a%v_gns(2,:) =     2.43722E+03   -1.96259E+03    2.30529E+03

BEFORE eci2eci: data_cl%Lev1a%r_leo(2,:) =    -3.64800E+06    5.05684E+06   -3.60806E+06
WARNING (from ropp_io_read_eumdata):  eci2eci called with SOFA = .TRUE., but SOFA library is unavailable. Will use Hoffman-Wellenhof formulae instead.
AFTER  eci2eci: data_cl%Lev1a%r_leo(2,:) =    -3.65789E+06    5.04647E+06   -3.61256E+06

BEFORE eci2eci: data_cl%Lev1a%v_leo(2,:) =     3.18962E+03   -2.21595E+03   -6.34154E+03
WARNING (from ropp_io_read_eumdata):  eci2eci called with SOFA = .TRUE., but SOFA library is unavailable. Will use Hoffman-Wellenhof formulae instead.
AFTER  eci2eci: data_cl%Lev1a%v_leo(2,:) =     3.20371E+03   -2.20688E+03   -6.33760E+03

as expected.

[Ian Culverwell]

plot_POD.gif

[Ian Culverwell]

These results are compared more closely, for all the (closed loop) data points in the file, here:

This figure compared the (x, y, z) co-ordinates of r_GNS in the true-of-date ECI frame defined by the simple formulae of Hoffman-Wellenhof et al, compared to the co-ordinates in the original J2000 ECI frame. Ditto for the equinox-based SOFA calculations, and the CIO-based SOFA calculations. And ditto ditto for v_GNS, r_LEO and v_LEO.

Generally, the position differences are around 10 km, and the velocity differences are around 10 m/s. I don't know if this is as expected. What probably is expected is that the radius |r| and speed |v| remain the same after these rotations, and this is indeed true to a high degree of accuracy:

Processing r_gns
max(abs(|r|_IERS_H-W - |r|_J2000)) =    7.4505806e-09
max(abs(|r|_SOFA_eqx - |r|_J2000)) =    7.4505806e-09
max(abs(|r|_SOFA_CIO - |r|_J2000)) =    7.4505806e-09

Processing v_gns
max(abs(|v|_IERS_H-W - |v|_J2000)) =    9.0949470e-13
max(abs(|v|_SOFA_eqx - |v|_J2000)) =    9.0949470e-13
max(abs(|v|_SOFA_CIO - |v|_J2000)) =    9.0949470e-13

Processing r_leo
max(abs(|r|_IERS_H-W - |r|_J2000)) =    1.8626451e-09
max(abs(|r|_SOFA_eqx - |r|_J2000)) =    1.8626451e-09
max(abs(|r|_SOFA_CIO - |r|_J2000)) =    2.7939677e-09

Processing v_leo
max(abs(|v|_IERS_H-W - |v|_J2000)) =    1.8189894e-12
max(abs(|v|_SOFA_eqx - |v|_J2000)) =    2.7284841e-12
max(abs(|v|_SOFA_CIO - |v|_J2000)) =    1.8189894e-12

(Timeseries of these quantities are also plotted in grey on the figure, but you can't distinguish them from zero by eye.)

[Ian Culverwell]

(I forgot to point out that the base code for this is contained in Hans's r5331 and r5332.)

[Ian Culverwell]

Hans Gleisner at DMI commented on these results:

Hi Ian,

I've now checked the test results you present in the 
ticket.

12 years of precession of the Earths rotational axis 
with respect to an inertial frame gives a difference of 
about 10 kilometers in LEO satellite position.

The expected differences in LEO satellite position 
between [SOFA=.false] and [SOFA=.true.,CIO=.false.] are a 
few hundred meters. This difference is due to the fact 
that with SOFA=.false. we use a simple textbook formula
that doesn't include the nutation. With SOFA=.true., the 
nutation is included, which explains the difference.

The option [SOFA=.true, CIO=.true.] is different as it 
uses an entirely different X-Y coordinate system. The 
differences in the X-Y plane with respect to [SOFA=.true, 
CIO=.false.] is on the order of 10 kilometers.

My conclusion is that your test results are consistent 
with what we expect.

        /Hans

Good.

[Ian Culverwell]

Note that the Cray intrinsic compiler ftn cannot cope with fortran file extension .for, so we need to include the following in sofa_configure_ftn_linux:

for file in $(find ${deps_src}/sofa -name "*.for") ; do
  echo cp $file $(echo $file |sed -es/'\.for'/'\.f'/)
       cp $file $(echo $file |sed -es/'\.for'/'\.f'/)
done
cat $mfile |sed -es/'\.for'/'\.f'/ > $HOME/makefile.temp
mv $HOME/makefile.temp $mfile

at r5636 and r5661.

[Ian Culverwell]

Note that the default behaviour, which does not use SOFA/CIO, does not affect the behaviour of eci2ecf or ecf2eci. There is an impact on the processing of EUM level 1a data, through a call to the new routine eci2eci, but even this does not(yet) take advantage of the SOFA/CIO routines, but uses the Hoffman-Wellenhof formulae instead.

[Ian Culverwell]

Passes automatic make tests OK, after altering the reference file at r5632 and r5633.

[Ian Culverwell]

At ROPP10.0 we should probably include a test of a big load of EUM lev1a data as part of the test folder. But for ROPP9.1, we can close the ticket.

[Ian Culverwell]

(For the record, the ic_sofa branch mods were merged into the ROPP91_prototype at r5628.)

[Ian Culverwell]

Reclosing ticket.