1
$\begingroup$

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?

Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ At physics Stack Exchange, implementational questions are considered off-topic. Your question might be better suited at scicomp.stackexchange.com or the astronomy Stack Exchange (although I don't know their policies). $\endgroup$ – Martin Dec 15 '15 at 9:49
  • $\begingroup$ Have you tried using CSPICE. There are functions for coordinate system transfers that are very easy to use. Check out the module github.com/AndrewAnnex/SpiceyPy . Use the function 'pxform' (for Cartesian) to create a rotation matrix between coordinate systems at specific times. This includes ITRF93 to J2000, which I have done before. The functions 'recrad' and 'radrec' transform between Cartesian and equatorial coordinates as well. There's a lot of good stuff in there. $\endgroup$ – tmwilson26 Dec 15 '15 at 13:19
  • $\begingroup$ Would Computational Science be a better home for this question? $\endgroup$ – Qmechanic♦ Dec 15 '15 at 14:03
  • $\begingroup$ Sorry for posting an off-topic question. As you can see, I'm pretty new here and a bit lost among these ~180 stack exchange sites... The question wasn't meant to be focused on the implementation, but rather on the procedure. I tried to make that clearer with the edit. Thanks anyway for the implementational hints, CSPice seems to be great. $\endgroup$ – Daniel Hampf Dec 18 '15 at 9:54
1
$\begingroup$

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
Improve this answer
$\endgroup$
2
  • $\begingroup$ +1 (but it won't get counted until I get more score, I believe). Thanks for the code, I find it very helpful. Unfortunately, I'm still not much closer to an agreement of the data (see edit). $\endgroup$ – Daniel Hampf Dec 18 '15 at 9:57
  • $\begingroup$ @DanielHampf I'm glad you were able to use it. The things that I work on don't require that kind of accuracy so I'll have to think about what else might be done. If I think of anything I'll let you know. Are you generating satellite orbital data from TLE's? I have more functions that can do this as well. There are a number of aberration corrections that you can choose which can be found here: naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/spkezr_c.html under the inputs to this function. Maybe something in here will suit your needs. Let me know if I can help any further. $\endgroup$ – tmwilson26 Dec 18 '15 at 13:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.