Skip to contents

Compute Magnetic Likelihood Maps

Source: R/geomag_map.R
geomag_map.Rd

Estimates a likelihood map for each stationary period based on observed magnetic data (intensity and inclination) compared to the expected values from the World Magnetic Model (WMM).

Usage

geomag_map(
  tag,
  compute_known = FALSE,
  sd_e_f = 0.009,
  sd_e_i = 2.6,
  sd_m_f = 0.014,
  sd_m_i = 3.5,
  ref_map = geomag_map_ref(tag),
  quiet = FALSE
)

Arguments

tag

A calibrated GeoPressureR tag object with map metadata set via tag_set_map().

compute_known

Logical. If TRUE, computes likelihood maps for known stationary periods; if FALSE, fixes the likelihood at the known location.

sd_e_f

Numeric. Standard deviation of observation noise of intensity error (single value or one per stap).

sd_e_i

Numeric. Standard deviation of observation noise of inclination error (single value or one per stap).

sd_m_f

Numeric. Standard deviation of stationary-period-specific noise of intensity error (single value or per stap).

sd_m_i

Numeric. Standard deviation of stationary-period-specific noise of inclination error (single value or per stap).

ref_map

A SpatRaster with layers intensity and inclination. By default this is computed by geomag_map_ref(). Provide it if you've already computed it to save computational time.

quiet

Logical. If TRUE, suppresses progress messages.

Value

A GeoPressureR tag object with likelihood maps added as:

  • tag$map_magnetic_intensity

  • tag$map_magnetic_inclination

  • tag$map_magnetic and updated parameters in tag$param.

Details

The tag must already contain calibrated magnetic data and map metadata created with tag_set_map(). geomag_map() uses the tag's pressure record to estimate altitude for the WMM reference map. Stationary-period identifiers must be consecutive integers from 1 to the number of stationary periods.

For an observed magnetic quantity \(y_j\) and WMM value \(x_s\) at location \(s\), the observation model is \(y_j = x_s + b + e_j\), with observation error \(e_j \sim N(0, \sigma_e^2)\) and stationary-period error \(b \sim N(0, \sigma_m^2)\). The likelihood integrates over \(b\) and is returned relative to its maximum value. When no valid observations are available, an uninformative all-one map is returned.

Intensity values and their standard deviations (sd_e_f, sd_m_f) are in Gauss. Inclination values and their standard deviations (sd_e_i, sd_m_i) are in degrees. Sensitivity analyses should vary both error components because larger values broaden the likelihood surface and reduce the spatial information contributed by magnetic data.

Examples

library(GeoPressureR)
withr::with_dir(system.file("extdata", package = "GeoMagR"), {
  tag <- tag_create("14DM", quiet = TRUE)
  tag <- tag_label(tag, quiet = TRUE)
  tag <- tag_set_map(tag,
    extent = c(-18, 23, 0, 50),
    scale = 2,
    known = data.frame(
      stap_id = c(1, -1),
      known_lon = 7.27,
      known_lat = 46.19
    )
  )
  tag <- geomag_calib(tag, quiet = TRUE)
  tag <- geomag_map(tag, quiet = TRUE)
  plot(tag, type = "map_magnetic_intensity")
  plot(tag, type = "map_magnetic_inclination")
})
#> ℹ Using raw magnetic data for calibration data
#> Calculating Magnetic Field map ■■■■■■■■■■■■                      36% | ETA:  2s
#> Calculating Magnetic Field map ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s