Magnetic Anomaly Maps
This file is best viewed in a Pluto notebook. To do so, from the MagNav.jl directory, run:
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)
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")
Overlay Perth mini-survey
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