World Magnetic Model#

The main utility of the World Magnetic Model (WMM) [WMM] is to provide magnetic declination for any desired location on the globe.

In addition to the magnetic declination, the WMM also provides the complete geometry of the field from 1 km below the World Geodetic System (WGS 84) [WGS84] ellipsoid surface to approximately 850 km above it. The magnetic field extends deep into the Earth and far out into space, but the WMM is not valid at these extremes.

Earth’s magnetic field is viewed as a magnetic dipole produced by a sphere of uniform magnetization. The “south” of the dipole lies on the northern Hemipshere and drifts westward at about 55 to 60 km per year, whereas the other pole, just outside the Antarctic Circle, drifts by about 10 to 15 km per year.

The strongest contribution to Earth’s magnetism is the magnetic field produced by the Earth’s liquid-iron outer core, called the “core field”. Magnetic minerals in the crust and upper mantle make a further contribution that can be locally significant. All these fields of “internal” origin and their large scale components are included in the WMM.

“External” magnetic fields, arising from electric currents in the upper atmosphere and near-Earth space, predominantly generated by solar wind, are time-varying and produce secondary magnetic fields, which are not represented in the WMM.

Current estimations include secondary thermal convection currents near Earth’s core-mantle producing local magnetic dipoles. These dipoles are superimposed to approximate the observed multipole nature of the magnetic field. The creation and decay of these inner convections are suspected to yield the magnetic drift.

The magnetic flux density (magnetic field) \(B\) can be expressed as the gradient of a potential, \(V\):

\[B = -\nabla V\]

This has a convenient solution expressed in spherical harmonics of \(n\) degrees:

\[V(r, \theta, \phi) = a \sum_{n=1}^k\Big(\frac{a}{r}\Big)^{n+1}\sum_{m=0}^{n}\big(g_n^m \cos(m\phi) + h_n^m \sin(m\phi)\big) P_n^m(\theta)\]

where \(a\) is Earth’s mean radius; \(g_n^m\) and \(h_n^m\) are Gaussian coefficients of degree \(n\) and order \(m\); \(r\), \(\theta\), and \(\phi\) are the geocentric radius, coelevation and longitude in spherical polar coordinates; and \(P_n^m(\theta)\) are the associated Legendre functions [Heiskanen].

The time-dependent Gauss coefficients are estimated empirically from a least-squares fit using satellite magnetic measurements [Langel]. These coefficients are provided by the NCEI Geomagnetic Modeling Team and British Geological Survey in a file with extension COF [WMM2020].

With degree \(n=1\) only dipoles are considered. For \(n=2\) the quadrupoles, \(n=3\) the octuploles, and so on. The method of this WMM expands to degree and order 12.

The secular variation (SV) is the yearly change of the core field, which is also accounted in the WMM by a linear model. Due to unpredictable changes in the core field, the values of the WMM coefficients are updated every five years (a lustrum). The most recent version is valid from 2020 to 2024.

The geomagnetic field vector B is described by 7 elements:

Element	Definition	Units	Range
Element	Definition	Units	Min	Max
X	Northerly intensity	nT	-17000	43000
Y	Easterly intensity	nT	-18000	17000
Z	Vertical intensity (Positive downwards)	nT	-67000	62000
H	Horizontal intensity	nT	0	43000
F	Total intensity	nT	23000	67000
I	Inclination angle (a.k.a. dip angle)	degree	-90	90
D	Declination angle (a.k.a. magnetic variation)	degree	-180	180

The quantities X, Y and Z are perpendicular vectors and can be used to determine the quantities F, I and D, and viceversa.

The vertical direction is perpendicular to the horizontal plane of the WGS 84 ellipsoid model of the Earth.

At a location on the plane of a chosen horizontal coordinate system, grivation is the angle between grid north and magnetic north, i.e., the angle measured clockwise from the direction parallel to the grid’s Northing axis to the horizontal component of the magnetic field at the observer’s location.

Grivation is useful for local surveys, where location is given by grid coordinates rather than by longitude and latitude. It is dependent on the map projection used to define the grid coordinates. In general, it is estimated as:

\[GV = D - C\]

where \(C\) is the “convergence-of-meridians” defined as the clockwise angle from the northward meridional arc to the grid Northing direction.

The class WMM contains a couple of methods to load and create a geomagnetic model from a set of given coefficients. To obtain the magnetic elements at a certain spot on the Earth, we call the method magnetic_field

>>> wmm = ahrs.utils.WMM()              # Create today's magnetic model
>>> wmm.magnetic_field(10.0, -20.0)     # Magnetic field at latitude = 10°, longitude = -20°
>>> wmm.D                               # Magnetic declination [degrees]
-9.122361367239034
>>> wmm.magnetic_field(10.0, -20.0, height=10.5)     # 10.5 km above sea level
>>> wmm.D
-9.128404039098971

By default, the class WMM will create a model for the day, when the object is being created. To ask for the values at a different date, we simply pass it as a parameter.

>>> wmm.magnetic_field(10.0, -20.0, height=10.5, date=datetime.date(2017, 5, 12))    # on 12th May, 2017
>>> wmm.D
-9.73078560629778

All main elements are computed at the same time and accessed independently

>>> wmm.X
30499.640469609083
>>> wmm.Y
-5230.267158472566
>>> wmm.Z
-1716.633311360368
>>> wmm.H
30944.850352270452
>>> wmm.F
30992.427998627096
>>> wmm.I
-3.1751692563622993
>>> wmm.GV
-9.73078560629778

or in a dictionary

>>> wmm.magnetic_elements
{'X': 30499.640469609083, 'Y': -5230.267158472566, 'Z': -1716.633311360368,
'H': 30944.850352270452, 'F': 30992.427998627096, 'I': -3.1751692563622993,
'D': -9.73078560629778, 'GV': -9.73078560629778}

Note

The model in this package includes coefficients for dates between 2015 and 2024 only. Models out of this timespan cannot be built.

The WMM was developed jointly by the National Centers for Environmental Information (NCEI, Boulder CO, USA) (formerly National Geophysical Data Center (NGDC)) and the British Geological Survey (BGS, Edinburgh, Scotland). As part of the regular update cycle of the World Magnetic Model both institutions have released the latest model on December 10th, 2019.

This script is based on the originally conceived one by Christopher Weiss (cmweiss@gmail.com) [Weiss], who adapted it from the geomagc software and World Magnetic Model of the NOAA Satellite and Information Service, National Geophysical Data Center [Chulliat].

License#

The WMM source code and binaries are in the public domain and not licensed or under copyright. The information and software may be used freely by the public. As required by 17 U.S.C. 403, third parties producing copyrighted works consisting predominantly of the material produced by U.S. government agencies must provide notice with such work(s) identifying the U.S. Government material incorporated and stating that such material is not subject to copyright protection.

References

[Chulliat]

Chulliat, A., W. Brown, P. Alken, C. Beggan, M. Nair, G. Cox, A. Woods, S. Macmillan, B. Meyer and M. Paniccia, The US/UK World Magnetic Model for 2020-2025: Technical Report, National Centers for Environmental Information, NOAA, doi:10.25923/ytk1-yx35, 2020. (https://www.ngdc.noaa.gov/geomag/WMM/data/WMM2020/WMM2020_Report.pdf)

[Heiskanen]

1. Heiskanen and H. Moritz. Physical Geodesy. TU Graz. 1993.

[Langel]

R. A. Langel and W. J. Hinze. The Magnetic Field of Earth’s Lithosphere: The Satellite Perspective. Cambridge University Press. 1998.

[Weiss]

Christopher Weiss’ GeoMag repository (cmweiss/geomag)

[Wertz]

James R. Wertz. Spacecraft Attitude Determination and Control. Kluwer Academics. 1978.

[WGS84]

World Geodetic System 1984. Its Definition and Relationships with Local Geodetic Systems. National Geospatial-Intelligence Agency (NGA) Standarization Document. 2014. ftp://ftp.nga.mil/pub2/gandg/website/wgs84/NGA.STND.0036_1.0.0_WGS84.pdf

[WMM]

The World Magnetic Model (https://www.ncei.noaa.gov/products/world-magnetic-model)

[WMM-TR]

Chulliat, A. et al. (2020). The US/UK World Magnetic Model for 2020-2025 : Technical Report. https://doi.org/10.25923/ytk1-yx35

[WMM2020]

WMM2020 Model values: NCEI Geomagnetic Modeling Team and British Geological Survey. 2019. World Magnetic Model 2020. NOAA National Centers for Environmental Information. doi: 10.25921/11v3-da71, 2020.

class ahrs.utils.wmm.WMM(date: date | int | float | None = None, latitude: float = 48.137154, longitude: float = 11.576124, height: float = 0.519, frame: str = 'NED')#

This class computes all elements of the World Magnetic Model (WMM) at any given point on Earth.

The main magnetic field \(B\) is a potential field defined, in geocentric spherical coordinates (longitude \(\lambda\), latitude \(\phi '\) and radius \(r\)), as the negative spatial gradient of a scalar potential at a time \(t\). This potential can be expanded in terms of spherical harmonics:

\[V(\lambda, \phi', r, t) = a\sum_{n=1}^{N}\Big(\frac{a}{r}\Big)^{n+1}\sum_{m=0}^{n}f(n, m, \lambda, t)P_n^m(\phi')\]

where

\[f(n, m, \lambda, t) = g_n^m(t) \cos(m\lambda) + h_n^m(t) \sin(m\lambda)\]

and the Schmidt semi-normalized associated Legendre functions \(P_n^m(\phi')\) are defined as:

\[\begin{split}P_n^m(\mu) = \left\{ \begin{array}{ll} \sqrt{2\frac{(n-m)!}{(n+m)!}}P_{n, m}(\mu) & \mathrm{if} \; m > 0 \\ P_{n, m}(\mu) & \mathrm{if} \; m = 0 \end{array} \right.\end{split}\]

Any object of this class is initialized with the corresponding epoch, determined by the given date. If no date is given, it is assumed for the day of the object’s creation.

Once the WMM object is created, the estimation of the geomagnetic elements is carried out with a call to the method magnetic_field giving the location on Earth at which the magnetic elements will be calculated. This location is given in decimal geodetic coordinates. See examples.

Every WMM object is created with a set of coefficients read from a COF file, defined by the desired working date of the model. The latest model available is WMM2020 corresponding to the lustrum 2020-2024.

This class can create models with dates between 2015 and 2024.

Parameters:

date (datetime.date, int or float, default: current day) – Date of desired magnetic field estimation.
latitude (float, default: None) – Latitude, in decimal degrees, in geodetic coordinates.
longitude (float, default: None) – Longitude, in decimal degrees, in geodetic coordinates.
height (float, default: 0.0) – Mean Sea Level Height, in kilometers.
frame (str, default: 'NED') – Local tangent plane coordinate frame. Valid options are right-handed 'NED' for North-East-Down and 'ENU' for East-North-Up.

Variables:

date (datetime.date, default: datetime.date.today()) – Desired date to estimate
date_dec (float) – Desired date to estimate as decimal
epoch (float) – Initial time of model in decimal years
model (str) – WMM Model identificator
modeldate (str) – Release date of WMM Model
wmm_filename (str) – COF File used to build Model
degree (int) – Degree of model
latitude (float) – Latitude, in decimal degrees, in geodetic coordinates
longitude (float) – Longitude in decimal degrees, in geodetic coordinates
height (float, default: 0.0) – Mean Sea Level Height in kilometers
X (float, default: None) – Northerly intensity, in nT
Y (float, default: None) – Easterly intensity, in nT
Z (float, default: None) – Vertical intensity, in nT
H (float, default: None) – Horizontal intensity, in nT
F (float, default: None) – Total intensity, in nT
I (float, default: None) – Inclination angle (dip), in degrees
D (float, default: None) – Declination angle, in degrees
GV (float, default: None) – Grivation, in degrees

Examples

The magnetic field can be computed at the creation of the WMM object by passing the main parameters to its constructor:

>>> wmm = ahrs.utils.WMM(datetime.date(2017, 5, 12), latitude=10.0, longitude=-20.0, height=10.5)
>>> wmm.magnetic_elements
{'X': 30499.640469609083, 'Y': -5230.267158472566, 'Z': -1716.633311360368,
'H': 30944.850352270452, 'F': 30992.427998627096, 'I': -3.1751692563622993,
'D': -9.73078560629778, 'GV': -9.73078560629778}

denormalize_coefficients(latitude: float) → None#

Recursively estimate associated Legendre polynomials and derivatives done in a recursive way as described by Michael Plett in [Wertz] for an efficient computation.

Given the Gaussian coefficients, it is possible to estimate the magnetic field at any latitude on Earth for a certain date.

First, it is assumed that \(P_n^m(x)\) are the Schmidt semi-normalized functions. A normalization is made so that the relative strength of terms of same degree \(n\) but order \(m\) are used by comparing their respective Gauss coefficients.

For \(m=0\) they are called Legendre Polynomials and can be computed recursively with:

\[P_n(x) = \frac{2n-1}{n} x P_{n-1}(x) - \frac{n-1}{n}P_{n-2}(x)\]

For \(m>0\) they are known as associated Legendre functions of degree \(n\) and order \(m\) and reduced to:

\[P_{nm}(x) = (1-t^2)^{m/2} \frac{d^m P_n(x)}{dt^m}\]

expressing the associated Legendre functions in terms of the Legendre polynomials of same degree \(n\).

A more general formula to estimate both polynomial and associated functions is given by:

\[P_{nm}(x) = 2^{-n}(1-x^2)^{m/2} \sum_{k=0}^{K}(-1)^k\frac{(2n-2k)!}{k!(n-k)!(n-m-2k)!}x^{n-m-2k}\]

where \(K\) is either \((n-m)/2\) or \((n-m-1)/2\), whichever is an integer. For a computational improvement, the terms are calculated recursively.

We have to denormalize the coefficients from Schmidt to Gauss. The Gauss functions \(P^{n, m}\) are related to Schmidt functions \(P_n^m\) as:

\[P_n^m = S_{n, m}P^{n, m}\]

where the factors \(S_{n, m}\) are combined with Gaussian coefficients to accelerate the computation, because they are independent of the geographic location. Thus, we denormalize the coefficients with:

\[\begin{split}\begin{array}{ll} g^{n,m} & = S_{n,m} g_n^m \\ h^{n,m} & = S_{n,m} h_n^m \end{array}\end{split}\]

The recursion for \(S_{n, m}\) is:

\[\begin{split}\begin{array}{rlr} S_{0,0} & = 1 & \\ S_{n,0} & = S_{n-1, 0} \frac{2n-1}{n} & n\geq 1 \\ S_{n,m} & = S_{n-1, m}\sqrt{\frac{(n-m+1)(\delta _m^1+1)}{n+m}} & m\geq 1 \end{array}\end{split}\]

where the Kronecker delta \(\delta_j^i\) is:

\[\begin{split}\delta_j^i = \left\{ \begin{array}{ll} 1 & \: i = j \\ 0 & \: \mathrm{otherwise} \end{array} \right.\end{split}\]

Similarly, \(P^{n, m}(x)\) can be recursively obtained:

\[\begin{split}\begin{array}{ll} P^{0,0} & = 1 \\ P^{n,n} & = \sin (x) P^{n-1, n-1} \\ P^{n,m} & = \cos (x) P^{n-1, m} - K^{n, m} P^{n-2, m} \end{array}\end{split}\]

where:

\[\begin{split}K^{n, m} = \left\{ \begin{array}{ll} \frac{(n-1)^2-m^2}{(2n-1)(2n-3)} & \: n>1 \\ 0 & \: n=1 \end{array} \right.\end{split}\]

Likewise, the gradient \(\frac{dP^{n, m}}{dx}\) is estimated as:

\[\begin{split}\begin{array}{llr} \frac{dP^{0, 0}}{dx} & = 1 & \\ \frac{dP^{n, n}}{dx} & = \sin (x) \frac{dP^{n-1, n-1}}{dx} + \cos (x) P^{n-1, n-1} & n\geq 1 \\ \frac{dP^{n, m}}{dx} & = \cos (x) \frac{dP^{n-1, m}}{dx} - \sin (x) P^{n-1, m} - K^{n, m} \frac{dP^{n-2, m}}{dx} & \end{array}\end{split}\]

Parameters:: latitude (float) – Latitude in spherical geocentric coordinates

property geodetic_vector: ndarray#

Magnetic field vector in geodetic coordinates

The magnetic field vector is returned as a 3D vector in geodetic coordinates.

Example

>>> wmm = WMM(datetime.date(2017, 5, 12), latitude=10.0, longitude=-20.0, height=10.5)
>>> wmm.geodetic_vector
array([ 30499.64046961,  -5230.26715847,  -1716.63331136])

get_properties(cof_file: str) → Dict[str, str | float]#

Return dictionary of WMM properties from COF file.

Three properties are read and returned in a dictionary:

epoch is the initial time \(t_0\) as a float.
model is a string of model used for the required lustrum.
modeldate is the release date of used magnetic model.

Parameters:: cof_file (str) – Path to COF file with the coefficients of the WMM
Returns:: properties – Dictionary with the three WMM properties.
Return type:: dictionary

Examples

>>> wmm = ahrs.WMM()
>>> wmm.get_properties('my_coefficients.COF')
{'model': 'WMM-2020', 'modeldate': '12/10/2019', 'epoch': 2020.0}

load_coefficients(cof_file: str) → None#

Load model coefficients from COF file.

The model coefficients, also referred to as Gauss coefficients, are listed in a COF file. These coefficients can be used to compute values of the fields elements and their annual rates of change at any location near the surface of the Earth.

The COF file has 6 columns:

n is the degree.
m is the order.
g are time-dependent Gauss coefficients of degree n and order m.
h are time-dependent Gauss coefficients of degree n and order m.
gd are secular variations of coefficient g.
hd are secular variations of coefficient h.

which constitute the model of the field. The first-order time derivatives are called secular terms. The units are nT for the main field, and nT/year for the secular variation.

The Gauss coefficients are defined for a time \(t\) as:

\[\begin{split}\begin{eqnarray} g_n^m(t) & = & g_n^m(t_0) + (t-t_0) \dot{g}_n^m(t_0) \\ h_n^m(t) & = & h_n^m(t_0) + (t-t_0) \dot{h}_n^m(t_0) \end{eqnarray}\end{split}\]

where time is given in decimal years and \(t_0\) corresponds to the epoch read from the corresponding COF file.

Parameters:: cof_file (str) – Path to COF file with the coefficients of the WMM

property magnetic_elements: Dict[str, float]#

Main geomagnetic elements in a dictionary

Element	Definition
X	Northerly intensity
Y	Easterly intensity
Z	Vertical intensity (Positive downwards)
H	Horizontal intensity
F	Total intensity
I	Inclination angle (a.k.a. dip angle)
D	Declination angle (a.k.a. magnetic variation)
GV	Grivation

Example

>>> wmm = WMM(datetime.date(2017, 5, 12), latitude=10.0, longitude=-20.0, height=10.5)
>>> wmm.magnetic_elements
{'X': 30499.640469609083, 'Y': -5230.267158472566, 'Z': -1716.633311360368,
'H': 30944.850352270452, 'F': 30992.427998627096, 'I': -3.1751692563622993,
'D': -9.73078560629778, 'GV': -9.73078560629778}

magnetic_field(latitude: float, longitude: float, height: float = 0.0, date: date | int | float = datetime.date(2024, 4, 7)) → None#

Calculate the geomagnetic field elements for a location on Earth.

The code includes comments with references to equation numbers corresponding to the ones in the official report.

Having the coefficients \(g^{n, m}\) and \(h^{n, m}\), we extrapolate them for the desired time \(t\) as:

\[\begin{split}\begin{array}{ll} g_n^m(t) & = g_n^m(t_0) + \Delta_t \dot{g}_n^m (t_0) \\ h_n^m(t) & = h_n^m(t_0) + \Delta_t \dot{h}_n^m (t_0) \end{array}\end{split}\]

where \(\Delta_t = t-t_0\) is the difference between the time \(t\) and the reference epoch model \(t_0\) (2020.0 for the newest version.)

The vector components of the main magnetic field \(B\) are then calculated with:

\[\begin{split}\begin{array}{ll} X' & = -\sum_{n=1}^N\Big(\frac{a}{r}\Big)^{n+2} \sum_{m=0}^n\big(g_n^m(t) \cos(m\lambda)+h_n^m(t)\sin(m\lambda)\big) \frac{dP_n^m(\sin \phi ')}{d\phi '} \\ Y' & = \frac{1}{\cos\phi '}\sum_{n=1}^N\Big(\frac{a}{r}\Big)^{n+2} \sum_{m=0}^n m\big(g_n^m(t) \sin(m\lambda)-h_n^m(t)\cos(m\lambda)\big)P_n^m(\sin \phi ') \\ Z' & = -\sum_{n=1}^N(n+1)\Big(\frac{a}{r}\Big)^{n+2} \sum_{m=0}^n\big(g_n^m(t) \cos(m\lambda)+h_n^m(t)\sin(m\lambda)\big)P_n^m(\sin \phi ') \end{array}\end{split}\]

Finally, the geomagnetic vector components are rotated into ellipsoidal reference frame.

\[\begin{split}\begin{array}{ll} X & = X'\cos(\phi ' - \phi) - Z' \sin(\phi ' - \phi) \\ Y & = Y' \\ Z & = X'\sin(\phi ' - \phi) + Z' \cos(\phi ' - \phi) \end{array}\end{split}\]

These components are used to compute the rest of the magnetic elements:

\[\begin{split}\begin{array}{ll} H & = \sqrt{X^2 + Y^2} \\ F & = \sqrt{H^2 + Z^2} \\ I & = \arctan(\frac{Z}{H}) \\ D & = \arctan(\frac{Y}{X}) \end{array}\end{split}\]

Note

The use of arctan2 yields a more precise result than arctan, because it estimates the angle exploring all quadrants.

For polar regions, where the declination changes drastically, the WMM defines two different grivations (one for each pole) defined as:

\[\begin{split}GV = \left\{ \begin{array}{ll} D-\lambda & \: \phi > 55 ° \\ D+\lambda & \: \phi < -55 ° \end{array} \right.\end{split}\]

Parameters:

latitude (float) – Latitude, in decimal degrees, in geodetic coordinates
longitude (float) – Longitude in decimal degrees, in geodetic coordinates
height (float, default: 0.0) – Mean Sea Level Height in kilometers
date (datetime.date, int or float, default: datetime.date.today()) – Desired date to estimate

reset_coefficients(date: date | int | float | None = None) → None#

Reset Gauss coefficients to given date.

Given the date, the corresponding coefficients are updated. Basic properties (epoch, release date, and model id) are read and updated in the current instance.

The two coefficient tables (arrays) are also updated, where the attribute c contains the Gaussian coefficients, while the attribute cd contains the secular Gaussian coefficients.

The lenght of the Gaussian coefficient array determines the degree \(n\) of the model. This property updates the value of attribute degree.

Parameters:: date (datetime.date, int or float, default: current day) – Date of desired magnetic field estimation.

reset_date(date: date | int | float) → None#

Set date to use with the model.

The WMM requires a date. This date can be given as an instance of datetime.date or as a decimalized date of the format YYYY.d.

If None is given it sets the date to the present day. In addition, the corresponding COF file is also set.

Please note that only coefficents from year 2015 and later are provided with this module.

Parameters:: date (datetime.date, int or float, default: current day) – Date of desired magnetic field estimation.

ahrs.utils.wmm.geodetic2spherical(lat: float, lon: float, h: float, a: float = 6378.137, b: float = 6356.7523142) → Tuple[float, float, float]#

Transform geodetic coordinates into spherical geocentric coordinates

The transformation cannot be a simple cylindric to spherical conversion, as we must also consider a planet’s ellipsoid form. With the aid of its pre-defined flatness and eccentricity, we can better approximate the values of the conversion.

In this function the Earth’s major and minor semi-axis are considered. However, we can convert the coordinates of different ellipsoidal bodies, by giving the dimensions of its semi-axes.

Notice that the longitude between both systems remains the same.

Parameters:

lat (float) – Latitude, in radians, of point in geodetic coordinates
lon (float) – Longitude, in radians, of point in geodetic coordinates
h (float) – Height, in kilometers, of point in geodetic coordinates
a (float, default: 6378.137) – Major semi-axis, in kilometers. Defaults to Earth’s equatorial radius
b (float, default: 6356.752314245) – Minor semi-axis, in kilometers. Defaults to Earth’s polar radius

Returns:

lat_spheric (float) – Latitude of point in spherical coordinates.
lon (float) – Longitue of point in spherical coordinates. Same as geodetic.
r (float) – Radial distance of point in spherical coordinates.