Magnetic Anomaly Maps

This file is best viewed in a Pluto notebook. To run it this way, from the MagNav.jl directory, do:

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

This is a reactive notebook, so feel free to change any parameters.

Import packages and DataFrames

The DataFrames listed below provide useful information about the flight data collected by Sander Geophysics Ltd. (SGL) & magnetic anomaly maps.

Dataframe	Description
`df_map`	map files relevant for SGL flights
`df_cal`	SGL calibration flight lines
`df_flight`	SGL flight files
`df_all`	all flight lines
`df_nav`	all navigation-capable flight lines
`df_event`	pilot-recorded in-flight events

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 Perth map

This is the Perth map (at 800 m) as provided by Sander Geophysics Ltd.

begin
    map_gxf = MagNav.ottawa_area_maps_gxf()*"/Perth_Mag.gxf"
    p_mapS_800 = map_gxf2h5(map_gxf,800)
    p1 = plot_map(p_mapS_800)
end

Display Perth map in Google Earth

Display the Perth map in Google Earth by uncommenting below to generate a KMZ file, then open in Google Earth.

# map2kmz(p_mapS_800,"Perth")

Overlaid on the Perth map are most of the flight lines used to generate the map, which were collected during Flight 1004 (see readme) & Flight 1005 (see readme). The full list of SGL flights is in df_flight.

begin
    xyz_1004 = get_XYZ(:Flt1004,df_flight;silent=true) # load flight data
    xyz_1005 = get_XYZ(:Flt1005,df_flight;silent=true) # load flight data
    lines = [4019,4018,4017,4016,4015,4012,4011,4010,4009,4008,4007,4004,
             4003,4002,4001,421,419,417,415,413,411,409,408,407,405,403,401]
    for line in lines
        xyz = line in xyz_1004.line ? xyz_1004 : xyz_1005
        ind = get_ind(xyz,line,df_nav) # get Boolean indices
        plot_path!(p1,xyz.traj,ind;show_plot=false,path_color=:black)
    end
    p1
end

Load previously processed maps

Eastern Ontario & NAMAD maps have been processed & saved previously. The full list of maps is in df_map.

begin # e_mask contains Boolean indices for "real" map data (not filled-in)
    e_mapS_395 = get_map(:Eastern_395,df_map) # load map data
    e_mask     = e_mapS_395.mask
end;

begin
    xx_lim = extrema(e_mapS_395.xx) .+ (-0.01,0.01)
    yy_lim = extrema(e_mapS_395.yy) .+ (-0.01,0.01)
    n_mapS_395 = upward_fft(map_trim(get_map(MagNav.namad),
                            xx_lim=xx_lim,yy_lim=yy_lim),e_mapS_395.alt)
end;

Plot all Ottawa area maps

3 maps are overlayed with 3 different color schemes.

begin
    clims = (-500,500)
    dpi   = 50
    p2 = plot_map(n_mapS_395;clims=clims,dpi=dpi,legend=false)    # 395 m
    plot_map!(p2,e_mapS_395;clims=clims,dpi=dpi,map_color=:magma) # 395 m
    plot_map!(p2,p_mapS_800;clims=clims,dpi=dpi,map_color=:gray)  # 800 m
    p2
end

Plot Eastern Ontario altitude map CDF

Most of the map data was collected at altitudes between 215 & 395 m.

begin # this map contains an additional drape (altitude) map
    e_mapS_drp = get_map(:Eastern_drape,df_map) # load map data
    e_alts  = e_mapS_drp.alt[e_mask]
    alt_avg = round(Int,mean(e_alts))
    alt_med = round(Int,median(e_alts))
    alt_val = 200:500
    alt_cdf = [sum(e_alts .< a) for a in alt_val] / sum(e_mask)
    p3 = plot(alt_val,alt_cdf,xlab="altitude [m]",ylab="fraction [-]",
              title="altitude map CDF",lab=false,dpi=200)
end

Plot Eastern Ontario altitude map

Minimal areas have map data collected at 395 m or higher (colored spots), so a level map can be generated at 395 m using almost entirely upward continuation & minimal downward continuation.

begin
    p4 = plot_map(e_mapS_395;dpi=dpi,legend=false,map_color=:gray)
    e_mapS_395_ = deepcopy(e_mapS_395)
    e_mapS_395_.map[e_mapS_drp.alt .<= 395] .= 0
    plot_map!(p4,e_mapS_395_;dpi=dpi)
    p4
end

Plot expanded Eastern Ontario map

The original map area is show with a black outline. During upward (or downward) continuation, the map is temporarily expanded with "wrapped" edges for a more accurate result.

begin
    using MagNav: get_step
    (map_map,px,py) = MagNav.map_expand(e_mapS_395.map,200)
    px_end    = size(map_map,2) - length(e_mapS_395.xx) - px
    py_end    = size(map_map,1) - length(e_mapS_395.yy) - py
    dx        = get_step(e_mapS_395.xx)
    dy        = get_step(e_mapS_395.yy)
    xx        = [(e_mapS_395.xx[1]-dx*px):dx:(e_mapS_395.xx[end]+dx*px_end);]
    yy        = [(e_mapS_395.yy[1]-dy*py):dy:(e_mapS_395.yy[end]+dy*py_end);]
    (lat,lon) = map_border(e_mapS_395;sort_border=true) # get map border
    p5 = plot_map(map_map,xx,yy;dpi=dpi)
    plot_path!(p5,lat,lon;path_color=:black)
    p5
end

Plot combined Eastern Ontario and NAMAD maps together

2 maps can be combined at the same altitude. Here, the Eastern Ontario map is better (higher resolution), but the NAMAD map can be used in the case of navigation near the Eastern Ontario map border. Note that only the southeast part of the NAMAD map is accurate here, as the map coverage ends & the remainder is filled-in (e.g., the northwest corner).

begin
    mapS_combined = map_combine(e_mapS_395,n_mapS_395)
    p6 = plot_map(mapS_combined;dpi=dpi,use_mask=false)
    plot_path!(p6,lat,lon;path_color=:black)
    p6
end

Compare map values on border of Eastern Ontario and NAMAD maps

Map values on 2 maps will be compared for the red path along the border.

begin
    ind_trim   = (rad2deg.(lon) .> -76) .& (rad2deg.(lat) .< 45)
    lon_trim   = lon[ind_trim]
    lat_trim   = lat[ind_trim]
    p7 = plot_map(mapS_combined;dpi=dpi,use_mask=false,map_color=:gray)
    plot_path!(p7,lat,lon;path_color=:black)
    plot_path!(p7,lat_trim,lon_trim;path_color=:red)
    p7
end

First the map_interpolate function is used to create map interpolation objects, which are then evaluated at the given latitudes & longitudes (red path above). The (standard deviation) error of the map values is approximately 100 nT in this case.

begin
    e_itp_mapS = map_interpolate(e_mapS_395)
    n_itp_mapS = map_interpolate(n_mapS_395)
    e_mapS_val = e_itp_mapS.(lat_trim,lon_trim)
    n_mapS_val = n_itp_mapS.(lat_trim,lon_trim)
    mapS_err   = round(Int,std(n_mapS_val - e_mapS_val)) # error
end

The trend of the map values agree, as expected.

begin
    p8 = plot(ylab="map value [nT]",dpi=200)
    plot!(p8,e_mapS_val,lab="Eastern")
    plot!(p8,n_mapS_val,lab="NAMAD")
end

Create a 3D map

All of the functionality in MagNav.jl can use 3D maps, which contain multiple stacked 2D maps with constant altitude spacing between each map level. There are 2 ways to create a 3D map, the first of which is to upward/downward continue to multiple altitudes using the upward_fft function. In this case, a filled Perth map at 800 m is stacked with the same map upward continued to 810 m.

begin
    p_mapS3D_1 = upward_fft(p_mapS_800, p_mapS_800.alt .+ [0,10]) # 3D map
    println(p_mapS3D_1.alt) # map altitude levels
    p9 = plot_map(p_mapS3D_1) # plot single map level (default is lowest)
end

The second way to create a 3D map is to combine multiple 2D maps using the map_combine function. In this case, the Eastern Ontario map (395 m) & Perth map (800 m) are stacked into a 3D map.

begin
    p_mapS3D_2 = map_combine([e_mapS_395,p_mapS_800]; # 3D map
                             n_levels = 2,
                             dx       = get_step(p_mapS_800.xx),
                             dy       = get_step(p_mapS_800.yy),
                             xx_lim   = extrema(p_mapS_800.xx),
                             yy_lim   = extrema(p_mapS_800.yy));
    println(p_mapS3D_2.alt) # map altitude levels
    p10 = plot_map(p_mapS3D_2) # plot single map level (default is lowest)
end