SPK Checker

This notebook processes an ESOC (Mission Analysis or Flight Dynamics) OEM trajectory file into a SPICE kernel (SPK).

In [1]:
!date # Notebook was last run
Tue Oct  9 16:17:45 CEST 2018
In [2]:
%matplotlib inline

import spiceypy
import os
import matplotlib as mpl
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import numpy as np

from datetime import datetime as dt
from numpy import round
In [3]:
spiceypy.furnsh('/Users/mcosta/JUICE/kernels/mk/juice_crema_3_2_ops_local.tm')
In [4]:
input_file = 'JUI-ESOC-MOC-RP-001_(WP-578)-4/JUICE_CReMA4d0_Backup_0001.oem'
spk_file = 'JUICE_CReMA4d0_Backup_0001_man.bsp'
In [5]:
spacecraft = 'JUICE'
spacecraft_id = '-28'
frame = 'J2000'
ltcorr = 'None'

Import data

Read in the Mission Analysis orbit data using read_table function in Pandas.

The input file has no header, we'll use the regular expression for 1 or more spaces as the column separator, we'll set the index column to be the first column (=0) and we'll parse the content of that first column using the date_parser function defined above.

We'll also manually specify the column names via the names attribute to the pd.read_table call.

In [6]:
def date_parser_function(et):
    et_list.append(float(et))
    return spiceypy.et2utc(float(et),'ISOC',6)

with open('oem_table.txt', 'w') as o:
    with open(input_file,'r') as i:
        data_flag = False
        previous_tdb = ''
        for line in i:
            if 'CENTER_NAME' in line:
                center = line.split('=')[-1].strip()
            if 'META_START' in line:
                data_flag = False
            if data_flag and line.strip() != '':
                line = line.split()
                tdb = spiceypy.str2et(line[0].replace('T',' ') + ' TDB')
                if tdb != previous_tdb:
                    o.write('{} {} {} {} {} {} {} {}\n'.format(tdb,
                                                         spiceypy.bodn2c(center),
                                                         line[1],line[2],line[3],
                                                         line[4],line[5],line[6])) 
                previous_tdb = tdb
            if 'META_STOP' in line:
                data_flag = True
In [7]:
et_list = []
trajectory_data = pd.read_table(
    'oem_table.txt', skiprows=0, header=None, sep=r"\s*", 
    engine='python', index_col=0, parse_dates=True, 
    date_parser = date_parser_function,
    names=['utc','c', 'x', 'y', 'z', 'vx', 'vy', 'vz'])

trajectory_data['r'] = np.sqrt(np.square(trajectory_data[['x','y','z']]).sum(axis=1))
trajectory_data['vr'] = np.sqrt(np.square(trajectory_data[['vx','vy','vz']]).sum(axis=1))
trajectory_data['et'] = et_list

trajectory_data = trajectory_data[['et','c','x', 'y', 'z', 'r', 'vx', 'vy', 'vz', 'vr']]
/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/pandas/io/parsers.py:2227: FutureWarning: split() requires a non-empty pattern match.
  yield pat.split(line.strip())
/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/pandas/io/parsers.py:2229: FutureWarning: split() requires a non-empty pattern match.
  yield pat.split(line.strip())
In [8]:
trajectory_data.head()
Out[8]:
et c x y z r vx vy vz vr
utc
2023-08-08 23:58:50.816911 7.448112e+08 399 6245.375540 -2200.531280 291.776086 6628.136000 3.736192 10.558979 -0.337786 11.205591
2023-08-09 00:01:39.588977 7.448114e+08 399 6751.816467 -387.493981 229.518665 6766.820229 2.261523 10.854867 -0.397030 11.095056
2023-08-09 00:04:33.685624 7.448115e+08 399 7018.350692 1498.127618 156.568483 7178.171528 0.833363 10.745385 -0.437319 10.786521
2023-08-09 00:07:43.895994 7.448117e+08 399 7050.670591 3503.631110 71.064415 7873.527601 -0.435858 10.302469 -0.458213 10.321860
2023-08-09 00:11:22.404057 7.448120e+08 399 6833.849642 5681.885000 -29.744862 8887.418233 -1.480150 9.622155 -0.461776 9.746279

Determine the number of significant figures/digits for the position and velocity vector data provided.

In [9]:
xyz_signifig = len(str(trajectory_data['x'].values[0]).split('.')[1])
vxvyvz_signifig = len(str( trajectory_data['vx'].values[0]).split('.')[1])

