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
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.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']]
In [8]:
trajectory_data.head()
Out[8]:
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))
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])
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]:
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]:
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]:
In [21]:
input_trajectory[['c','x', 'y', 'z', 'vx', 'vy', 'vz']].to_csv('filtered_input.data', sep=" ", header=False)
In [22]:
! commnt -r $spk_file
In [23]:
! brief $spk_file -c -t -s
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]:
In [29]:
trajectory_data.head()
Out[29]:
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]:
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]: