From 62231eef45c901573e4ee34ac23f4319a77c96d9 Mon Sep 17 00:00:00 2001 From: Christopher Rabotin Date: Sat, 5 Aug 2023 15:24:52 -0600 Subject: [PATCH 1/4] Add RIC export to traj parquet But I don't think that's the right approach --- python/nyx_space/plots/traj.py | 102 +++++++++++++++++++++++++++- src/md/trajectory/traj.rs | 30 +++++++- tests/python/test_mission_design.py | 57 +++++++++++++--- 3 files changed, 176 insertions(+), 13 deletions(-) diff --git a/python/nyx_space/plots/traj.py b/python/nyx_space/plots/traj.py index a362eceb..6b27d170 100644 --- a/python/nyx_space/plots/traj.py +++ b/python/nyx_space/plots/traj.py @@ -27,7 +27,6 @@ colors, finalize_plot, body_color, - nyx_tpl, ) @@ -205,10 +204,9 @@ def plot_orbit_elements( Args: dfs (pandas.DataFrame): The data frame containing the trajectory (or a list thereof) title (str): The title of the plot + names (List[str]): The names of each data frame html_out (str, optional): The path to save the HTML to. Defaults to None. copyright (str, optional): The copyright to display on the plot. Defaults to None. - fig (plotly.graph_objects.Figure, optional): The figure to add the trajectory to. Defaults to None. - center (str, optional): The name of the center object, e.g. `Luna` (Moon). Defaults to "Earth". show (bool, optional): Whether to show the plot. Defaults to True. If set to false, the figure will be returned. """ @@ -273,3 +271,101 @@ def plot_orbit_elements( fig.show() else: return fig + + +def plot_ric_traj( + dfs, + title, + names=[], + ric=True, + vertical=False, + html_out=None, + copyright=None, + show=True, +): + """ + Plot the trajectories in the Radial, In-track, Cross-track frame + + Args: + dfs (pandas.DataFrame): The data frame containing the trajectory (or a list thereof), must include the RIC components (`x_ric (km)`, `y_ric (km)`, `z_ric (km)`). + title (str): The title of the plot + names (List[str]): The names of each data frame + ric (bool): Set to False to plot the Cartesian X, Y, and Z components in the frame of the data, defaults to `True`. + vertical (bool): Set to True to plot X, Y, and Z of each data frame side by side (vertical layout) instead of one above another (horizontal layout). + html_out (str, optional): The path to save the HTML to. Defaults to None. + copyright (str, optional): The copyright to display on the plot. Defaults to None. + show (bool, optional): Whether to show the plot. Defaults to True. If set to false, the figure will be returned. + """ + + if not isinstance(dfs, list): + dfs = [dfs] + + for df in dfs: + pd_ok_epochs = [] + for epoch in df["Epoch:Gregorian UTC"]: + epoch = epoch.replace("UTC", "").strip() + if "." not in epoch: + epoch += ".0" + pd_ok_epochs += [epoch] + df["Epoch"] = pd.to_datetime(pd_ok_epochs) + + if not isinstance(names, list): + names = [names] + + if len(names) > 1: + assert len(names) == len(dfs), "Number of names must match number of dataframes" + + if ric: + columns = [ + "x_ric (km)", + "y_ric (km)", + "z_ric (km)", + ] + else: + columns = [ + "x (km)", + "y (km)", + "z (km)", + ] + + if vertical: + fig = make_subplots(rows=1, cols=3, subplot_titles=columns) + else: + fig = make_subplots(rows=3, cols=1, subplot_titles=columns) + + row_i = 0 + col_i = 0 + + for col in columns: + for k, df in enumerate(dfs): + try: + name = f"{names[k]} {col}" + except IndexError: + name = col + + try: + fig.add_trace( + go.Scatter(x=df["Epoch"], y=df[col], name=name), + row=row_i + 1, + col=col_i + 1, + ) + except KeyError: + raise KeyError( + f"Rebuild the trajectory and export the RIC frame: missing `{col}` for df `{names[k]}`" + ) + if vertical: + col_i += 1 + else: + row_i += 1 + + finalize_plot(fig, title, copyright=copyright) + + if html_out: + with open(html_out, "w") as f: + f.write(fig.to_html()) + print(f"Saved HTML to {html_out}") + + if show: + fig.show() + else: + return fig diff --git a/src/md/trajectory/traj.rs b/src/md/trajectory/traj.rs index a28ffe9f..bf90e648 100644 --- a/src/md/trajectory/traj.rs +++ b/src/md/trajectory/traj.rs @@ -23,7 +23,7 @@ use crate::errors::NyxError; use crate::io::watermark::pq_writer; use crate::linalg::allocator::Allocator; use crate::linalg::DefaultAllocator; -use crate::md::prelude::{GuidanceMode, StateParameter}; +use crate::md::prelude::{Frame, GuidanceMode, StateParameter}; use crate::md::EventEvaluator; use crate::time::{Duration, Epoch, TimeSeries, TimeUnits, Unit}; use arrow::array::{Array, Float64Builder, StringBuilder}; @@ -475,6 +475,17 @@ where hdrs.push(field.to_field(more_meta.clone())); } + // Add the RIC frame headers + for coord in ["x", "y", "z"] { + let mut meta = HashMap::new(); + meta.insert("unit".to_string(), "km".to_string()); + + let field = Field::new(format!("{coord}_ric (km)"), DataType::Float64, false) + .with_metadata(meta); + + hdrs.push(field); + } + if let Some(events) = events.as_ref() { for event in events { let field = Field::new(format!("{event}"), DataType::Float64, false); @@ -529,6 +540,23 @@ where record.push(Arc::new(data.finish())); } } + + // Add the RIC frame headers + for coord_no in 0..3 { + let mut data = Float64Builder::new(); + for s in &states { + // Convert to RIC. + let dcm_inertial2ric = s + .orbit() + .dcm_from_traj_frame(Frame::RIC) + .unwrap() + .transpose(); + let s_pos = dcm_inertial2ric * s.orbit().radius(); + data.append_value(s_pos[coord_no]); + } + record.push(Arc::new(data.finish())); + } + info!("Serialized {} states", states.len()); // Add all of the evaluated events diff --git a/tests/python/test_mission_design.py b/tests/python/test_mission_design.py index 6ce18727..7350b665 100644 --- a/tests/python/test_mission_design.py +++ b/tests/python/test_mission_design.py @@ -1,19 +1,22 @@ import logging -from pathlib import Path import pickle +import sys +from pathlib import Path from timeit import timeit -from nyx_space.cosmic import Spacecraft, Orbit, SrpConfig, Cosm +import pandas as pd +from nyx_space.cosmic import Cosm, Orbit, Spacecraft, SrpConfig from nyx_space.mission_design import ( - propagate, - two_body, - StateParameter, Event, SpacecraftDynamics, + StateParameter, TrajectoryLoader, + propagate, + two_body, ) -from nyx_space.time import Duration, Unit, Epoch, TimeSeries -from nyx_space.monte_carlo import generate_orbits, StateParameter +from nyx_space.monte_carlo import StateParameter, generate_orbits +from nyx_space.plots.traj import plot_ric_traj +from nyx_space.time import Duration, Epoch, TimeSeries, Unit def test_propagate(): @@ -94,10 +97,46 @@ def test_propagate(): traj = traj.resample(Unit.Second * 25.0) # Export this trajectory with additional metadata and the events + # Base path + root = Path(__file__).joinpath("../../../").resolve() + outpath = root.joinpath("output_data/") + + # Note: Python interface only supports strings for paths, not Path objects. traj.to_parquet( - "lofi_with_events.parquet", metadata={"test key": "test value"}, events=[event] + str(outpath.joinpath("./lofi_with_events.parquet")), + metadata={"dynamics": str(dynamics["lofi"])}, + events=[event], + ) + + # Propagate until the lo-fi epoch to compare + _, traj_hifi = propagate( + sc, + dynamics["hifi"], + epoch=rslt_apo.epoch, + ) + traj_hifi.to_parquet( + str(outpath.joinpath("./hifi_with_events.parquet")), + metadata={"dynamics": str(dynamics["hifi"])}, + events=[event], ) + # Plot both trajectories in RIC frame + if sys.platform != "win32": + # We must reload these as dataframes + traj_lofi_df = pd.read_parquet(outpath.joinpath("./lofi_with_events.parquet")) + traj_hifi_df = pd.read_parquet(outpath.joinpath("./hifi_with_events.parquet")) + + plot_ric_traj( + [traj_lofi_df, traj_hifi_df], + "LoFi vs HiFi", + ["LoFi", "HiFi"], + html_out=outpath.joinpath("./md_ric_lofi_hifi_diff.html"), + show=False, + ) + + # Resample the trajectory at fixed step size of 25 seconds + traj = traj.resample(Unit.Second * 25.0) + # Also export this ground track cosm = Cosm.de438() iau_earth = cosm.frame("IAU Earth") @@ -260,4 +299,4 @@ def test_merge_traj(): if __name__ == "__main__": - test_merge_traj() + test_propagate() From 55cf63ffa922a911b4e9714e5ca2d0c6df55f561 Mon Sep 17 00:00:00 2001 From: Christopher Rabotin Date: Sun, 6 Aug 2023 22:30:05 -0600 Subject: [PATCH 2/4] Continued work on RIC frame Need to refine Python API --- python/nyx_space/plots/traj.py | 23 ++-- src/md/trajectory/traj.rs | 194 ++++++++++++++++++++++++---- tests/python/test_mission_design.py | 4 +- 3 files changed, 184 insertions(+), 37 deletions(-) diff --git a/python/nyx_space/plots/traj.py b/python/nyx_space/plots/traj.py index 6b27d170..410dbd6a 100644 --- a/python/nyx_space/plots/traj.py +++ b/python/nyx_space/plots/traj.py @@ -273,10 +273,10 @@ def plot_orbit_elements( return fig -def plot_ric_traj( +def plot_traj_errors( dfs, title, - names=[], + names, ric=True, vertical=False, html_out=None, @@ -284,10 +284,10 @@ def plot_ric_traj( show=True, ): """ - Plot the trajectories in the Radial, In-track, Cross-track frame + Plot the provided trajectory data frames in the Radial, In-track, Cross-track frame of the first trajectory Args: - dfs (pandas.DataFrame): The data frame containing the trajectory (or a list thereof), must include the RIC components (`x_ric (km)`, `y_ric (km)`, `z_ric (km)`). + dfs (pandas.DataFrame): The data frame containing at least two trajectory data frames title (str): The title of the plot names (List[str]): The names of each data frame ric (bool): Set to False to plot the Cartesian X, Y, and Z components in the frame of the data, defaults to `True`. @@ -297,8 +297,16 @@ def plot_ric_traj( show (bool, optional): Whether to show the plot. Defaults to True. If set to false, the figure will be returned. """ - if not isinstance(dfs, list): - dfs = [dfs] + # TODO: This must read in an OrbitTraj or an ScTraj and a sample rate. + # For each state, build the RIC DCM of the first trajectory and convert the other trajectory states to that frame. + # And then plot that. + # You can't use the data frames as-is because the epochs may not be aligned (and that would cause weird plotting artifacts). + # You also can't convert everything to RIC beforehand because each object would be in its own RIC frame, so I and C will be zero for all. + + if isinstance(dfs, list) or len(dfs) < 2: + raise ValueError("Please provide at least two trajectories to compare") + + dfs = [dfs] for df in dfs: pd_ok_epochs = [] @@ -312,8 +320,7 @@ def plot_ric_traj( if not isinstance(names, list): names = [names] - if len(names) > 1: - assert len(names) == len(dfs), "Number of names must match number of dataframes" + assert len(names) == len(dfs), "Number of names must match number of dataframes" if ric: columns = [ diff --git a/src/md/trajectory/traj.rs b/src/md/trajectory/traj.rs index bf90e648..dc043e95 100644 --- a/src/md/trajectory/traj.rs +++ b/src/md/trajectory/traj.rs @@ -475,17 +475,6 @@ where hdrs.push(field.to_field(more_meta.clone())); } - // Add the RIC frame headers - for coord in ["x", "y", "z"] { - let mut meta = HashMap::new(); - meta.insert("unit".to_string(), "km".to_string()); - - let field = Field::new(format!("{coord}_ric (km)"), DataType::Float64, false) - .with_metadata(meta); - - hdrs.push(field); - } - if let Some(events) = events.as_ref() { for event in events { let field = Field::new(format!("{event}"), DataType::Float64, false); @@ -541,22 +530,6 @@ where } } - // Add the RIC frame headers - for coord_no in 0..3 { - let mut data = Float64Builder::new(); - for s in &states { - // Convert to RIC. - let dcm_inertial2ric = s - .orbit() - .dcm_from_traj_frame(Frame::RIC) - .unwrap() - .transpose(); - let s_pos = dcm_inertial2ric * s.orbit().radius(); - data.append_value(s_pos[coord_no]); - } - record.push(Arc::new(data.finish())); - } - info!("Serialized {} states", states.len()); // Add all of the evaluated events @@ -616,6 +589,173 @@ where Ok(traj) } + + /// Export the difference in RIC from of this trajectory compare to the "other" trajectory in parquet format. + /// + /// # Notes + /// + The RIC frame accounts for the transport theorem by performing a finite differencing of the RIC frame. + pub fn ric_diff_to_parquet>( + &self, + other: &Self, + path: P, + cfg: ExportCfg, + ) -> Result> { + let tick = Epoch::now().unwrap(); + info!("Exporting trajectory to parquet file..."); + + // Grab the path here before we move stuff. + let path_buf = cfg.actual_path(path); + + // Build the schema + let mut hdrs = vec![ + Field::new("Epoch:Gregorian UTC", DataType::Utf8, false), + Field::new("Epoch:Gregorian TAI", DataType::Utf8, false), + Field::new("Epoch:TAI (s)", DataType::Float64, false), + ]; + + // Add the RIC headers + for coord in ["x", "y", "z"] { + let mut meta = HashMap::new(); + meta.insert("unit".to_string(), "km".to_string()); + + let field = Field::new(format!("delta_{coord}_ric (km)"), DataType::Float64, false) + .with_metadata(meta); + + hdrs.push(field); + } + + let more_meta = Some(vec![( + "Frame".to_string(), + format!("{}", self.states[0].frame()), + )]); + + let mut cfg = cfg.clone(); + + let mut fields = match cfg.fields { + Some(fields) => fields, + None => S::export_params(), + }; + + // Check that we can retrieve this information + fields.retain(|param| match self.first().value(*param) { + Ok(_) => true, + Err(_) => { + warn!("Removed unavailable field `{param}` from trajectory export",); + false + } + }); + + for field in &fields { + hdrs.push(field.to_field(more_meta.clone())); + } + + // Build the schema + let schema = Arc::new(Schema::new(hdrs)); + let mut record: Vec> = Vec::new(); + + // Ensure the times match. + cfg.start_epoch = if self.first().epoch() > other.first().epoch() { + Some(self.first().epoch()) + } else { + Some(other.first().epoch()) + }; + + cfg.end_epoch = if self.last().epoch() > other.last().epoch() { + Some(other.last().epoch()) + } else { + Some(self.last().epoch()) + }; + + // Build the states iterator + let step = cfg.step.unwrap_or_else(|| 1.minutes()); + let self_states = self + .every_between(step, cfg.start_epoch.unwrap(), cfg.end_epoch.unwrap()) + .collect::>(); + + let other_states = other + .every_between(step, cfg.start_epoch.unwrap(), cfg.end_epoch.unwrap()) + .collect::>(); + + // Build all of the records + + // Epochs (both match for self and others) + let mut utc_epoch = StringBuilder::new(); + let mut tai_epoch = StringBuilder::new(); + let mut tai_s = Float64Builder::new(); + for s in &self_states { + utc_epoch.append_value(format!("{}", s.epoch())); + tai_epoch.append_value(format!("{:x}", s.epoch())); + tai_s.append_value(s.epoch().to_tai_seconds()); + } + record.push(Arc::new(utc_epoch.finish())); + record.push(Arc::new(tai_epoch.finish())); + record.push(Arc::new(tai_s.finish())); + + // Add the RIC data + for coord_no in 0..3 { + let mut data = Float64Builder::new(); + for (ii, self_state) in self_states.iter().enumerate() { + let self_orbit = self_state.orbit(); + let other_orbit = other_states[ii].orbit(); + let dcm_inertial2ric = self_orbit + .dcm_from_traj_frame(Frame::RIC) + .unwrap() + .transpose(); + + // TODO: Build the dcm time derivative + + let self_ric = dcm_inertial2ric * self_orbit.radius(); + let other_ric = dcm_inertial2ric * other_orbit.radius(); + + let delta_ric = (self_ric - other_ric)[coord_no]; + + data.append_value(delta_ric); + } + record.push(Arc::new(data.finish())); + } + + // Add all of the fields + for field in fields { + let mut data = Float64Builder::new(); + for (ii, self_state) in self_states.iter().enumerate() { + let self_val = self_state.value(field).unwrap(); + let other_val = other_states[ii].value(field).unwrap(); + data.append_value(self_val - other_val); + } + record.push(Arc::new(data.finish())); + } + + info!("Serialized {} states differences", self_states.len()); + + // Serialize all of the devices and add that to the parquet file too. + let mut metadata = HashMap::new(); + metadata.insert( + "Purpose".to_string(), + "Trajectory difference data".to_string(), + ); + if let Some(add_meta) = cfg.metadata { + for (k, v) in add_meta { + metadata.insert(k, v); + } + } + + let props = pq_writer(Some(metadata)); + + let file = File::create(&path_buf)?; + let mut writer = ArrowWriter::try_new(file, schema.clone(), props).unwrap(); + + let batch = RecordBatch::try_new(schema, record)?; + writer.write(&batch)?; + writer.close()?; + + // Return the path this was written to + let tock_time = Epoch::now().unwrap() - tick; + info!( + "Trajectory written to {} in {tock_time}", + path_buf.display() + ); + Ok(path_buf) + } } impl ops::Add for Traj diff --git a/tests/python/test_mission_design.py b/tests/python/test_mission_design.py index 7350b665..976273f9 100644 --- a/tests/python/test_mission_design.py +++ b/tests/python/test_mission_design.py @@ -15,7 +15,7 @@ two_body, ) from nyx_space.monte_carlo import StateParameter, generate_orbits -from nyx_space.plots.traj import plot_ric_traj +from nyx_space.plots.traj import plot_traj_errors from nyx_space.time import Duration, Epoch, TimeSeries, Unit @@ -126,7 +126,7 @@ def test_propagate(): traj_lofi_df = pd.read_parquet(outpath.joinpath("./lofi_with_events.parquet")) traj_hifi_df = pd.read_parquet(outpath.joinpath("./hifi_with_events.parquet")) - plot_ric_traj( + plot_traj_errors( [traj_lofi_df, traj_hifi_df], "LoFi vs HiFi", ["LoFi", "HiFi"], From 77d198c7203c59521753f321eee0be3226b8de9a Mon Sep 17 00:00:00 2001 From: Christopher Rabotin Date: Tue, 8 Aug 2023 17:56:12 -0600 Subject: [PATCH 3/4] Fix RIC and update plotting scripts + Plots now use the same color for each data frame + RIC error plots now available using the ric_diff_to_parquet function on a spacecraft trajectory + Error plots for orbital elements now available using path to parquet or OEM file --- python/nyx_space/analysis/__init__.py | 23 ++++++ python/nyx_space/analysis/traj.py | 87 ++++++++++++++++++++ python/nyx_space/plots/__init__.py | 3 +- python/nyx_space/plots/traj.py | 83 +++++++++++-------- src/md/trajectory/traj.rs | 96 ++++++++++++++++++---- src/python/mission_design/sc_trajectory.rs | 16 ++++ tests/python/test_mission_design.py | 26 +++--- tests/python/test_orbit_determination.py | 27 +++++- 8 files changed, 296 insertions(+), 65 deletions(-) create mode 100644 python/nyx_space/analysis/__init__.py create mode 100644 python/nyx_space/analysis/traj.py diff --git a/python/nyx_space/analysis/__init__.py b/python/nyx_space/analysis/__init__.py new file mode 100644 index 00000000..5333bb35 --- /dev/null +++ b/python/nyx_space/analysis/__init__.py @@ -0,0 +1,23 @@ +""" +Nyx, blazing fast astrodynamics +Copyright (C) 2023 Christopher Rabotin + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU Affero General Public License as published +by the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU Affero General Public License for more details. + +You should have received a copy of the GNU Affero General Public License +along with this program. If not, see . +""" + +from .traj import diff_traj_parquet + +__all__ = [ + "diff_traj_parquet", +] diff --git a/python/nyx_space/analysis/traj.py b/python/nyx_space/analysis/traj.py new file mode 100644 index 00000000..ee891b2f --- /dev/null +++ b/python/nyx_space/analysis/traj.py @@ -0,0 +1,87 @@ +""" +Nyx, blazing fast astrodynamics +Copyright (C) 2023 Christopher Rabotin + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU Affero General Public License as published +by the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU Affero General Public License for more details. + +You should have received a copy of the GNU Affero General Public License +along with this program. If not, see . +""" + +import pandas as pd +from nyx_space.time import Unit, TimeSeries +from nyx_space.mission_design import TrajectoryLoader +from nyx_space.monte_carlo import StateParameter + + +def diff_traj_parquet(path1: str, path2: str, step=Unit.Minute * 1) -> pd.DataFrame: + """ + Load both trajectory files from parquet, sample them at the provided step, and return an error data frame that can be plotted with the typical tools. + """ + + oe_list_str = [ + "sma", + "ecc", + "inc", + "raan", + "aop", + "ta", + "tlong", + "aol", + "x", + "y", + "z", + "vx", + "vy", + "vz", + ] + + oe_list = [StateParameter(name) for name in oe_list_str] + + def loader(path: str) -> TrajectoryLoader: + # Try to be somewhat clever + if path.lower().endswith(".oem"): + # Load as OEM and build the parquet file + return TrajectoryLoader( + str(path), "oem", path[:-3] + "parquet" + ).to_orbit_traj() + else: + return TrajectoryLoader(str(path), "parquet").to_orbit_traj() + + traj1 = loader(path1) + traj2 = loader(path2) + + err_data = {"Epoch:Gregorian UTC": []} + for oe in oe_list: + err_data[str(oe)] = [] + + # Sample to build the data. + + for epoch in TimeSeries( + traj1.first().epoch, traj1.last().epoch, step, inclusive=True + ): + try: + state1 = traj1.at(epoch) + state2 = traj2.at(epoch) + except: + # Missing data at either of these, let's ignore + pass + else: + err_data["Epoch:Gregorian UTC"] += [str(epoch)] + for oe in oe_list: + val1 = state1.value_of(oe) + val2 = state2.value_of(oe) + err_data[str(oe)] += [val1 - val2] + + # Build the data frame + err_df = pd.DataFrame(err_data, columns=list(err_data.keys())) + + return err_df diff --git a/python/nyx_space/plots/__init__.py b/python/nyx_space/plots/__init__.py index 534c2874..e379b932 100644 --- a/python/nyx_space/plots/__init__.py +++ b/python/nyx_space/plots/__init__.py @@ -18,13 +18,14 @@ from .gauss_markov import plot_gauss_markov from .od import plot_covar, plot_estimates, plot_measurements -from .traj import plot_traj, plot_ground_track +from .traj import plot_traj, plot_ground_track, plot_traj_errors __all__ = [ "plot_gauss_markov", "plot_covar", "plot_estimates", "plot_traj", + "plot_traj_errors", "plot_ground_track", "plot_measurements", ] diff --git a/python/nyx_space/plots/traj.py b/python/nyx_space/plots/traj.py index 410dbd6a..f280ede1 100644 --- a/python/nyx_space/plots/traj.py +++ b/python/nyx_space/plots/traj.py @@ -244,6 +244,8 @@ def plot_orbit_elements( row_i = 0 col_i = 0 + color_values = list(colors.values()) + for col in columns: for k, df in enumerate(dfs): try: @@ -251,8 +253,12 @@ def plot_orbit_elements( except IndexError: name = col + # Build the color for this data frame + color = color_values[k % len(color_values)] + color = f"rgb({int(color[0])}, {int(color[1])}, {int(color[2])})" + fig.add_trace( - go.Scatter(x=df["Epoch"], y=df[col], name=name), + go.Scatter(x=df["Epoch"], y=df[col], name=name, marker_color=color), row=row_i + 1, col=col_i + 1, ) @@ -276,37 +282,34 @@ def plot_orbit_elements( def plot_traj_errors( dfs, title, - names, + names=[], ric=True, vertical=False, + velocity=False, html_out=None, copyright=None, show=True, ): """ - Plot the provided trajectory data frames in the Radial, In-track, Cross-track frame of the first trajectory + Plot the trajectory error data in the Radial, In-track, Cross-track frame. Args: - dfs (pandas.DataFrame): The data frame containing at least two trajectory data frames + dfs (List[pandas.DataFrame]): The list of trajectory data frames title (str): The title of the plot - names (List[str]): The names of each data frame - ric (bool): Set to False to plot the Cartesian X, Y, and Z components in the frame of the data, defaults to `True`. - vertical (bool): Set to True to plot X, Y, and Z of each data frame side by side (vertical layout) instead of one above another (horizontal layout). - html_out (str, optional): The path to save the HTML to. Defaults to None. - copyright (str, optional): The copyright to display on the plot. Defaults to None. - show (bool, optional): Whether to show the plot. Defaults to True. If set to false, the figure will be returned. + names (List[str]): The name of each trajectory + ric (bool): Set to False to plot Cartesian errors instead of RIC. + vertical (bool): Set True to plot components vertically instead of stacked. + velocity (bool): Set True to also plot velocity errors instead of the position errors. + html_out (str, optional): Path to save HTML plot to. + copyright (str, optional): Copyright notice. + show (bool, optional): Whether to show the plot. + + Returns: + plotly.graph_objects.Figure: The Figure instance """ - # TODO: This must read in an OrbitTraj or an ScTraj and a sample rate. - # For each state, build the RIC DCM of the first trajectory and convert the other trajectory states to that frame. - # And then plot that. - # You can't use the data frames as-is because the epochs may not be aligned (and that would cause weird plotting artifacts). - # You also can't convert everything to RIC beforehand because each object would be in its own RIC frame, so I and C will be zero for all. - - if isinstance(dfs, list) or len(dfs) < 2: - raise ValueError("Please provide at least two trajectories to compare") - - dfs = [dfs] + if not isinstance(dfs, list): + dfs = [dfs] for df in dfs: pd_ok_epochs = [] @@ -320,20 +323,32 @@ def plot_traj_errors( if not isinstance(names, list): names = [names] - assert len(names) == len(dfs), "Number of names must match number of dataframes" - - if ric: - columns = [ - "x_ric (km)", - "y_ric (km)", - "z_ric (km)", - ] + if velocity: + if ric: + columns = [ + "delta_vx_ric (km/s)", + "delta_vy_ric (km/s)", + "delta_vz_ric (km/s)", + ] + else: + columns = [ + "vx (km/s)", + "vy (km/s)", + "vz (km/s)", + ] else: - columns = [ - "x (km)", - "y (km)", - "z (km)", - ] + if ric: + columns = [ + "delta_x_ric (km)", + "delta_y_ric (km)", + "delta_z_ric (km)", + ] + else: + columns = [ + "x (km)", + "y (km)", + "z (km)", + ] if vertical: fig = make_subplots(rows=1, cols=3, subplot_titles=columns) @@ -358,7 +373,7 @@ def plot_traj_errors( ) except KeyError: raise KeyError( - f"Rebuild the trajectory and export the RIC frame: missing `{col}` for df `{names[k]}`" + f"Rebuild the trajectory and export the RIC frame: missing `{col}`" ) if vertical: col_i += 1 diff --git a/src/md/trajectory/traj.rs b/src/md/trajectory/traj.rs index dc043e95..7ca06ce6 100644 --- a/src/md/trajectory/traj.rs +++ b/src/md/trajectory/traj.rs @@ -26,6 +26,7 @@ use crate::linalg::DefaultAllocator; use crate::md::prelude::{Frame, GuidanceMode, StateParameter}; use crate::md::EventEvaluator; use crate::time::{Duration, Epoch, TimeSeries, TimeUnits, Unit}; +use crate::utils::dcm_finite_differencing; use arrow::array::{Array, Float64Builder, StringBuilder}; use arrow::datatypes::{DataType, Field, Schema}; use arrow::record_batch::RecordBatch; @@ -624,6 +625,20 @@ where hdrs.push(field); } + for coord in ["x", "y", "z"] { + let mut meta = HashMap::new(); + meta.insert("unit".to_string(), "km/s".to_string()); + + let field = Field::new( + format!("delta_v{coord}_ric (km/s)"), + DataType::Float64, + false, + ) + .with_metadata(meta); + + hdrs.push(field); + } + let more_meta = Some(vec![( "Frame".to_string(), format!("{}", self.states[0].frame()), @@ -640,11 +655,14 @@ where fields.retain(|param| match self.first().value(*param) { Ok(_) => true, Err(_) => { - warn!("Removed unavailable field `{param}` from trajectory export",); + warn!("Removed unavailable field `{param}` from RIC export",); false } }); + // Disallowed fields + fields.retain(|param| param != &StateParameter::GuidanceMode); + for field in &fields { hdrs.push(field.to_field(more_meta.clone())); } @@ -672,6 +690,56 @@ where .every_between(step, cfg.start_epoch.unwrap(), cfg.end_epoch.unwrap()) .collect::>(); + let self_states_pre = self + .every_between( + step, + cfg.start_epoch.unwrap() + step - 1.seconds(), + cfg.end_epoch.unwrap(), + ) + .collect::>(); + + let self_states_post = self + .every_between( + step, + cfg.start_epoch.unwrap() + 1.seconds(), + cfg.end_epoch.unwrap(), + ) + .collect::>(); + + // Build the list of 6x6 RIC DCM + // Assuming identical rate just before and after the first DCM and for the last DCM + let mut inertial2ric_dcms = Vec::with_capacity(self_states.len()); + for ii in 0..self_states.len() { + let dcm_cur = self_states[ii] + .orbit() + .dcm_from_traj_frame(Frame::RIC) + .unwrap() + .transpose(); + + let dcm_pre = if ii == 0 { + dcm_cur + } else { + self_states_pre[ii - 1] + .orbit() + .dcm_from_traj_frame(Frame::RIC) + .unwrap() + .transpose() + }; + + let dcm_post = if ii == self_states.len() { + dcm_cur + } else { + self_states_post[ii] + .orbit() + .dcm_from_traj_frame(Frame::RIC) + .unwrap() + .transpose() + }; + + let dcm6x6 = dcm_finite_differencing(dcm_pre, dcm_cur, dcm_post); + inertial2ric_dcms.push(dcm6x6); + } + let other_states = other .every_between(step, cfg.start_epoch.unwrap(), cfg.end_epoch.unwrap()) .collect::>(); @@ -692,24 +760,22 @@ where record.push(Arc::new(tai_s.finish())); // Add the RIC data - for coord_no in 0..3 { + for coord_no in 0..6 { let mut data = Float64Builder::new(); - for (ii, self_state) in self_states.iter().enumerate() { - let self_orbit = self_state.orbit(); - let other_orbit = other_states[ii].orbit(); - let dcm_inertial2ric = self_orbit - .dcm_from_traj_frame(Frame::RIC) - .unwrap() - .transpose(); + for (ii, other_state) in other_states.iter().enumerate() { + let mut self_orbit = *self_states[ii].orbit(); + let mut other_orbit = *other_state.orbit(); - // TODO: Build the dcm time derivative + // Grab the DCM + let dcm_inertial2ric = inertial2ric_dcms[ii]; - let self_ric = dcm_inertial2ric * self_orbit.radius(); - let other_ric = dcm_inertial2ric * other_orbit.radius(); + // Rotate both into the "self" RIC + self_orbit.rotate_by(dcm_inertial2ric); + other_orbit.rotate_by(dcm_inertial2ric); - let delta_ric = (self_ric - other_ric)[coord_no]; - - data.append_value(delta_ric); + data.append_value( + (self_orbit.to_cartesian_vec() - other_orbit.to_cartesian_vec())[coord_no], + ); } record.push(Arc::new(data.finish())); } diff --git a/src/python/mission_design/sc_trajectory.rs b/src/python/mission_design/sc_trajectory.rs index 7c92da5a..f82e95ea 100644 --- a/src/python/mission_design/sc_trajectory.rs +++ b/src/python/mission_design/sc_trajectory.rs @@ -184,6 +184,22 @@ impl SpacecraftTraj { Ok(Self { inner: conv_traj }) } + fn ric_diff_to_parquet( + &self, + other: &Self, + path: String, + cfg: Option, + ) -> Result { + match self.inner.ric_diff_to_parquet( + &other.inner, + path, + cfg.unwrap_or_else(|| ExportCfg::default()), + ) { + Ok(path) => Ok(format!("{}", path.to_str().unwrap())), + Err(e) => Err(NyxError::CustomError(e.to_string())), + } + } + fn __add__(&self, rhs: &Self) -> Result { let inner = (self.inner.clone() + rhs.inner.clone())?; diff --git a/tests/python/test_mission_design.py b/tests/python/test_mission_design.py index 976273f9..10ce2f42 100644 --- a/tests/python/test_mission_design.py +++ b/tests/python/test_mission_design.py @@ -15,7 +15,7 @@ two_body, ) from nyx_space.monte_carlo import StateParameter, generate_orbits -from nyx_space.plots.traj import plot_traj_errors +from nyx_space.plots import plot_traj_errors from nyx_space.time import Duration, Epoch, TimeSeries, Unit @@ -114,26 +114,30 @@ def test_propagate(): dynamics["hifi"], epoch=rslt_apo.epoch, ) - traj_hifi.to_parquet( - str(outpath.joinpath("./hifi_with_events.parquet")), - metadata={"dynamics": str(dynamics["hifi"])}, - events=[event], - ) # Plot both trajectories in RIC frame if sys.platform != "win32": - # We must reload these as dataframes - traj_lofi_df = pd.read_parquet(outpath.joinpath("./lofi_with_events.parquet")) - traj_hifi_df = pd.read_parquet(outpath.joinpath("./hifi_with_events.parquet")) + diff_path = str(outpath.joinpath("./lofi_hifi_ric_diff.parquet")) + traj.ric_diff_to_parquet(traj_hifi, diff_path) + + traj_diff_ric = pd.read_parquet(diff_path) plot_traj_errors( - [traj_lofi_df, traj_hifi_df], + traj_diff_ric, "LoFi vs HiFi", - ["LoFi", "HiFi"], html_out=outpath.joinpath("./md_ric_lofi_hifi_diff.html"), show=False, ) + plot_traj_errors( + traj_diff_ric, + "LoFi vs HiFi", + html_out=outpath.joinpath("./md_ric_lofi_hifi_diff_velocity.html"), + vertical=True, + velocity=True, + show=False, + ) + # Resample the trajectory at fixed step size of 25 seconds traj = traj.resample(Unit.Second * 25.0) diff --git a/tests/python/test_orbit_determination.py b/tests/python/test_orbit_determination.py index 0a053f4e..120459f3 100644 --- a/tests/python/test_orbit_determination.py +++ b/tests/python/test_orbit_determination.py @@ -24,6 +24,7 @@ plot_residual_histogram, ) from nyx_space.plots.traj import plot_orbit_elements +from nyx_space.analysis import diff_traj_parquet def test_filter_arc(): @@ -172,6 +173,15 @@ def test_filter_arc(): show=False, ) + # Plot the errors between the nominal and estimation + err_df = diff_traj_parquet(traj_file, rslt_path) + plot_orbit_elements( + err_df, + "Error in orbital elements", + html_out=str(outpath.joinpath("./od_vs_ref_error_elements.html")), + show=False, + ) + def test_one_way_msr(): """ @@ -202,9 +212,17 @@ def test_one_way_msr(): gs = devices[1] print(f"Using {gs}") cosm = Cosm.de438() - data = {"epoch": [], "range (km)": [], "doppler (km/s)": [], "azimuth (deg)": [], "elevation (deg)": []} + data = { + "epoch": [], + "range (km)": [], + "doppler (km/s)": [], + "azimuth (deg)": [], + "elevation (deg)": [], + } # Start by building a time series - ts = TimeSeries(traj.first().epoch, traj.last().epoch, step=Unit.Minute*30, inclusive=True) + ts = TimeSeries( + traj.first().epoch, traj.last().epoch, step=Unit.Minute * 30, inclusive=True + ) # And iterate over it for epoch in ts: orbit = traj.at(epoch).orbit @@ -235,12 +253,13 @@ def test_one_way_msr(): assert abs(doppler_km_s - -0.2498238312640348) < 0.1 # Azimuth and elevation - + az_deg, el_deg = devices[0].compute_azimuth_elevation(end_sc.orbit, cosm) assert abs(az_deg - 128.66181520071825) < 1e-10 assert abs(el_deg - 27.904687635388676) < 1e-10 + def test_pure_prediction(): # Initialize logging FORMAT = "%(levelname)s %(name)s %(asctime)-15s %(filename)s:%(lineno)d %(message)s" @@ -304,4 +323,4 @@ def test_pure_prediction(): if __name__ == "__main__": - test_one_way_msr() + test_filter_arc() From 4d5268b275aba9b64d61fbfb75584c4472236f90 Mon Sep 17 00:00:00 2001 From: Christopher Rabotin Date: Wed, 9 Aug 2023 22:18:36 -0600 Subject: [PATCH 4/4] Clippy --- src/md/trajectory/traj.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/md/trajectory/traj.rs b/src/md/trajectory/traj.rs index 7ca06ce6..f0140467 100644 --- a/src/md/trajectory/traj.rs +++ b/src/md/trajectory/traj.rs @@ -644,7 +644,7 @@ where format!("{}", self.states[0].frame()), )]); - let mut cfg = cfg.clone(); + let mut cfg = cfg; let mut fields = match cfg.fields { Some(fields) => fields,