print ("Position Vectors Significance: {}".format(xyz_signifig))
print ("Velocity Vectors Significance: {}".format(vxvyvz_signifig))
Position Vectors Significance: 5
Velocity Vectors Significance: 6

Examine sample of input data

Define start and end times for the sample of data you wish to examine. The range of times in the input data are listed below.

In [10]:
print(trajectory_data.index[0])
print(trajectory_data.index[-1])
2023-08-08 23:58:50.816911
2036-08-21 14:13:15.244477
In [11]:
#print(trajectory_data.loc['2030-10-01 15'])
#print(trajectory_data.loc['2030-10-01 17'])

#sample_start = trajectory_data.index.get_loc('2030-10-01 15:30:50.735136') 
#sample_end = trajectory_data.index.get_loc('2030-10-01 15:59:16.630431')  

sample_start = 0
sample_end = -1
#sample_end = min(1000, trajectory_data.shape[0])
In [12]:
mpl.rcParams['figure.figsize'] = (26.0, 10.0)

title = 'Input Data State Vector Components'
ax = trajectory_data[sample_start:sample_end][['x','y','z']].plot(title=title)
ax.set_ylabel("km");
plt.savefig('output/{}.png'.format('_'.join(title.split())))
In [13]:
title = 'Input Data State Vector Magnitude'
ax = trajectory_data[sample_start:sample_end][['r']].plot(title=title, marker='.')
ax.set_ylabel("km");
plt.savefig('output/{}.png'.format('_'.join(title.split())))
In [14]:
title = 'Input Data Velocity Vector Components'
ax = trajectory_data[sample_start:sample_end][['vx','vy','vz']].plot(title=title)
ax.set_ylabel("km/s");
plt.savefig('output/{}.png'.format('_'.join(title.split())))
In [15]:
title = 'Input Data Velocity Vector Magnitude'
ax = trajectory_data[sample_start:sample_end][['vr']].plot(title=title)
ax.set_ylabel("km/s");
plt.savefig('output/{}.png'.format('_'.join(title.split())))

Check the Step of the input for for interpolation analysis

In [16]:
data_sample1 = trajectory_data.head(len(trajectory_data.index)-1)
data_sample2 = trajectory_data.tail(len(trajectory_data.index)-1)

if data_sample1.shape==data_sample2.shape:
    cols = data_sample2.columns.tolist()
    diff_df = pd.DataFrame(columns=cols)
    for col in data_sample2.columns.tolist():
        diff_df[col] = data_sample2[col].values - data_sample1[col].values
    diff_df.index = data_sample2.index
else:
    print('The input and spice-sampled data were not of the same shape!')
    
diff_df.tail()
Out[16]:
et c x y z r vx vy vz vr
utc
2036-08-21 13:49:58.813058 350.734368 0 -480.113174 -122.512680 370.035407 -2.803030 0.214785 -0.144987 0.233654 0.001789
2036-08-21 13:55:48.966470 350.153412 0 -395.253308 -170.411624 443.595094 -3.492612 0.263869 -0.128877 0.188630 0.002181
2036-08-21 14:01:38.473797 349.507327 0 -295.164489 -211.590288 499.743228 -4.010357 0.302905 -0.107786 0.136273 0.002441
2036-08-21 14:07:27.209618 348.735821 0 -183.754658 -244.408404 536.222004 -4.342506 0.330265 -0.082472 0.078526 0.002565
2036-08-21 14:13:15.244477 348.034859 0 -65.501544 -267.733775 551.940800 -4.486966 0.344977 -0.053927 0.017640 0.002570
In [17]:
title = 'Delta ET Differences (Sample)'
ax = diff_df[sample_start:sample_end][['et']].plot(title=title, marker='.')
ax.set_ylabel("seconds");
plt.savefig('output/{}.png'.format('_'.join(title.split())))
In [18]:
diff_df.loc['2032-10-15 10:08:23.999992':'2032-10-15 23:12:28.800018']
Out[18]:
et c x y z r vx vy vz vr
utc
2032-10-15 13:29:46.904788 38787.231017 0 55770.430340 25013.448362 8558.543296 53094.524334 -1.380600e-02 0.000151 0.000394 -1.235935e-02
2032-10-15 13:29:47.336788 0.432000 0 0.618226 0.278660 0.095425 0.588726 0.000000e+00 0.000000 0.000000 0.000000e+00
2032-10-15 13:29:47.768788 0.432000 0 0.618227 0.278660 0.095425 0.588727 -1.000000e-06 0.000000 0.000000 -9.027741e-07
2032-10-15 13:29:48.200788 0.432000 0 0.618226 0.278660 0.095426 0.588726 0.000000e+00 0.000000 0.000000 0.000000e+00
2032-10-15 13:29:48.632788 0.432000 0 0.618226 0.278660 0.095425 0.588726 0.000000e+00 0.000000 0.000000 0.000000e+00
2032-10-15 13:29:49.064788 0.432000 0 0.618227 0.278660 0.095425 0.588727 0.000000e+00 0.000000 0.000000 0.000000e+00
2032-10-15 13:29:49.496788 0.432000 0 0.618226 0.278660 0.095425 0.588726 0.000000e+00 0.000000 0.000000 0.000000e+00
2032-10-15 13:29:49.928788 0.432000 0 0.618226 0.278660 0.095426 0.588726 0.000000e+00 0.000000 0.000000 0.000000e+00
2032-10-15 13:29:50.360788 0.432000 0 0.618226 0.278660 0.095425 0.588726 0.000000e+00 0.000000 0.000000 0.000000e+00
2032-10-15 22:19:50.846662 31800.485875 0 45335.225266 20529.216133 7036.588711 43180.546143 -1.103600e-02 0.001194 0.000837 -9.349004e-03

Re-process input

In [19]:
et_times = trajectory_data['et'].tolist()

input_trajectory = pd.DataFrame(columns=('et','c','x', 'y', 'z', 'vx', 'vy', 'vz'))

min_step = 760.3

for i, epoch in enumerate(et_times):
    if i == 0:
        step = 0
    else: 
        step += et_times[i]-et_times[i-1]
        
    if step >= min_step:
        input_trajectory.loc[i] = trajectory_data.iloc[i]
        step = 0

input_trajectory.set_index('et', inplace=True)
In [20]:
input_trajectory.tail()
Out[20]:
c x y z vx vy vz
et
1.156208e+09 503.0 2904.614803 204.468294 -1165.530481 -0.625042 0.806967 -1.451362
1.156209e+09 503.0 1789.804050 968.277157 -2402.565766 -1.433874 0.605188 -0.836294
1.156210e+09 503.0 59.671711 1404.612723 -2823.729450 -1.750954 0.198419 0.060081
1.156211e+09 503.0 -1691.026268 1362.775874 -2283.506425 -1.471783 -0.275710 0.934833
1.156212e+09 503.0 -2861.557239 858.261282 -970.132696 -0.690224 -0.657360 1.493390
In [21]:
input_trajectory[['c','x', 'y', 'z', 'vx', 'vy', 'vz']].to_csv('filtered_input.data', sep=" ", header=False)

Quick-check the SPK

The SPICE toolkit app brief can be used to see a few details of the resulting kernel, and commnt can be used to view the setup used in producing the kernel.

In [22]:
! commnt -r $spk_file
 
CREATION DATE           =         2018-05-17
ORIGIN                  =         ESOC-OPS/GFA
AUTHOR                  =         A. BOUTONNET
INTERPLANETARY EPH      =         JPL DE432
JUPITER EPH             =         IMCCE L2
 
This kernel is the end-to-end trajectory used to prepare the CReMA 4.0
 
 
 
In [23]:
! brief $spk_file -c -t -s
 
BRIEF -- Version 4.0.0, September 8, 2010 -- Toolkit Version N0066
 
 
Summary for: JUICE_CReMA4d0_Backup_0001_man.bsp
 
Bodies                         Start of Interval (ET)          End of Interval (ET)
-------                        -----------------------------   -----------------------------
-28 JUICE w.r.t. 399 EARTH     2023 AUG 08 23:59:59.999        2023 AUG 12 19:18:05.214
-28 JUICE w.r.t. 10 SUN        2023 AUG 12 19:18:05.214        2024 AUG 04 10:54:27.489
-28 JUICE w.r.t. 399 EARTH     2024 AUG 04 10:54:27.489        2024 AUG 09 03:58:14.400
-28 JUICE w.r.t. 301 MOON      2024 AUG 09 03:58:14.400        2024 AUG 09 15:59:01.156
-28 JUICE w.r.t. 399 EARTH     2024 AUG 09 15:59:01.156        2024 AUG 11 09:11:09.185
-28 JUICE w.r.t. 10 SUN        2024 AUG 11 09:11:09.185        2025 AUG 26 20:15:54.701
-28 JUICE w.r.t. 299 VENUS     2025 AUG 26 20:15:54.701        2025 AUG 28 23:13:56.921
-28 JUICE w.r.t. 10 SUN        2025 AUG 28 23:13:56.921        2026 DEC 05 09:46:49.285
-28 JUICE w.r.t. 299 VENUS     2026 DEC 05 09:46:49.285        2026 DEC 07 12:50:17.289
-28 JUICE w.r.t. 10 SUN        2026 DEC 07 12:50:17.289        2028 APR 20 02:51:49.030
-28 JUICE w.r.t. 399 EARTH     2028 APR 20 02:51:49.030        2028 APR 22 03:35:17.776
-28 JUICE w.r.t. 10 SUN        2028 APR 22 03:35:17.776        2029 DEC 27 09:41:06.404
-28 JUICE w.r.t. 399 EARTH     2029 DEC 27 09:41:06.404        2029 DEC 29 10:25:04.333
-28 JUICE w.r.t. 10 SUN        2029 DEC 29 10:25:04.333        2032 MAY 23 15:29:42.545
-28 JUICE w.r.t. 599 JUPITER   2032 MAY 23 15:29:42.545        2032 AUG 11 08:44:22.118
-28 JUICE w.r.t. 503 GANYMEDE  2032 AUG 11 08:44:22.118        2032 AUG 11 10:39:12.327
-28 JUICE w.r.t. 599 JUPITER   2032 AUG 11 10:39:12.327        2033 APR 11 19:20:37.206
-28 JUICE w.r.t. 503 GANYMEDE  2033 APR 11 19:20:37.206        2033 APR 11 21:52:19.022
-28 JUICE w.r.t. 599 JUPITER   2033 APR 11 21:52:19.022        2033 JUN 08 00:02:48.985
-28 JUICE w.r.t. 503 GANYMEDE  2033 JUN 08 00:02:48.985        2033 JUN 08 02:31:13.409
-28 JUICE w.r.t. 599 JUPITER   2033 JUN 08 02:31:13.409        2033 JUL 06 14:44:53.270
-28 JUICE w.r.t. 503 GANYMEDE  2033 JUL 06 14:44:53.270        2033 JUL 06 17:13:02.515
-28 JUICE w.r.t. 599 JUPITER   2033 JUL 06 17:13:02.515        2033 JUL 28 00:03:05.420
-28 JUICE w.r.t. 504 CALLISTO  2033 JUL 28 00:03:05.420        2033 JUL 28 04:09:41.380
-28 JUICE w.r.t. 599 JUPITER   2033 JUL 28 04:09:41.380        2033 AUG 13 16:26:40.774
-28 JUICE w.r.t. 504 CALLISTO  2033 AUG 13 16:26:40.774        2033 AUG 13 20:33:06.014
-28 JUICE w.r.t. 599 JUPITER   2033 AUG 13 20:33:06.014        2033 AUG 30 08:48:03.379
-28 JUICE w.r.t. 504 CALLISTO  2033 AUG 30 08:48:03.379        2033 AUG 30 12:54:09.988
-28 JUICE w.r.t. 599 JUPITER   2033 AUG 30 12:54:09.988        2033 SEP 16 01:11:17.292
-28 JUICE w.r.t. 504 CALLISTO  2033 SEP 16 01:11:17.292        2033 SEP 16 05:17:18.584
-28 JUICE w.r.t. 599 JUPITER   2033 SEP 16 05:17:18.584        2034 JAN 02 11:18:40.883
-28 JUICE w.r.t. 504 CALLISTO  2034 JAN 02 11:18:40.883        2034 JAN 02 15:26:57.847
-28 JUICE w.r.t. 599 JUPITER   2034 JAN 02 15:26:57.847        2034 MAR 26 21:35:40.422
-28 JUICE w.r.t. 504 CALLISTO  2034 MAR 26 21:35:40.422        2034 MAR 27 01:43:15.696
-28 JUICE w.r.t. 599 JUPITER   2034 MAR 27 01:43:15.696        2034 APR 12 14:01:15.510
-28 JUICE w.r.t. 504 CALLISTO  2034 APR 12 14:01:15.510        2034 APR 12 18:09:27.941
-28 JUICE w.r.t. 599 JUPITER   2034 APR 12 18:09:27.941        2034 APR 29 06:22:53.224
-28 JUICE w.r.t. 504 CALLISTO  2034 APR 29 06:22:53.224        2034 APR 29 10:31:04.487
-28 JUICE w.r.t. 599 JUPITER   2034 APR 29 10:31:04.487        2034 MAY 15 22:43:19.289
-28 JUICE w.r.t. 504 CALLISTO  2034 MAY 15 22:43:19.289        2034 MAY 16 02:51:43.797
-28 JUICE w.r.t. 599 JUPITER   2034 MAY 16 02:51:43.797        2034 JUN 01 14:48:57.742
-28 JUICE w.r.t. 504 CALLISTO  2034 JUN 01 14:48:57.742        2034 JUN 01 18:58:29.953
-28 JUICE w.r.t. 599 JUPITER   2034 JUN 01 18:58:29.953        2034 JUN 28 10:57:34.254
-28 JUICE w.r.t. 504 CALLISTO  2034 JUN 28 10:57:34.254        2034 JUN 28 15:07:34.484
-28 JUICE w.r.t. 599 JUPITER   2034 JUN 28 15:07:34.484        2034 JUL 15 03:21:55.574
-28 JUICE w.r.t. 504 CALLISTO  2034 JUL 15 03:21:55.574        2034 JUL 15 07:31:40.029
-28 JUICE w.r.t. 599 JUPITER   2034 JUL 15 07:31:40.029        2034 JUL 26 05:23:47.570
-28 JUICE w.r.t. 502 EUROPA    2034 JUL 26 05:23:47.570        2034 JUL 26 06:49:36.006
-28 JUICE w.r.t. 599 JUPITER   2034 JUL 26 06:49:36.006        2034 AUG 09 10:43:21.455
-28 JUICE w.r.t. 502 EUROPA    2034 AUG 09 10:43:21.455        2034 AUG 09 12:08:35.125
-28 JUICE w.r.t. 599 JUPITER   2034 AUG 09 12:08:35.125        2034 AUG 22 02:01:51.530
-28 JUICE w.r.t. 504 CALLISTO  2034 AUG 22 02:01:51.530        2034 AUG 22 05:51:32.713
-28 JUICE w.r.t. 599 JUPITER   2034 AUG 22 05:51:32.713        2034 OCT 03 03:41:26.388
-28 JUICE w.r.t. 503 GANYMEDE  2034 OCT 03 03:41:26.388        2034 OCT 03 06:00:47.172
-28 JUICE w.r.t. 599 JUPITER   2034 OCT 03 06:00:47.172        2034 OCT 26 10:24:42.024
-28 JUICE w.r.t. 504 CALLISTO  2034 OCT 26 10:24:42.024        2034 OCT 26 14:59:35.863
-28 JUICE w.r.t. 599 JUPITER   2034 OCT 26 14:59:35.863        2034 NOV 12 02:38:07.105
-28 JUICE w.r.t. 504 CALLISTO  2034 NOV 12 02:38:07.105        2034 NOV 12 07:14:04.586
-28 JUICE w.r.t. 599 JUPITER   2034 NOV 12 07:14:04.586        2035 JAN 16 11:30:36.534
-28 JUICE w.r.t. 503 GANYMEDE  2035 JAN 16 11:30:36.534        2035 JAN 16 16:17:34.178
-28 JUICE w.r.t. 599 JUPITER   2035 JAN 16 16:17:34.178        2035 JAN 29 08:08:22.469
-28 JUICE w.r.t. 503 GANYMEDE  2035 JAN 29 08:08:22.469        2035 JAN 29 12:22:15.248
-28 JUICE w.r.t. 599 JUPITER   2035 JAN 29 12:22:15.248        2035 FEB 17 11:27:59.517
-28 JUICE w.r.t. 504 CALLISTO  2035 FEB 17 11:27:59.517        2035 FEB 17 19:32:39.137
-28 JUICE w.r.t. 599 JUPITER   2035 FEB 17 19:32:39.137        2035 MAY 04 06:54:53.594
-28 JUICE w.r.t. 504 CALLISTO  2035 MAY 04 06:54:53.594        2035 MAY 04 15:02:51.431
-28 JUICE w.r.t. 599 JUPITER   2035 MAY 04 15:02:51.431        2035 NOV 09 18:09:03.474
-28 JUICE w.r.t. 503 GANYMEDE  2035 NOV 09 18:09:03.474        2036 AUG 21 14:14:24.427

Check SPK against input data

In [24]:
spiceypy.furnsh(spk_file)

Set times based on input data

These times will be in ephemerides time (et) format.

In [25]:
et_times = trajectory_data['et'].tolist()
utc_times = trajectory_data.index.tolist()
centers = trajectory_data['c'].tolist()
sample_size = len(et_times)

states_x = trajectory_data['x']
states_y = trajectory_data['y']
states_z = trajectory_data['z']

states_vx = trajectory_data['vx']
states_vy = trajectory_data['vy']
states_vz = trajectory_data['vz']

Compare sample of input (mission analysis) and output (generated SPK) data

For the times defined above build a Pandas dataframe with spacecraft states using spice.spkezr.
Set the time column as the index and map it from et to datetime objects.

In [26]:
spice_trajectory = pd.DataFrame(columns=('utc','et','x', 'y', 'z', 'vx', 'vy', 'vz'))

for i, epoch in enumerate(et_times):
    try:
        x, y, z, vx, vy, vz = spiceypy.spkezr(spacecraft, float(epoch), frame, ltcorr, spiceypy.bodc2n(centers[i]))[0]
        spice_trajectory.loc[i] =  [ utc_times[i], float(epoch), x, y, z, vx, vy, vz]
    except:
         spice_trajectory.loc[i] =  [ utc_times[i], float(epoch), states_x[i], states_y[i], states_z[i],
                                      states_vx[i], states_vy[i], states_vz[i]]
In [27]:
spice_trajectory.set_index('utc', inplace=True)

spice_trajectory['c'] = trajectory_data['c'].tolist()
spice_trajectory['r'] = np.sqrt(np.square(spice_trajectory[['x','y','z']]).sum(axis=1))
spice_trajectory['vr'] = np.sqrt(np.square(spice_trajectory[['vx','vy','vz']]).sum(axis=1))

spice_trajectory = spice_trajectory[['et','c', 'x', 'y', 'z', 'r', 'vx', 'vy', 'vz', 'vr']]

Examine the heads of the SPICE-based dataframe and read-in Mission Analysis dataframe. The latter has had its time index converted to et.

In [28]:
spice_trajectory.head()
Out[28]:
et c x y z r vx vy vz vr
utc
2023-08-08 23:58:50.816911 7.448112e+08 399 6245.376773 -2200.527808 291.775975 6628.136004 3.736189 10.558980 -0.337786 11.205591
2023-08-09 00:01:39.588977 7.448114e+08 399 6751.817335 -387.489808 229.518512 6766.820851 2.261520 10.854867 -0.397030 11.095055
2023-08-09 00:04:33.685624 7.448115e+08 399 7018.350833 1498.129428 156.568409 7178.172042 0.833362 10.745385 -0.437319 10.786521
2023-08-09 00:07:43.895994 7.448117e+08 399 7050.670562 3503.631799 71.064384 7873.527882 -0.435858 10.302469 -0.458213 10.321860
2023-08-09 00:11:22.404057 7.448120e+08 399 6833.849760 5681.884236 -29.744825 8887.417835 -1.480150 9.622155 -0.461776 9.746279
In [29]:
trajectory_data.head()
Out[29]:
et c x y z r vx vy vz vr
utc
2023-08-08 23:58:50.816911 7.448112e+08 399 6245.375540 -2200.531280 291.776086 6628.136000 3.736192 10.558979 -0.337786 11.205591
2023-08-09 00:01:39.588977 7.448114e+08 399 6751.816467 -387.493981 229.518665 6766.820229 2.261523 10.854867 -0.397030 11.095056
2023-08-09 00:04:33.685624 7.448115e+08 399 7018.350692 1498.127618 156.568483 7178.171528 0.833363 10.745385 -0.437319 10.786521
2023-08-09 00:07:43.895994 7.448117e+08 399 7050.670591 3503.631110 71.064415 7873.527601 -0.435858 10.302469 -0.458213 10.321860
2023-08-09 00:11:22.404057 7.448120e+08 399 6833.849642 5681.885000 -29.744862 8887.418233 -1.480150 9.622155 -0.461776 9.746279

Do an element-by-element comparison of the 2 dataframes if they have the same shape.

In [30]:
spice_sample = spice_trajectory.head(sample_size)
data_sample = trajectory_data.head(sample_size)
if spice_sample.shape==data_sample.shape:
    cols = spice_sample.columns.tolist()
    diff_df = pd.DataFrame(columns=cols)
    for col in spice_sample.columns.tolist():
        diff_df[col] = spice_sample[col].values - data_sample[col].values
    diff_df.index = spice_sample.index
else:
    print('The input and spice-sampled data were not of the same shape.!')
In [31]:
diff_df.shape
Out[31]:
(50492, 10)

Plot the differences.

In [32]:
title = 'State Vector Components Max Differences (Sample)'
ax = diff_df[['x','y','z']].plot(title=title)
ax.set_ylabel("km");
plt.savefig('output/{}.png'.format('_'.join(title.split())))
In [33]:
title='State Vector Magnitude Max Differences (Sample)'
ax = diff_df[['r']].plot(title=title)
ax.set_ylabel("km");
plt.savefig('output/{}.png'.format('_'.join(title.split())))
In [34]:
title='Velocity Vector Components Max. Differences (Sample)'
ax = diff_df[['vx','vy','vz']].plot(title=title)
ax.set_ylabel("km/s");
plt.savefig('output/{}.png'.format('_'.join(title.split())))

Examine input (OEM) and output (SPK) data for maximum difference

Note: this may take some time to complete. It is done row-by-row to avoid having to load large (500,000+ rows) dataframes into memory.

In [35]:
def et_to_datetime(et, scale='TDB'):
    """
    convert a SPICE ephemerides epoch (TBD seconds) to a python datetime
    object. The default time scale returned will be TDB but can be set
    to any of the accepted SPICE time scales.

    Args:
        et (float): SPICE ephemerides sceonds (TBD)
        scale (str, optional): time scale of output time (default: TDB)

    Returns:
        datetime: python datetime
    """
    t = spiceypy.timout(et, 'YYYY-MON-DD HR:MN:SC.### ::{}'.format(scale), 41)
    return dt.strptime(t, '%Y-%b-%d %H:%M:%S.%f')


max_diff = [0, 0, 0, 0, 0, 0]
max_diff_times = [0, 0, 0, 0, 0, 0]
max_diff_et = [0, 0, 0, 0, 0, 0]
max_diff_i = [0, 0, 0, 0, 0, 0]

for i, epoch in enumerate(et_times):
    try:
        state0 = spiceypy.spkezr(spacecraft, epoch, frame, ltcorr, spiceypy.bodc2n(centers[i]))[0]
    except:
        state0 = np.array((states_x[i], states_y[i], states_z[i], states_vx[i], states_vy[i], states_vz[i]))
        
    et, c, x, y, z, r, vx, vy, vz, vr = trajectory_data.values[i]
    state1 = [x, y, z, vx, vy, vz]
    gt = abs(state1-state0) > max_diff
    for j, bool in enumerate(gt):
        if bool:
            max_diff[j] = abs(state1-state0)[j]
            max_diff_times[j] = et_to_datetime(epoch)
            max_diff_et[j] = epoch
            max_diff_i[j] = i          
        else:
            pass
        
for i, val in enumerate(max_diff[:3]):
    max_diff[i]=val.round(xyz_signifig)
for i, val in enumerate(max_diff[3:]):
    max_diff[i+3]=val.round(vxvyvz_signifig)   
    
orbit_diff = pd.DataFrame.from_records(
    [max_diff_times, max_diff_et, max_diff, max_diff_i], 
    columns=('x [km]', 'y [km]', 'z [km]', 'vx [km/s]', 'vy [km/s]', 'vz [km/s]'), 
    index=['time', 'et', 'max diff', 'index'])
orbit_diff
Out[35]:
x [km] y [km] z [km] vx [km/s] vy [km/s] vz [km/s]
time 2035-05-04 03:35:03.992000 2035-02-17 08:10:53.705000 2035-02-17 08:10:53.705000 2035-05-04 03:35:03.992000 2035-02-17 08:10:53.705000 2035-02-17 08:10:53.705000
et 1.11513e+09 1.10858e+09 1.10858e+09 1.11513e+09 1.10858e+09 1.10858e+09
max diff 1133.48 723.663 342.722 0.128705 0.09824 0.046764
index 7237 7014 7014 7237 7014 7014