Using the MagNav Package with Real SGL Flight Data

julia> using Pluto
julia> Pluto.run() # select & open notebook

Import packages and DataFrames

begin
    cd(@__DIR__)
    # uncomment line below to use local MagNav.jl (downloaded folder)
    # using Pkg; Pkg.activate("../"); Pkg.instantiate()
    using MagNav
    using CSV, DataFrames
    using Plots: plot, plot!
    using Random: seed!
    using Statistics: mean, median, std
    seed!(33); # for reproducibility
    include("dataframes_setup.jl"); # setup DataFrames
end;

Load flight and map data

begin
    flight = :Flt1006 # select flight, full list in df_flight
    xyz    = get_XYZ(flight,df_flight) # load flight data
end;

df_flight

typeof(xyz)

MagNav.XYZ20{Int64, Float64}

fieldnames(MagNav.XYZ20)

(:info, :traj, :ins, :flux_a, :flux_b, :flux_c, :flux_d, :flight, :line, :year, :doy, :utm_x, :utm_y, :utm_z, :msl, :baro, :diurnal, :igrf, :mag_1_c, :mag_1_lag, :mag_1_dc, :mag_1_igrf, :mag_1_uc, :mag_2_uc, :mag_3_uc, :mag_4_uc, :mag_5_uc, :mag_6_uc, :ogs_mag, :ogs_alt, :ins_wander, :ins_roll, :ins_pitch, :ins_yaw, :roll_rate, :pitch_rate, :yaw_rate, :ins_acc_x, :ins_acc_y, :ins_acc_z, :lgtl_acc, :ltrl_acc, :nrml_acc, :pitot_p, :static_p, :total_p, :cur_com_1, :cur_ac_hi, :cur_ac_lo, :cur_tank, :cur_flap, :cur_strb, :cur_srvo_o, :cur_srvo_m, :cur_srvo_i, :cur_heat, :cur_acpwr, :cur_outpwr, :cur_bat_1, :cur_bat_2, :vol_acpwr, :vol_outpwr, :vol_bat_1, :vol_bat_2, :vol_res_p, :vol_res_n, :vol_back_p, :vol_back_n, :vol_gyro_1, :vol_gyro_2, :vol_acc_p, :vol_acc_n, :vol_block, :vol_back, :vol_srvo, :vol_cabt, :vol_fan, :aux_1, :aux_2, :aux_3)

begin
    map_name = :Eastern_395 # select map, full list in df_map
    df_options = df_nav[(df_nav.flight   .== flight  ) .&
                        (df_nav.map_name .== map_name),:]
end

begin
    line = df_options.line[1] # select flight line (row) from df_options
    ind  = get_ind(xyz,line,df_nav) # get Boolean indices
    # ind = get_ind(xyz;lines=[line]) # alternative
end;

begin
    TL_i   = 6 # select first calibration box of 1006.04
    TL_ind = get_ind(xyz;tt_lim=[df_cal.t_start[TL_i],df_cal.t_end[TL_i]])
end;

Baseline plots of scalar and vector (fluxgate) magnetometer data

begin
    show_plot    = false
    save_plot    = false
    detrend_data = true
end;

b1 = plot_mag(xyz;ind,show_plot,save_plot,detrend_data,
              use_mags = [:mag_1_uc,:mag_3_uc,:mag_4_uc,:mag_5_uc])

b2 = plot_mag(xyz;ind,show_plot,save_plot,detrend_data,
              use_mags = [:flux_d]) # try changing to :flux_a, :flux_b, :flux_c

b3 = plot_mag(xyz;ind,show_plot,save_plot,
              use_mags = [:comp_mags])

b4 = plot_mag_c(xyz,xyz;ind,show_plot,save_plot,detrend_data,
                ind_comp  = TL_ind,
                use_mags  = [:mag_1_uc],
                use_vec   = :flux_d,
                plot_diff = true)

b5 = plot_frequency(xyz;ind,show_plot,save_plot,detrend_data,
                    field     = :mag_1_uc,
                    freq_type = :PSD)

b6 = plot_frequency(xyz;ind,show_plot,save_plot,detrend_data,
                    field     = :mag_1_uc,
                    freq_type = :spec)

Tolles-Lawson calibration and compensation

begin
    λ       = 0.025   # ridge parameter for ridge regression
    use_vec = :flux_d # selected vector (flux) magnetometer
    flux    = getfield(xyz,use_vec) # load Flux D data
    TL_d_1  = create_TL_coef(flux,xyz.mag_1_uc,TL_ind;λ=λ) # Flux D & Mag 1
    TL_d_2  = create_TL_coef(flux,xyz.mag_2_uc,TL_ind;λ=λ) # Flux D & Mag 2
    TL_d_3  = create_TL_coef(flux,xyz.mag_3_uc,TL_ind;λ=λ) # Flux D & Mag 3
    TL_d_4  = create_TL_coef(flux,xyz.mag_4_uc,TL_ind;λ=λ) # Flux D & Mag 4
    TL_d_5  = create_TL_coef(flux,xyz.mag_5_uc,TL_ind;λ=λ) # Flux D & Mag 5
end;

begin
    A = create_TL_A(flux,ind)
    mag_1_sgl = xyz.mag_1_c[ind]
    mag_1_uc  = xyz.mag_1_uc[ind]
    mag_2_uc  = xyz.mag_2_uc[ind]
    mag_3_uc  = xyz.mag_3_uc[ind]
    mag_4_uc  = xyz.mag_4_uc[ind]
    mag_5_uc  = xyz.mag_5_uc[ind]
end;

begin
    mag_1_c = mag_1_uc - detrend(A*TL_d_1;mean_only=true)
    mag_2_c = mag_2_uc - detrend(A*TL_d_2;mean_only=true)
    mag_3_c = mag_3_uc - detrend(A*TL_d_3;mean_only=true)
    mag_4_c = mag_4_uc - detrend(A*TL_d_4;mean_only=true)
    mag_5_c = mag_5_uc - detrend(A*TL_d_5;mean_only=true)
end;

begin
    traj = get_traj(xyz,ind) # trajectory (GPS) struct
    ins  = get_ins( xyz,ind;N_zero_ll=1) # INS struct, "zero" lat/lon to match first `traj` data point
    mapS = get_map(map_name,df_map) # load map data
    # get map values & map interpolation function
    (map_val,itp_mapS) = get_map_val(mapS,traj;return_itp=true)
    map_val += (xyz.diurnal + xyz.igrf)[ind] # add in diurnal & core (IGRF)
end;

begin
    println("Mag 1: ",round(std(map_val-mag_1_c),digits=2))
    println("Mag 2: ",round(std(map_val-mag_2_c),digits=2))
    println("Mag 3: ",round(std(map_val-mag_3_c),digits=2))
    println("Mag 4: ",round(std(map_val-mag_4_c),digits=2))
    println("Mag 5: ",round(std(map_val-mag_5_c),digits=2))
end

Navigation

(P0,Qd,R) = create_model(traj.dt,traj.lat[1];
                         init_pos_sigma = 0.1,
                         init_alt_sigma = 1.0,
                         init_vel_sigma = 1.0,
                         meas_var       = 5^2, # increase if mag_use is bad
                         fogm_sigma     = 3,
                         fogm_tau       = 180);

begin
    mag_use = mag_1_c # selected magnetometer, modify & see what happens
    (crlb_out,ins_out,filt_out) = run_filt(traj,ins,mag_use,itp_mapS,:ekf;
                                           P0,Qd,R,core=true)
end;

begin
    t0  = traj.tt[1]/60    # [min]
    tt  = traj.tt/60 .- t0 # [min]
    dpi = 200
end;

begin
    p1 = plot(xlab="time [min]",ylab="magnetic field [nT]",legend=:topleft,dpi=dpi)
    plot!(p1,tt,detrend(mag_1_uc ),lab="SGL raw Mag 1" ,color=:cyan,lw=2)
    plot!(p1,tt,detrend(mag_1_sgl),lab="SGL comp Mag 1",color=:blue,lw=2)
    plot!(p1,tt,detrend(mag_1_c  ),lab="MIT comp Mag 1",color=:red ,lw=2,ls=:dash)
    # plot!(p1,tt,detrend(mag_2_c  ),lab="MIT comp Mag 2",color=:purple) # bad
    plot!(p1,tt,detrend(mag_3_c  ),lab="MIT comp Mag 3",color=:green)
    plot!(p1,tt,detrend(mag_4_c  ),lab="MIT comp Mag 4",color=:black)
    plot!(p1,tt,detrend(mag_5_c  ),lab="MIT comp Mag 5",color=:orange)
    # png(p1,"comp_prof_1") # to save figure
end

begin
    p2 = plot_map(mapS;map_color=:gray) # map background
    plot_filt!(p2,traj,ins,filt_out;show_plot=false) # overlay GPS, INS, & filter
    plot!(p2,legend=:topleft) # move as needed
end

begin
    p3 = plot(xlab="time [min]",ylab="error [m]",legend=:topleft,dpi=dpi)
    plot!(p3,tt,ins_out.n_err,lab="northing")
    plot!(p3,tt,ins_out.e_err,lab="easting")
end

(p4,p5) = plot_filt_err(traj,filt_out,crlb_out;show_plot,save_plot);

p4

p5

p6 = plot_mag_map(traj,mag_use,itp_mapS)

In-flight events

begin
    p7 = plot(xlab="time [min]",ylab="magnetic field [nT]",dpi=dpi)
    plot!(p7,tt,mag_1_uc,lab="mag_1_uc")
    # plot!(p7,tt,mag_2_uc,lab="mag_2_uc") # bad
    plot!(p7,tt,mag_3_uc,lab="mag_3_uc")
    plot!(p7,tt,mag_4_uc,lab="mag_4_uc")
    plot!(p7,tt,mag_5_uc,lab="mag_5_uc")
    plot_events!(p7,flight,df_event;t0=t0,t_units=:min)
    display(p7)
end