Initial scaffolding build for geophysical research

.gitignore

@@ -43,6 +43,7 @@ batch_process.nu
 # Data output
 examples/geonav-sim/ukf
 examples/geonav-sim/pf
+data/output/*
 
 # Ignore ruff
 .ruff_cache/
@@ -51,5 +52,3 @@ papers/anomaly_ukf/data/combined_statistics_error.csv
 # mdBook output
 book/book/
 
-# Generated performance plots
-data/output/*.png

analysis/src/analysis/__init__.py

@@ -10,7 +10,7 @@
 from pathlib import Path
 
 from analysis.preprocess import preprocess_data
-from analysis.plotting import plot_performance
+from analysis.plotting import plot_performance, plot_relative_performance
 from haversine import haversine_vector, Unit
 
 __version__ = "0.1.0"
@@ -85,12 +85,46 @@ def main() -> None:
         default="data/output"
     )
 
+    geoperformance = command.add_parser(
+        "geoperformance", help="Generate geophysical performance plots."
+    )
+    geoperformance.add_argument(
+        "-p",
+        "--processed",
+        type=str,
+        help="Input directory containing the processed navigation result CSV files.",
+    )
+    geoperformance.add_argument(
+        "-r",
+        "--reference",
+        type=str,
+        help="Input directory containing the reference GPS CSV files.",
+        default="data/input"
+    )
+    geoperformance.add_argument(
+        "-d",
+        "--degraded",
+        type=str,
+        help="Input directory containing the degraded navigation result CSV files.",
+    )
+    geoperformance.add_argument(
+        "-o",
+        "--output",
+        type=str,
+        help="Output directory for the geophysical performance plots.",
+        default="data/output"
+    )
+
     args = parser.parse_args()
 
     if args.command == "preprocess":
         preprocess_data(args)
     elif args.command == "performance":
         performance_analysis(args)
+    elif args.command == "geoperformance":
+        geophysical_performance_analysis(args)
+    else:
+        parser.print_help()
 
 
 def performance_analysis(args):
@@ -158,6 +192,52 @@ def performance_analysis(args):
     summary_df.to_csv(summary_file)
     print("Performance analysis completed.")
 
+def geophysical_performance_analysis(args):
+    """Generate geophysical performance plots."""
+    input_dir = args.processed
+    print(f"Generating geophysical performance plots from data in: {input_dir}")
+
+    datasets = list(Path(input_dir).glob("*.csv"))
+    print(f"Found {len(datasets)} datasets to process.")
+
+    print(f"Comparing to reference data in: {args.reference}")
+    references = list(Path(args.reference).glob("*.csv"))
+    print(f"Found {len(references)} reference datasets.")
+
+    print(f"Comparing to degraded data in: {args.degraded}")
+    degradeds = list(Path(args.degraded).glob("*.csv"))
+    print(f"Found {len(degradeds)} degraded datasets.")
+
+    output_path = Path(args.output)
+    output_path.mkdir(parents=True, exist_ok=True)
+    print(f"Saving geophysical performance plots to: {args.output}")
+
+    reference_path = Path(args.reference)
+    degraded_path = Path(args.degraded)
+
+    for dataset in datasets:
+        geo = read_csv(dataset, parse_dates=True, index_col=0)
+        try:
+            reference_file = reference_path / dataset.name
+            nav = read_csv(reference_file, parse_dates=True, index_col=0)
+        except FileNotFoundError:
+            print(f"Reference file for {dataset.name} not found in {reference_path}. Skipping.")
+            continue
+        try:
+            degraded_file = degraded_path / dataset.name
+            degraded_nav = read_csv(degraded_file, parse_dates=True, index_col=0)
+        except FileNotFoundError:
+            print(f"Degraded file for {dataset.name} not found in {degraded_path}. Skipping.")
+            continue
+        output_plot = output_path / f"{dataset.stem}_geophysical_performance.png"
+        print(f"Processing dataset {dataset} ({len(nav)}) with reference {reference_file.name} ({len(nav)}) and degraded {degraded_file.name} ({len(degraded_nav)})")
+        try:
+            plot_relative_performance(geo, degraded_nav, nav, output_plot)
+        except Exception as e:
+            print(f"Error plotting geophysical performance for {dataset.name}, possible dimension mismatch or missing data: {e}")
+            continue
+
+    print("Geophysical performance analysis completed.")
 
 if __name__ == "__main__":
     main()

analysis/src/analysis/plotting.py

@@ -138,7 +138,7 @@ def plot_performance(nav: DataFrame, gps: DataFrame, output_path: Path | str):
         linestyle="--",
     )
     ax.set_xlim(left=0)
-    ax.set_ylim((0, 50))
+    #ax.set_ylim((0, 50))
     ax.set_xlabel("Time (s)")
     ax.set_ylabel("2D Haversine Error (m)")
     ax.set_title(
@@ -150,6 +150,82 @@ def plot_performance(nav: DataFrame, gps: DataFrame, output_path: Path | str):
     plt.close(fig)
 
 
+def plot_relative_performance(
+    geo: DataFrame,
+    deg: DataFrame, # degraded/baseline navigation
+    nav: DataFrame, # GPS-aided truth
+    output_path: Path | str,
+):
+    """Plot relative performance: geophysical-aided nav w.r.t. degraded nav.
+
+    This mirrors the geonav-sim notebook visualization by plotting the
+    difference of the cumulative root-mean-squared error series:
+
+        root_mean_geo_cum_error - root_mean_deg_cum_error
+
+    Positive values are filled red, negative values green.
+    """
+    output_path = Path(output_path)
+
+    distance_traveled = haversine_vector(
+        nav[["latitude", "longitude"]].to_numpy()[:-1, :],
+        nav[["latitude", "longitude"]].to_numpy()[1:, :],
+        Unit.METERS,
+    )
+    distance_traveled = np.hstack(([0], distance_traveled))
+    distance_traveled = np.nancumsum(distance_traveled)
+
+    geo_error = haversine_vector(
+        geo[["latitude", "longitude"]].to_numpy(dtype=np.float64, copy=False),
+        nav[["latitude", "longitude"]].to_numpy(),
+        Unit.METERS,
+    )
+
+    deg_error = haversine_vector(
+        deg[["latitude", "longitude"]].to_numpy(),
+        nav[["latitude", "longitude"]].to_numpy(),
+        Unit.METERS,
+    )
+
+    time = (nav.index - nav.index[0]).total_seconds() / 3600
+
+    err_diff = geo_error - deg_error
+    geo_rmse = np.sqrt(np.nanmean(geo_error**2))
+    deg_rmse = np.sqrt(np.nanmean(deg_error**2))
+    # General errors
+    fig, ax = plt.subplots(1, 1, figsize=(24, 6), layout="tight")
+    ax.plot(time, err_diff, label="Error Difference", color="black", linewidth=1)
+    ax.fill_between(
+        time,
+        err_diff,
+        0,
+        where=err_diff > 0,
+        alpha=0.4,
+        color="red",
+        label="Geo INS > Degraded INS",
+    )
+    ax.fill_between(
+        time,
+        err_diff,
+        0,
+        where=err_diff < 0,
+        alpha=0.4,
+        color="green",
+        label="Degraded INS > Geo INS",
+    )
+    ax.set_xlim(left=0)
+    # ax[0].set_ylim((-0.1, 0.1))
+    ax.set_xlabel("Time (h)")
+    ax.set_ylabel("Distance (m)")
+    ax.set_title(
+        f"Geophysical Navigation Performance | RMSE difference: {geo_rmse - deg_rmse:0.2f}",
+        fontsize=16,
+    )
+    ax.grid()
+    ax.legend()
+    fig.savefig(output_path, dpi=300)
+    plt.close(fig)
+
 
 def plot_street_map(
     nav: DataFrame,

conf/ekf_both.toml

@@ -0,0 +1,43 @@
+# EKF with both gravity and magnetic anomaly measurements
+input = "data/input"
+output = "data/output/ekf/both"
+mode = "closed-loop"
+seed = 42
+parallel = true
+
+[logging]
+level = "info"
+file = "log/ekf_both.log"
+
+[closed_loop]
+filter = "ekf"
+
+[geophysical]
+# Gravity anomaly measurements
+gravity_resolution = "one_minute"
+gravity_bias = 0.0
+gravity_noise_std = 100.0
+
+# Magnetic anomaly measurements
+magnetic_resolution = "two_minutes"
+magnetic_bias = 0.0
+magnetic_noise_std = 150.0
+
+# Measurement frequency (seconds)
+geo_frequency_s = 1.0
+
+[gnss_degradation]
+seed = 42
+
+[gnss_degradation.scheduler]
+kind = "fixed_interval"
+interval_s = 60.0
+phase_s = 0.0
+
+[gnss_degradation.fault]
+kind = "degraded"
+rho_pos = 0.99
+sigma_pos_m = 15.0
+rho_vel = 0.95
+sigma_vel_mps = 5.0
+r_scale = 15.0

conf/ekf_degraded.toml

@@ -4,6 +4,10 @@ mode = "closed-loop"
 seed = 42
 generate_plot = true
 
+[logging]
+level = "info"
+file = "log/ekf_degraded.log"
+
 [closed_loop]
 filter = "ekf"
 

conf/ekf_grav.toml

@@ -0,0 +1,37 @@
+# EKF with gravity anomaly measurements only
+input = "data/input"
+output = "data/output/ekf/grav"
+mode = "closed-loop"
+seed = 42
+parallel = true
+
+[logging]
+level = "info"
+file = "log/ekf_grav.log"
+
+[closed_loop]
+filter = "ekf"
+
+[geophysical]
+# Gravity anomaly measurements
+gravity_resolution = "one_minute"
+gravity_bias = 0.0
+gravity_noise_std = 100.0
+# gravity_map_file will be auto-detected if not specified
+
+# Magnetic anomaly disabled for this config
+# magnetic_resolution = ""
+# magnetic_bias = 0.0
+# magnetic_noise_std = 150.0
+
+# Measurement frequency (seconds)
+geo_frequency_s = 1.0
+
+[gnss_degradation]
+seed = 42
+
+[gnss_degradation.scheduler]
+kind = "pass_through"
+
+[gnss_degradation.fault]
+kind = "none"

conf/ekf_mag.toml

@@ -0,0 +1,37 @@
+# EKF with magnetic anomaly measurements only
+input = "data/input"
+output = "data/output/ekf/mag"
+mode = "closed-loop"
+seed = 42
+parallel = true
+
+[logging]
+level = "info"
+file = "log/ekf_mag.log"
+
+[closed_loop]
+filter = "ekf"
+
+[geophysical]
+# Magnetic anomaly measurements
+magnetic_resolution = "two_minutes"
+magnetic_bias = 0.0
+magnetic_noise_std = 150.0
+# magnetic_map_file will be auto-detected if not specified
+
+# Gravity anomaly disabled for this config
+# gravity_resolution = ""
+# gravity_bias = 0.0
+# gravity_noise_std = 100.0
+
+# Measurement frequency (seconds)
+geo_frequency_s = 1.0
+
+[gnss_degradation]
+seed = 42
+
+[gnss_degradation.scheduler]
+kind = "pass_through"
+
+[gnss_degradation.fault]
+kind = "none"

conf/ekf_truth.toml

@@ -4,6 +4,10 @@ mode = "closed-loop"
 seed = 42
 generate_plot = true
 
+[logging]
+level = "info"
+file = "log/ekf_truth.log"
+
 [closed_loop]
 filter = "ekf"