Skip to content

Download NASA dataset and trace flight by Kalman Filter#684

Draft
ThreeMonth03 wants to merge 1 commit intosolvcon:masterfrom
ThreeMonth03:dev_nasa_data
Draft

Download NASA dataset and trace flight by Kalman Filter#684
ThreeMonth03 wants to merge 1 commit intosolvcon:masterfrom
ThreeMonth03:dev_nasa_data

Conversation

@ThreeMonth03
Copy link
Collaborator

Introduction

In this pull request, I download NASA fligth dataset. It includes observation data of flight, which is measured by Inertial Measurement Unit(IMU), Commercial LiDAR and ground truth.

Furthermore, I predict the position, velocity and direction of Center of Navigation(CON) with IMU data only, and the following output is the performance of baseline.

(base) threemonth03@threemonthwin:~/Downloads/modmesh$ PYTHONPATH=/home/threemonth03/Downloads/modmesh/modmesh /home/threemonth03/devenv/flavors/prime/usr/bin/python3 -m track.track_flight --no-earth-rotation
Kalman Filter Prediction Result:
final_con_pos_in_ecef = [-1330901.12447895 -5268174.02295982  3297021.64236977]
final_gt_con_pos_in_ecef = [-1387989.68147338 -5267172.39836541  3309312.46171108]
final_pos_error_ecef = [ 57088.55699443  -1001.62459441 -12290.81934131]
final_error_m = 58405.22948859141
final_quat_con_ecef = [-0.157174   -0.47091926  0.86318824  0.09185558]
final_gt_quat_con_ecef = [ 0.17784711  0.46345816 -0.86384828 -0.08569185]
final_attitude_error_deg = 2.6168453892731907

Last but not least, I also adjust the state with centrifugal force model and Coriolis force model, which can make our estimation better.

(base) threemonth03@threemonthwin:~/Downloads/modmesh$ make track
PYTHONPATH=/home/threemonth03/Downloads/modmesh/modmesh /home/threemonth03/devenv/flavors/prime/usr/bin/python3 -m track.track_flight
Kalman Filter Prediction Result:
final_con_pos_in_ecef = [-1336224.91707701 -5279162.81365069  3298133.01314272]
final_gt_con_pos_in_ecef = [-1387989.68147338 -5267172.39836541  3309312.46171108]
final_pos_error_ecef = [ 51764.76439637 -11990.41528528 -11179.44856836]
final_error_m = 54298.62762554346
final_quat_con_ecef = [-0.157174   -0.47091926  0.86318824  0.09185558]
final_gt_quat_con_ecef = [ 0.17784711  0.46345816 -0.86384828 -0.08569185]
final_attitude_error_deg = 2.6168453892731907

Implementation Method

First, I ask GPT-5.3-Codex to survey how to transform the attitude and how to interpolate the data, which is implemented in ThreeMonth03/nasa_data. After finding a decent baseline, I refactor code and rename variable bymyself. Last, I fill the Sphinx-style comment by GPT-5.3-Codex.

Copy link
Collaborator Author

@ThreeMonth03 ThreeMonth03 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@yungyuc Please review this pull request when you are available.

Comment on lines +137 to +147
.PHONY: track
NASA_DATASET_DIR := download/Flight1_Catered_Dataset-20201013
NASA_DLC_CSV := $(NASA_DATASET_DIR)/Data/dlc.csv
NASA_TRUTH_CSV := $(NASA_DATASET_DIR)/Data/truth.csv

$(NASA_DLC_CSV) $(NASA_TRUTH_CSV) &:
@mkdir -p download
PYTHONPATH=$(MODMESH_ROOT)/modmesh $(WHICH_PYTHON) -m track.dataset

track: $(NASA_DLC_CSV) $(NASA_TRUTH_CSV)
PYTHONPATH=$(MODMESH_ROOT)/modmesh $(WHICH_PYTHON) -m track.track_flight
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We can reproduce the best experiment result by the command make track.

Comment on lines +13 to +95
class KalmanFilter:
"""
Minimal linear Kalman filter implementation using NumPy.
"""

def __init__(self, dim_x, dim_z):
"""
Initialize filter dimensions and default matrices.

:param dim_x: State dimension.
:type dim_x: int
:param dim_z: Measurement dimension.
:type dim_z: int
:raises ValueError: If any dimension is non-positive.
"""
if dim_x <= 0 or dim_z <= 0:
raise ValueError("dim_x and dim_z must be positive")

self.dim_x = dim_x
self.dim_z = dim_z

self.x = np.zeros(dim_x, dtype=np.float64)
self.F = np.eye(dim_x, dtype=np.float64)
self.P = np.eye(dim_x, dtype=np.float64)
self.Q = np.eye(dim_x, dtype=np.float64)
self.H = np.zeros((dim_z, dim_x), dtype=np.float64)
self.R = np.eye(dim_z, dtype=np.float64)

def predict(self, u=None, B=None):
"""
Perform linear prediction step.

Implements:
``x = F x + B u`` and ``P = F P F^T + Q``.

:param u: Optional control input.
:type u: numpy.ndarray or None
:param B: Optional control-input matrix.
:type B: numpy.ndarray or None
:return: Predicted state and covariance.
:rtype: tuple[numpy.ndarray, numpy.ndarray]
:raises ValueError: If ``u`` is provided but ``B`` is ``None``.
"""
self.x = self.F @ self.x
if u is not None:
if B is None:
raise ValueError("B is required when u is provided")
self.x = self.x + B @ u
self.P = self.F @ self.P @ self.F.T + self.Q
return self.x, self.P

def update(self, z, H=None, R=None):
"""
Perform linear measurement update step.

:param z: Measurement vector. If ``None``, update is skipped.
:type z: numpy.ndarray or None
:param H: Optional measurement matrix overriding ``self.H``.
:type H: numpy.ndarray or None
:param R: Optional measurement covariance overriding ``self.R``.
:type R: numpy.ndarray or None
:return: Updated state and covariance.
:rtype: tuple[numpy.ndarray, numpy.ndarray]
"""
if z is None:
return self.x, self.P

z = np.asarray(z, dtype=np.float64).reshape(-1)
h = self.H if H is None else np.asarray(H, dtype=np.float64)
r = self.R if R is None else np.asarray(R, dtype=np.float64)

y = z - h @ self.x
s = h @ self.P @ h.T + r
k = self.P @ h.T @ np.linalg.inv(s)

self.x = self.x + k @ y

identity = np.eye(self.dim_x, dtype=np.float64)
i_kh = identity - k @ h
self.P = i_kh @ self.P @ i_kh.T + k @ r @ k.T

return self.x, self.P

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I make a numpy based Kalman Filter to predict data now. I will use the modmesh.KalmanFilter in the later pull request.

Comment on lines +33 to +44
@staticmethod
def gravity_ecef(pos_ecef):
"""
Compute central gravity acceleration in ECEF frame.

:param pos_ecef: Position in ECEF coordinates, shape ``(3,)``.
:type pos_ecef: numpy.ndarray
:return: Gravity acceleration in m/s^2, shape ``(3,)``.
:rtype: numpy.ndarray
"""
r = np.linalg.norm(pos_ecef)
return -Earth.MU_EARTH / (r ** 3) * pos_ecef
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The adjustment of gravity is necessary. If we don't consider gravity, the error is huge.

PYTHONPATH=/home/threemonth03/Downloads/modmesh/modmesh /home/threemonth03/devenv/flavors/prime/usr/bin/python3 -m track.track_flight
Kalman Filter Prediction Result:
final_con_pos_in_ecef = [-2112100.18722017 -7828473.02794707  4914963.02325346]
final_gt_con_pos_in_ecef = [-1387989.68147338 -5267172.39836541  3309312.46171108]
final_pos_error_ecef = [ -724110.50574679 -2561300.62958166  1605650.56154238]
final_error_m = 3108490.094146941
final_quat_con_ecef = [-0.157174   -0.47091926  0.86318824  0.09185558]
final_gt_quat_con_ecef = [ 0.17784711  0.46345816 -0.86384828 -0.08569185]
final_attitude_error_deg = 2.6168453892731907

@yungyuc yungyuc added the array Multi-dimensional array implementation label Mar 1, 2026
@yungyuc
Copy link
Member

yungyuc commented Mar 1, 2026

Good job. The code will need a very different organization than modmesh.track, which should be reserved for features more general. The code specific to the NASA problem should be outside the namespace modmesh.

We will need more discussions before moving forward the code review. I am converting the PR to draft before the conclusion is drawn. But feel free to change code any time.

@yungyuc yungyuc marked this pull request as draft March 1, 2026 14:10
dv_ecef = dcm_ecef_imu @ imu_dvel_in_imu
accel_ecef = dv_ecef / dt
accel_ecef = accel_ecef + Earth.gravity_ecef(con_pos_in_ecef)
if use_earth_rotation:
Copy link
Collaborator

@j8xixo12 j8xixo12 Mar 1, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ThreeMonth03 Could you please try to compensate the centrifugal and tangential acceleration induced by IMU lever arm effect ?

You can find the equation in my note, i think this could reduce the final position error.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

array Multi-dimensional array implementation

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants