Section # 1 - First, we store the data into an array after initial cleaning (here, it is done manually and weird looking data points are removed by visual inspection). Further, the data is plotted to visualize the given information.
In [3]:
import json
import numpy as np
import matplotlib.pyplot as plt
# load ephemeris data from the gps.json file
with open('gps.json') as f: data_f = json.load(f)
# Matrix [7,33]
data = np.array([[sub['gnss.sample.t_j2000'] for sub in data_f], [sub['gnss.sample.rx_ecef'] for sub in data_f], [sub['gnss.sample.ry_ecef'] for sub in data_f],[sub['gnss.sample.rz_ecef'] for sub in data_f], \
[sub['gnss.sample.vx_ecef'] for sub in data_f],[sub['gnss.sample.vy_ecef'] for sub in data_f],[sub['gnss.sample.vz_ecef'] for sub in data_f]])
time_UTC = np.array([sub['t'] for sub in data_f])
# Rows: [t,x,y,z,vx,vy,vz]'
# Columns: [1,2,...., 33] datapoints
data_modified = np.delete(data,16,1) # data point deleted because all measurements are 0
time_UTC = np.delete(time_UTC,16)
time = data_modified[0,:]
x_gps, y_gps, z_gps = data_modified[1,:], data_modified[2,:], data_modified[3,:]
vx_gps, vy_gps, vz_gps = data_modified[4,:], data_modified[5,:], data_modified[6,:]
# Plot raw GPS position measurements
fig, axs = plt.subplots(3)
plt.suptitle("GPS Position measurements")
axs[0].plot(time,x_gps,'r*')
axs[0].set(ylabel='x coord.')
axs[1].plot(time,y_gps,'b*')
axs[1].set(ylabel='y coord.')
axs[2].plot(time,z_gps,'g*')
axs[2].set(xlabel='time',ylabel='z coord.')
for ax in axs.flat:ax.label_outer()
# Plot raw GPS velocity measurements
fig, axs = plt.subplots(3)
plt.suptitle("GPS Velocity measurements")
axs[0].plot(time,vx_gps,'r*')
axs[0].set(ylabel='x velocity')
axs[1].plot(time,vy_gps,'b*')
axs[1].set(ylabel='y velocity')
axs[2].plot(time,vz_gps,'g*')
axs[2].set(xlabel='time',ylabel='velocity')
for ax in axs.flat:ax.label_outer()
In [4]:
# loading orekit library for orbital analysis
import orekit
orekit.initVM()
from orekit.pyhelpers import setup_orekit_curdir
setup_orekit_curdir()
Section # 2 - Save Measurements in different frames
In [5]:
from org.hipparchus.geometry.euclidean.threed import Vector3D
from org.orekit.frames import FramesFactory
from org.orekit.utils import IERSConventions, Constants, PVCoordinates, PVCoordinatesProvider, AbsolutePVCoordinates
from orekit.pyhelpers import datetime_to_absolutedate
from dateutil import parser
ITRF = FramesFactory.getITRF(IERSConventions.IERS_2010, True)
GCRF = FramesFactory.getGCRF()
EME = FramesFactory.getEME2000()
pos_gcrf = []
vel_gcrf = []
pos_itrf = []
vel_itrf = []
pos_eme = []
vel_eme = []
for i in range(32):
pos_ind = Vector3D(data_modified[1:4,i].tolist())
pos_itrf.append(pos_ind)
vel_ind = Vector3D(data_modified[4:7,i].tolist())
vel_itrf.append(vel_ind)
pv_itrf = PVCoordinates(pos_ind, vel_ind)
time_index = datetime_to_absolutedate(parser.parse(time_UTC[i]))
pv_gcrf = ITRF.getTransformTo(GCRF, time_index).transformPVCoordinates(pv_itrf)
pos_gcrf.append(pv_gcrf.getPosition())
vel_gcrf.append(pv_gcrf.getVelocity())
pv_eme = ITRF.getTransformTo(EME, time_index).transformPVCoordinates(pv_itrf)
pos_eme.append(pv_eme.getPosition())
vel_eme.append(pv_eme.getVelocity())
Transfer saved measurements to pandas dataframe type
In [6]:
import pandas as pd
position_rawGPS = pd.DataFrame(columns=['x', 'y', 'z', 'vx', 'vy', 'vz'])
position_GCRF = pd.DataFrame(columns=['x', 'y', 'z', 'vx', 'vy', 'vz'])
position_EME = pd.DataFrame(columns=['x', 'y', 'z', 'vx', 'vy', 'vz'])
for i in range(32):
time = data_modified[0,i]
pos = data_modified[1:4,i]
vel = data_modified[4:7,i]
position_rawGPS.loc[parser.parse(time_UTC[i])] = np.concatenate((pos,vel),axis=0)
pos_gcrf_ind = np.array(pos_gcrf[i].toArray())
vel_gcrf_ind = np.array(vel_gcrf[i].toArray())
position_GCRF.loc[parser.parse(time_UTC[i])] = np.concatenate((pos_gcrf_ind,vel_gcrf_ind),axis=0)
pos_eme_ind = np.array(pos_eme[i].toArray())
vel_eme_ind = np.array(vel_eme[i].toArray())
position_EME.loc[parser.parse(time_UTC[i])] = np.concatenate((pos_eme_ind,vel_eme_ind),axis=0)
print('Raw GPS Measurements')
display(position_rawGPS)
print('GPS Measurements in GCRF Frame')
display(position_GCRF)
print('GPS Measurements in EME2000 Frame')
display(position_EME)
Section # 3 - Make an initial guess based on three evenly spaced measurements by the Gibbs IOD method
In [7]:
from org.orekit.estimation.iod import IodGibbs
from org.orekit.utils import Constants as orekit_constants
from org.orekit.time import TimeScalesFactory
utc = TimeScalesFactory.getUTC()
# select evenly spaced three datapoints from the measurements to get an initial orbit starting point
measurement_1 = pos_gcrf[0]
date_1 = datetime_to_absolutedate(parser.parse(time_UTC[0]))
measurement_2 = pos_gcrf[10]
date_2 = datetime_to_absolutedate(parser.parse(time_UTC[10]))
measurement_3 = pos_gcrf[20]
date_3 = datetime_to_absolutedate(parser.parse(time_UTC[20]))
iod_gibbs = IodGibbs(orekit_constants.EIGEN5C_EARTH_MU)
orbit_first_guess = iod_gibbs.estimate(GCRF, measurement_1, date_1, measurement_2, date_2, measurement_3, date_3)
from org.orekit.propagation.analytical import KeplerianPropagator
kepler_propagator_iod = KeplerianPropagator(orbit_first_guess)
display(orbit_first_guess)
Section # 4 - Assumptions about noise in the measurements and hyperparameters for the estimations
In [8]:
sigma_position = 100e3 # Noise (in terms of standard deviation of gaussian distribution) of input position data in meters
sigma_velocity = 100.0 # Noise of input velocity data in meters per second
# Estimator parameters
estimator_position_scale = 1.0 # m
estimator_convergence_thres = 1e-2
estimator_max_iterations = 35
estimator_max_evaluations = 35
# Orbit propagator parameters
prop_min_step = 0.001 # s
prop_max_step = 300.0 # s
prop_position_error = 10.0 # m
Section # 5 - Initialize the estimator with the above hyperparemeters
In [9]:
from org.orekit.propagation.conversion import DormandPrince853IntegratorBuilder
integratorBuilder = DormandPrince853IntegratorBuilder(prop_min_step, prop_max_step, prop_position_error)
from org.orekit.propagation.conversion import NumericalPropagatorBuilder
from org.orekit.orbits import PositionAngle
propagatorBuilder = NumericalPropagatorBuilder(orbit_first_guess,
integratorBuilder, PositionAngle.TRUE, estimator_position_scale)
from org.hipparchus.linear import QRDecomposer
matrix_decomposer = QRDecomposer(1e-11)
from org.hipparchus.optim.nonlinear.vector.leastsquares import GaussNewtonOptimizer
optimizer = GaussNewtonOptimizer(matrix_decomposer, False)
from org.orekit.estimation.leastsquares import BatchLSEstimator
estimator = BatchLSEstimator(optimizer, propagatorBuilder)
estimator.setParametersConvergenceThreshold(estimator_convergence_thres)
estimator.setMaxIterations(estimator_max_iterations)
estimator.setMaxEvaluations(estimator_max_evaluations)
Section # 6 - Feeding only position measurements into estimator since we don't trust velocity
In [10]:
from orekit.pyhelpers import datetime_to_absolutedate
from org.orekit.estimation.measurements import Position, ObservableSatellite
observableSatellite = ObservableSatellite(0) # Propagator index = 0
for i in range(len(pos_gcrf)):
orekit_position = Position(
datetime_to_absolutedate(parser.parse(time_UTC[i])),
pos_gcrf[i],
sigma_position,
1.0, # Base weight
observableSatellite
)
estimator.addMeasurement(orekit_position)
Obtain the estimated propagated satellite state from the estimator
In [11]:
estimatedPropagatorArray = estimator.estimate()
Section # 7 - Propagate the estimator between specified bounds
In [12]:
dt = 60.0 # propagate at 1 minute intervals
date_start = datetime_to_absolutedate(parser.parse('2020-09-01T00:00:00.000000'))
date_end = datetime_to_absolutedate(parser.parse('2020-09-03T00:00:00.000000'))
# First propagating in ephemeris mode
estimatedPropagator = estimatedPropagatorArray[0]
estimatedInitialState = estimatedPropagator.getInitialState()
display(estimatedInitialState.getOrbit())
estimatedPropagator.resetInitialState(estimatedInitialState)
estimatedPropagator.setEphemerisMode()
estimatedPropagator.propagate(date_start, date_end)
bounded_propagator = estimatedPropagator.getGeneratedEphemeris()
Section # 8 - Save results for data analysis and visualization
In [14]:
from orekit.pyhelpers import absolutedate_to_datetime
import pandas as pd
pv_bls_df = pd.DataFrame(columns=['x', 'y', 'z', 'vx', 'vy', 'vz'])
pv_iod_df = pd.DataFrame(columns=['x', 'y', 'z', 'vx', 'vy', 'vz'])
pv_bls_list = []
pv_iod_list = []
date_list = []
date_current = date_start
while date_current.compareTo(date_end) <= 0:
datetime_current = absolutedate_to_datetime(date_current)
spacecraftState = bounded_propagator.propagate(date_current)
pv_bls = spacecraftState.getPVCoordinates(GCRF)
pos_bls = np.array(pv_bls.getPosition().toArray())
vel_bls = np.array(pv_bls.getVelocity().toArray())
pv_bls_df.loc[datetime_current] = np.concatenate((pos_bls, vel_bls))
pv_bls_list.append(np.concatenate((pos_bls, vel_bls),axis=None))
pv_iod_ind = kepler_propagator_iod.getPVCoordinates(date_current, GCRF)
pos_iod = np.array(pv_iod_ind.getPosition().toArray())
vel_iod = np.array(pv_iod_ind.getVelocity().toArray())
pv_iod_df.loc[datetime_current] = np.concatenate((pos_iod,vel_iod))
pv_iod_list.append(np.concatenate((pos_iod, vel_iod),axis=None))
date_list.append(date_current)
date_current = date_current.shiftedBy(dt)
Section # 9 - Tranform estimated trajectory to EME2000 and ITRF frame
In [15]:
pos_bls_array = np.zeros([len(date_list),3]) # shape - [len(date_list) X (x,y,z)]
vel_bls_array = np.zeros_like(pos_bls_array)
pos_iod_array = np.zeros([len(date_list),3])
vel_iod_array = np.zeros_like(pos_iod_array)
for j in range(len(date_list)):
pos_bls_array[j,:] = pv_bls_list[j][0:3]
vel_bls_array[j,:] = pv_bls_list[j][3:6]
pos_iod_array[j,:] = pv_iod_list[j][0:3]
vel_iod_array[j,:] = pv_iod_list[j][3:6]
In [16]:
ITRF = FramesFactory.getITRF(IERSConventions.IERS_2010, True)
GCRF = FramesFactory.getGCRF()
EME = FramesFactory.getEME2000()
est_pos_itrf = []
est_vel_itrf = []
est_pos_eme = []
est_vel_eme = []
for i in range(len(date_list)):
pos_ind = Vector3D(pos_bls_array[i,:].tolist())
vel_ind = Vector3D(vel_bls_array[i,:].tolist())
pv_gcrf = PVCoordinates(pos_ind, vel_ind)
pv_itrf = GCRF.getTransformTo(ITRF, date_list[i]).transformPVCoordinates(pv_gcrf)
pv_eme = GCRF.getTransformTo(EME, date_list[i]).transformPVCoordinates(pv_gcrf)
est_pos_itrf.append(pv_itrf.getPosition())
est_vel_itrf.append(pv_itrf.getVelocity())
est_pos_eme.append(pv_eme.getPosition())
est_vel_eme.append(pv_eme.getVelocity())
In [17]:
estimate_pv_EME2000 = pd.DataFrame(columns=['x', 'y', 'z', 'vx', 'vy', 'vz'])
estimate_pv_ITRF = pd.DataFrame(columns=['x', 'y', 'z', 'vx', 'vy', 'vz'])
for i in range(len(date_list)):
pos_itrf_ind = np.array(est_pos_itrf[i].toArray())
vel_itrf_ind = np.array(est_vel_itrf[i].toArray())
estimate_pv_ITRF.loc[absolutedate_to_datetime(date_list[i])] = np.concatenate((pos_itrf_ind,vel_itrf_ind),axis=0)
pos_eme_ind = np.array(est_pos_eme[i].toArray())
vel_eme_ind = np.array(est_pos_eme[i].toArray())
estimate_pv_EME2000.loc[absolutedate_to_datetime(date_list[i])] = np.concatenate((pos_eme_ind,vel_eme_ind),axis=0)
print('Estimate state in ITRF')
display(estimate_pv_ITRF)
print('Estimate state in EME2000')
display(estimate_pv_EME2000)
Section # 10 - PLOTTING ESTIMATE AND MEASUREMENTS IN EME2000 FRAME
In [18]:
import plotly.graph_objects as go
# position_EME
fig = go.Figure(data=[
go.Scatter(x=position_EME.index, y=position_EME['x'], name='X coord. GPS - measurements', mode='markers'),
go.Scatter(x=estimate_pv_EME2000.index, y=estimate_pv_EME2000['x'], mode='lines', name='Batch least squares X-approximation (EME2000 frame)')])
fig.show()
fig = go.Figure(data=[
go.Scatter(x=position_EME.index, y=position_EME['y'], name='Y coord. GPS - measurements', mode='markers'),
go.Scatter(x=estimate_pv_EME2000.index, y=estimate_pv_EME2000['y'], mode='lines', name='Batch least squares Y-approximation (EME2000 frame)')])
fig.show()
fig = go.Figure(data=[
go.Scatter(x=position_EME.index, y=position_EME['z'], name='Z coord. GPS - measurements', mode='markers'),
go.Scatter(x=estimate_pv_EME2000.index, y=estimate_pv_EME2000['z'], mode='lines', name='Batch least squares Z-approximation (EME2000 frame)')])
fig.show()
PLOTTING ESTIMATE AND MEASUREMENTS IN ITRF FRAME (original measurements in this frame)
In [19]:
import plotly.graph_objects as go
fig = go.Figure(data=[
go.Scatter(x=position_rawGPS.index, y=position_rawGPS['x'], name='X coord. GPS - measurements', mode='markers'),
go.Scatter(x=estimate_pv_ITRF.index, y=estimate_pv_ITRF['x'], mode='lines', name='Batch least squares X-approximation (ITRF frame)')])
fig.show()
fig = go.Figure(data=[
go.Scatter(x=position_rawGPS.index, y=position_rawGPS['y'], name='Y coord. GPS - measurements', mode='markers'),
go.Scatter(x=estimate_pv_ITRF.index, y=estimate_pv_ITRF['y'], mode='lines', name='Batch least squares Y-approximation (ITRF frame)')])
fig.show()
fig = go.Figure(data=[
go.Scatter(x=position_rawGPS.index, y=position_rawGPS['z'], name='Z coord. GPS - measurements', mode='markers'),
go.Scatter(x=estimate_pv_ITRF.index, y=estimate_pv_ITRF['z'], mode='lines', name='Batch least squares Z-approximation (ITRF frame)')])
fig.show()
3D Orbit Visulaization of the Determined Orbit
In [20]:
# extracting points used for finding initial orbit for plotting
import math
i = math.ceil(len(position_GCRF.index) / 3)-1
points_for_iod = position_GCRF.iloc[::i, :]
points_for_iod = points_for_iod[:-1]
display(points_for_iod)
In [26]:
measurement_1 = pos_gcrf[0]
date_1 = datetime_to_absolutedate(parser.parse(time_UTC[0]))
measurement_2 = pos_gcrf[10]
date_2 = datetime_to_absolutedate(parser.parse(time_UTC[10]))
measurement_3 = pos_gcrf[20]
date_3 = datetime_to_absolutedate(parser.parse(time_UTC[20]))
fig_data = data=[go.Scatter3d(x=pv_iod_df['x'], y=pv_iod_df['y'], z=pv_iod_df['z'], mode='lines', name='IOD solution'),
go.Scatter3d(x=pv_bls_df['x'], y=pv_bls_df['y'], z=pv_bls_df['z'], mode='lines', name='Batch least squares solution'),
go.Scatter3d(x=position_GCRF['x'], y=position_GCRF['y'], z=position_GCRF['z'], mode='markers', name='Measurements in GCRF Frame'),
go.Scatter3d(x=points_for_iod['x'], y=points_for_iod['y'], z=points_for_iod['z'], mode='markers', name='Measurements used for IOD')]
scene=dict(aspectmode='data',
)
layout = dict(
scene=scene)
fig = go.Figure(data=fig_data,
layout=layout)
fig.show()
In [28]:
residuals_df = pd.DataFrame(columns=['bls_minus_measurements_norm', 'iod_minus_measurements_norm'])
for timestamp, pv_gcrf in position_GCRF.iterrows():
date_current = datetime_to_absolutedate(timestamp)
pv_bls = bounded_propagator.getPVCoordinates(date_current, GCRF)
pos_bls = np.array(pv_bls.getPosition().toArray())
pv_iod = kepler_propagator_iod.getPVCoordinates(date_current, GCRF)
pos_iod = np.array(pv_iod.getPosition().toArray())
pv_measurements = position_GCRF.loc[timestamp]
pos_measurements = pv_measurements[['x', 'y', 'z']]
bls_minus_measurements = np.linalg.norm(pos_bls - pos_measurements)
iod_minus_measurements = np.linalg.norm(pos_iod - pos_measurements)
residuals_df.loc[timestamp] = [
np.linalg.norm(pos_bls - pos_measurements),
np.linalg.norm(pos_iod - pos_measurements)
]
display(residuals_df)
In [30]:
fig = go.Figure(data=[
go.Scatter(x=residuals_df.index, y=residuals_df['bls_minus_measurements_norm'],
name='BLS - measurements',
mode='markers+lines'),
go.Scatter(x=residuals_df.index, y=residuals_df['iod_minus_measurements_norm'],
name='IOD - measurements',
mode='markers+lines')
])
fig.show()
Section # 11 - Export ephemeris data to OEM file in the CCSDS format
In [151]:
from datetime import datetime
# export ephemris data to OEM in the CCSDS format
def export_OEM(estimated_pv_df, date_list)->str:
"""Export ephemerides in CCSDS OEM format.
Parameters
----------
estimated_pv_df : Estimation results dataframe to export
date_list : dates/times for each estimate datapoint
Returns
-------
Ephemerides in OEM format.
"""
# frame = Frame.ICRF if (cfg.prop_inertial_frame == Frame.GCRF) else cfg.prop_inertial_frame
oem_header = f"""CCSDS_OEM_VERS = 2.0
CREATION_DATE = {datetime.utcnow().strftime("%Y-%m-%dT%H:%M:%S")}
ORIGINATOR = PRIT CHOVATIYA
META_START
OBJECT_NAME = DOVE by PLANET
OBJECT_ID = 0f4e
CENTER_NAME = EARTH
REF_FRAME = EME2000
TIME_SYSTEM = UTC
START_TIME = {date_list[0]}
STOP_TIME = {date_list[-1]}
META_STOP
"""
eph_data = []
for i in range(len(date_list)):
utc = date_list[i]
pv_dat = estimated_pv_df.loc[absolutedate_to_datetime(utc)].astype(str)
eph_data.append(f"{utc}\t\t{pv_dat[0][0:15]}\t\t{pv_dat[1][0:15]}\t\t{pv_dat[2][0:15]}\t\t{pv_dat[3][0:15]}\t\t{pv_dat[4][0:15]}\t\t{pv_dat[5][0:15]}")
oem_data = oem_header + "\n".join(eph_data)
return(oem_data)
In [142]:
with open("oem.txt", "w") as fp:
fp.write(export_OEM(estimate_pv_EME2000, date_list)) # writing estimates ephem in EME2000 frame
fp.close()
In [ ]: