Conversion of satellite coordinates from ITRF to J2000

Question

I have coordinates of various satellites in two coordinate systems:

1. Cartesian coordinates in the international terrestrial reference frame (ITRF)
2. RA / DEC in J2000 epoch, as derived from plate solving using a star catalogue

I believe both positions are accurate to about 1 arcsecond and times are accurate to 1 ms, and I would like to verify that. Thus I will need to convert one position to the other system. I decided to try to convert (1) to (2), i.e. cartesian to RA/DEC.

That's what I did (code snippet below):

1. I converted the cartesian coordinates to spherical ones in the same coordinate system, using the Greenwich sidereal time for the hour angle / RA calculation
2. Converted epoch of observation system to J2000 epoch using pyephem
3. Corrected the annual velocity aberration using pyastronomy
4. Actually as a first step, I corrected for light travel time between satellite and observer, i.e. I took the cartesian coordinates from a few ms earlier than I have my J2000 coordinates for (which I got from an image)

Here the (python) code:

vectorX = x - siteX
vectorY = y - siteY
vectorZ = z - siteZ

# Distance
self.cur_range = sqrt(vectorX ** 2 + vectorY**2 + vectorZ**2)

# Declination
dec = atan(vectorZ / sqrt(vectorX**2 + vectorY**2)) * 180. / pi

# Right Ascension
angle_1 = atan2(vectorY, vectorX) * 180. / pi           
if angle_1 < 0:
    angle_1 = 360 + angle_1

LST = self.observer.sidereal_time()
GST = LST - self.observer.lon

RA = angle_1 + GST * 180. / pi

# conversion to J2000
coord_current = ephem.Equatorial(RA / 180. * pi, dec / 180. * pi, epoch=time)
coord_J2000 = ephem.Equatorial(coord_current, epoch = ephem.J2000)

RA = (coord_J2000.ra * 180. / pi) / 15.
dec = coord_J2000.dec * 180. / pi

# Annual aberration
if pyastro_available:
    delta_RA, delta_dec = correct_aberration(time, RA, dec)
    RA -= delta_RA
    dec -= delta_dec

return RA, dec

(The code as shown here is not complete, but it should give an idea of what I do)

By and large it works, but I keep getting deviations of 10 to 20 arcseconds. The deviations are rather systematic (in size and direction) for a given observation, but can be quite different for the next satellite or another day.

By no means I insist on doing the conversion manually- I'd be happy to use a function from pyephem / skyfield / NOVAS / pyastronomy / SOFA / anything, but I am not familiar enough with any of the packages to know how exactly it can be done.

Any ideas on how this can be done correctly?

---

Update: Thanks for the CSpice recommendation. I installed it and it works very well in replacement for my steps 1 and 2. Still I'm left with deviations, and still I wonder about the corrections I have to apply. I've got light travel time, annual aberration and satellite velocity aberration. Any ideas what else needs to be considered?

Answer 1

This question will likely get closed, but I think I am in a position to help you out, so I'll post some example code to get you going in the right direction.

Give the following a shot if you can install the SpiceyPy module (documentation here). The first function downloads all of the necessary CSPICE kernels from NASA's servers. If you end up using it, there is a special (under 'ea_latest') kernel for the latest high precision earth data, its updated somewhat regularly, so make sure to download that often (I set up an auto-download script to do this). Where ever you try to run this file from, create a sub-folder called "kernels", and these files will be downloaded automatically if you don't have them when you run the script. Thats what the first, long function does. The rest of the functions are pretty short, so you can see how easy it is to use this. I hope this helps.

import spiceypy.wrapper as spw
import numpy as np
import datetime
import os
import urllib

