From 9c22b15dce05b167895d31d281999679f1e56e88 Mon Sep 17 00:00:00 2001 From: Christopher Rabotin Date: Sun, 8 Dec 2024 00:40:29 -0700 Subject: [PATCH] Reenable obstruction in tracking data --- examples/04_lro_od/dsn-network.yaml | 6 +-- examples/04_lro_od/main.rs | 2 +- examples/04_lro_od/plot_od_rslt.py | 66 +++++++---------------------- examples/04_lro_od/requirements.txt | 38 +++++++++++++++++ src/od/ground_station/trk_device.rs | 10 ++++- 5 files changed, 66 insertions(+), 56 deletions(-) create mode 100644 examples/04_lro_od/requirements.txt diff --git a/examples/04_lro_od/dsn-network.yaml b/examples/04_lro_od/dsn-network.yaml index 5f6bc446..4484ed62 100644 --- a/examples/04_lro_od/dsn-network.yaml +++ b/examples/04_lro_od/dsn-network.yaml @@ -15,7 +15,7 @@ DSS-65 Madrid: white_noise: mean: 0.0 sigma: 50.0e-6 # 5 cm/s - light_time_correction: false + light_time_correction: true latitude_deg: 40.427222 longitude_deg: 4.250556 height_km: 0.834939 @@ -43,7 +43,7 @@ DSS-34 Canberra: white_noise: mean: 0.0 sigma: 50.0e-6 # 5 cm/s - light_time_correction: false + light_time_correction: true measurement_types: - range_km - doppler_km_s @@ -68,7 +68,7 @@ DSS-13 Goldstone: white_noise: mean: 0.0 sigma: 50.0e-6 # 5 cm/s - light_time_correction: false + light_time_correction: true measurement_types: - range_km - doppler_km_s diff --git a/examples/04_lro_od/main.rs b/examples/04_lro_od/main.rs index c948a0e3..6a92101b 100644 --- a/examples/04_lro_od/main.rs +++ b/examples/04_lro_od/main.rs @@ -253,7 +253,7 @@ fn main() -> Result<(), Box> { // Increase the initial covariance to account for larger deviation. initial_estimate, // Until https://github.com/nyx-space/nyx/issues/351, we need to specify the SNC in the acceleration of the Moon J2000 frame. - SNC3::from_diagonal(10 * Unit::Minute, &[1e-11, 1e-11, 1e-11]), + SNC3::from_diagonal(10 * Unit::Minute, &[1e-12, 1e-12, 1e-12]), ); // We'll set up the OD process to reject measurements whose residuals are mover than 4 sigmas away from what we expect. diff --git a/examples/04_lro_od/plot_od_rslt.py b/examples/04_lro_od/plot_od_rslt.py index eeafa63d..176d53f0 100644 --- a/examples/04_lro_od/plot_od_rslt.py +++ b/examples/04_lro_od/plot_od_rslt.py @@ -4,9 +4,13 @@ import numpy as np import plotly.graph_objects as go import plotly.express as px +import click -if __name__ == "__main__": - df = pl.read_parquet("output_data/ekf_rng_dpl_az_el_odp.parquet") + +@click.command +@click.option("-p", "--path", type=str, default="./04_lro_od_results.parquet") +def main(path: str): + df = pl.read_parquet(path) df = ( df.with_columns(pl.col("Epoch (UTC)").str.to_datetime("%Y-%m-%dT%H:%M:%S%.f")) @@ -49,68 +53,25 @@ # Add scatter points fig_qq.add_trace( go.Scatter( - x=qq[0][0], - y=qq[0][1], - mode='markers', - name='Sample Data', - marker=dict(color='blue') + x=qq[0][0], y=qq[0][1], mode="markers", name="Residuals ratios (QQ)", marker=dict(color="blue") ) ) # Add the theoretical line fig_qq.add_trace( - go.Scatter( - x=x_qq, - y=y_qq, - mode='lines', - name='Theoretical Normal', - line=dict(color='red') - ) + go.Scatter(x=x_qq, y=y_qq, mode="lines", name="Theoretical Normal", line=dict(color="red")) ) # Update layout fig_qq.update_layout( - title='Normal Q-Q Plot', - xaxis_title='Theoretical Quantiles', - yaxis_title='Sample Quantiles', + title="Normal Q-Q Plot", + xaxis_title="Theoretical Quantiles", + yaxis_title="Sample Quantiles", ) # Show QQ plot fig_qq.show() - # Create histogram with normal distribution overlay - hist_fig = px.histogram( - df, - x="Residual ratio", - color="Tracker", - marginal="rug", - hover_data=df.columns, - ) - - # Calculate normal distribution parameters - mean = residual_ratio.mean() - std = residual_ratio.std() - x_range = np.linspace(residual_ratio.min(), residual_ratio.max(), 100) - y_normal = stats.norm.pdf(x_range, mean, std) - - # Scale the normal distribution to match histogram height - max_hist_height = 100 - scaling_factor = max_hist_height / max(y_normal) - y_normal_scaled = y_normal * scaling_factor - - # Add normal distribution curve - hist_fig.add_trace( - go.Scatter( - x=x_range, - y=y_normal_scaled, - name='Normal Distribution', - line=dict(color='red', width=2), - ) - ) - - # Show histogram with normal overlay - hist_fig.show() - px.histogram( df, x="Residual ratio", @@ -163,7 +124,6 @@ y=["Sigma Vx (RIC) (km/s)", "Sigma Vy (RIC) (km/s)", "Sigma Vz (RIC) (km/s)"], ).show() - raise AssertionError("stop") # Load the RIC diff. for fname, errname in [ ("04_lro_od_truth_error", "OD vs Flown"), @@ -219,3 +179,7 @@ ], title=f"Velocity error with {errname} ({fname})", ).show() + + +if __name__ == "__main__": + main() diff --git a/examples/04_lro_od/requirements.txt b/examples/04_lro_od/requirements.txt new file mode 100644 index 00000000..c7556494 --- /dev/null +++ b/examples/04_lro_od/requirements.txt @@ -0,0 +1,38 @@ +asttokens==2.4.1 +attrs==23.2.0 +click==8.1.7 +decorator==5.1.1 +executing==2.0.1 +iniconfig==2.0.0 +ipython==8.25.0 +jedi==0.19.1 +matplotlib-inline==0.1.7 +maturin==1.6.0 +numpy==1.26.4 +packaging==24.1 +pandas==2.1.4 +parso==0.8.4 +pexpect==4.9.0 +pip==24.0 +plotly==5.16.1 +pluggy==1.5.0 +polars==0.20.31 +prompt-toolkit==3.0.47 +ptyprocess==0.7.0 +pure-eval==0.2.2 +pyarrow==13.0.0 +pygments==2.18.0 +pytest==7.2.2 +python-dateutil==2.9.0.post0 +python-slugify==8.0.4 +pytz==2024.1 +ruff==0.5.0 +scipy==1.11.4 +six==1.16.0 +stack-data==0.6.3 +tenacity==8.3.0 +text-unidecode==1.3 +traitlets==5.14.3 +typing-extensions==4.12.2 +tzdata==2024.1 +wcwidth==0.2.13 diff --git a/src/od/ground_station/trk_device.rs b/src/od/ground_station/trk_device.rs index b2b31e05..9d973614 100644 --- a/src/od/ground_station/trk_device.rs +++ b/src/od/ground_station/trk_device.rs @@ -84,6 +84,14 @@ impl TrackingDevice for GroundStation { self.name, self.elevation_mask_deg, aer_t0.elevation_deg, aer_t1.elevation_deg ); return Ok(None); + } else if aer_t0.is_obstructed() || aer_t1.is_obstructed() { + debug!( + "{} obstruction at t0={}, t1={} -- no measurement", + self.name, + aer_t0.is_obstructed(), + aer_t1.is_obstructed() + ); + return Ok(None); } // Noises are computed at the midpoint of the integration time. @@ -128,7 +136,7 @@ impl TrackingDevice for GroundStation { action: "computing AER", })?; - if aer.elevation_deg >= self.elevation_mask_deg { + if aer.elevation_deg >= self.elevation_mask_deg && !aer.is_obstructed() { // Only update the noises if the measurement is valid. let noises = self.noises(rx.orbit.epoch, rng)?;