def get_ephem_kernels():

    # The ephemeris data for the moon and the sun is now downloaded from JPL in the form of
    # CSPICE kernels.  We use CSPICE to load all of this data into our environment so that
    # it can be called into our program whenever it is needed.

    print 'Retrieving Ephemeris Data'

    # This is the current site for the leapsecond data.  This may change in the future as more leapseconds are added.
    # check the path up to '.../lsk/' to find the most recent kernel
    ls_in = 'http://naif.jpl.nasa.gov/pub/naif/generic_kernels/lsk/naif0011.tls'
    ls_file = './kernels/ls.tls'

    # This is the ephemeris data, using DE421.  I use this because it goes continuously up to the year 2050.
    # There are other versions that we could consider using at this point, but they will only provide very minor
    # corrections to the sun and moon positions.
    de_in = 'http://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/a_old_versions/de421.bsp'
    de_file = './kernels/de.bsp'


    # These three files introduce physical constants, and additional parameters related
    # to the lunar reference frame.  At the moment, they are not used, but may be needed in order
    # to match the results of the VIIRS code.
    pc_in = 'http://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/pck00010.tpc'
    pc_file = './kernels/pck.tpc'

    mn_pa_de_in = 'http://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/moon_pa_de421_1900-2050.bpc'
    mn_pa_de_file = './kernels/mn_pa_de.bpc'

    mn_pa_in = 'http://naif.jpl.nasa.gov/pub/naif/generic_kernels/fk/satellites/moon_assoc_pa.tf'
    mn_pa_file = './kernels/mn_pa.tf'

    ea_pa_in = 'http://naif.jpl.nasa.gov/pub/naif/generic_kernels/fk/planets/'
    ea_pa_file = './kernels/ea_pa.tf'

    ea_pck_in = 'http://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/earth_070425_370426_predict.bpc'
    ea_pck_file = './kernels/ea_pck_predict.bpc'

    ea_pck_hist = 'http://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/earth_720101_070426.bpc'
    ea_pck_hist_file = './kernels/ea_pck_hist.bpc'

    ea_latest = 'http://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/earth_latest_high_prec.bpc'
    ea_latest_file = './kernels/ea_latest.bpc'

    # The first thing that we figure out is whether or not the files exist on our machine.
    # If they do, it will pass these statements. If not, it will download the necessary files.
    if os.path.isfile(ls_file) == False:
        urllib.urlretrieve(ls_in,ls_file)
    if os.path.isfile(de_file) == False:
        urllib.urlretrieve(de_in,de_file)
    if os.path.isfile(pc_file) == False:
        urllib.urlretrieve(pc_in,pc_file)
    if os.path.isfile(mn_pa_de_file) == False:
        urllib.urlretrieve(mn_pa_de_in,mn_pa_de_file)
    if os.path.isfile(mn_pa_file) == False:
        urllib.urlretrieve(mn_pa_in,mn_pa_file)
    if os.path.isfile(ea_pa_file) == False:
        urllib.urlretrieve(ea_pa_in,ea_pa_file)
    if os.path.isfile(ea_pck_file) == False:
        urllib.urlretrieve(ea_pck_in,ea_pck_file)
    if os.path.isfile(ea_latest_file) == False:
        urllib.urlretrieve(ea_latest,ea_latest_file)
    if os.path.isfile(ea_pck_hist_file) == False:
        urllib.urlretrieve(ea_pck_hist,ea_pck_hist_file)


    # This uses CSPICE to load the data so that it can be called later.  The data is loaded throughout
    # the whole environment so it can be called within any function or class object.
    spw.furnsh(ls_file)
    spw.furnsh(de_file)
    spw.furnsh(pc_file)
    spw.furnsh(mn_pa_de_file)
    spw.furnsh(mn_pa_file)
    spw.furnsh(ea_pa_file)
    spw.furnsh(ea_pck_file)
    spw.furnsh(ea_latest_file)
    spw.furnsh(ea_pck_hist_file)

def eq_to_cart(ra,dec,deg='yes'):

    '''
    This transforms equatorial coordinates to cartesian coordinates in whatever reference frame you are in.
    This just gives a unit vector however, and any range will need to be supplied by you
    '''

    if deg == 'yes':
        ra = np.deg2rad(ra)
        dec = np.deg2rad(dec)

    pv_out = spw.radrec(1.0,ra,dec)

    return pv_out

def cart_to_eq(pv,deg='yes'):

    '''
    This transforms from cartesian to equatorial in your current reference frame.  The range is in your
    current units
    '''

    eq_data = spw.recrad(pv)

    rang = eq_data[0]
    ra = eq_data[1]
    dec = eq_data[2]

    if deg == 'yes':
        ra = np.rad2deg(ra)
        dec = np.rad2deg(dec)

    return rang,ra,dec

def change_ref_frame(pv,time_ref,current='J2000',new='ITRF93'):

    '''
    This function will change the reference frame between two chosen frames using a rotation matrix
    at a given time
    '''
    new_time = spw.str2et(str(time_ref))
    pv = np.transpose(np.matrix(pv))
    rot_mat = np.matrix(spw.pxform(current,new,new_time))

    pv_rot = (rot_mat*pv).tolist()
    for i in range(len(pv_rot)):
        pv_rot[i] = pv_rot[i][0]

    return pv_rot

if __name__ == '__main__':

    get_ephem_kernels()

    ra = 50.12
    dec = -30.34

    time_ref = datetime.datetime.strptime('2015-12-15 12:00:00','%Y-%m-%d %H:%M:%S')

    pv = eq_to_cart(ra,dec)

    rang,ra,dec = cart_to_eq(pv)

    new_pv = change_ref_frame(pv,time_ref,current='J2000',new='ITRF93')

    print new_pv