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 liquidiron 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 nearEarth space, predominantly generated by solar wind, are timevarying and produce secondary magnetic fields, which are not represented in the WMM.
Current estimations include secondary thermal convection currents near Earth’s coremantle 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\):
This has a convenient solution expressed in spherical harmonics of \(n\) degrees:
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 timedependent Gauss coefficients are estimated empirically from a leastsquares 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  

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:
where \(C\) is the “convergenceofmeridians” 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 20202025: Technical Report, National Centers for Environmental Information, NOAA, doi:10.25923/ytk1yx35, 2020. (https://www.ngdc.noaa.gov/geomag/WMM/data/WMM2020/WMM2020_Report.pdf) 
[Heiskanen] 

[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 (https://github.com/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 GeospatialIntelligence 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.ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml) 
[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/11v3da71, 2020. 

class
ahrs.utils.wmm.
WMM
(date: Union[datetime.date, int, float] = None, latitude: float = None, longitude: float = None, height: float = 0.0)¶ World Magnetic Model
It is mainly used to compute 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 seminormalized associated Legendre functions \(P_n^m(\phi')\) are defined as:
\[\begin{split}P_n^m(\mu) = \left\{ \begin{array}{ll} \sqrt{2\frac{(nm)!}{(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 20202024.
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.
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) → NoReturn¶ 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 seminormalized 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{2n1}{n} x P_{n1}(x)  \frac{n1}{n}P_{n2}(x)\]For \(m>0\) they are known as associated Legendre functions of degree \(n\) and order \(m\) and reduced to:
\[P_{nm}(x) = (1t^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}(1x^2)^{m/2} \sum_{k=0}^{K}(1)^k\frac{(2n2k)!}{k!(nk)!(nm2k)!}x^{nm2k}\]where \(K\) is either \((nm)/2\) or \((nm1)/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_{n1, 0} \frac{2n1}{n} & n\geq 1 \\ S_{n,m} & = S_{n1, m}\sqrt{\frac{(nm+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^{n1, n1} \\ P^{n,m} & = \cos (x) P^{n1, m}  K^{n, m} P^{n2, m} \end{array}\end{split}\]where:
\[\begin{split}K^{n, m} = \left\{ \begin{array}{ll} \frac{(n1)^2m^2}{(2n1)(2n3)} & \: 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^{n1, n1}}{dx} + \cos (x) P^{n1, n1} & n\geq 1 \\ \frac{dP^{n, m}}{dx} & = \cos (x) \frac{dP^{n1, m}}{dx}  \sin (x) P^{n1, m}  K^{n, m} \frac{dP^{n2, m}}{dx} & \end{array}\end{split}\]Parameters: latitude (float) – Latitude in spherical geocentric coordinates

get_properties
(cof_file: str) → Dict[str, Union[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': 'WMM2020', 'modeldate': '12/10/2019', 'epoch': 2020.0}

load_coefficients
(cof_file: str) → NoReturn¶ 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 timedependent Gauss coefficients of degreen
and orderm
.h
are timedependent Gauss coefficients of degreen
and orderm
.gd
are secular variations of coefficientg
.hd
are secular variations of coefficienth
.
which constitute the model of the field. The firstorder time derivatives are called secular terms. The units are
nT
for the main field, andnT/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) + (tt_0) \dot{g}_n^m(t_0) \\ h_n^m(t) & = & h_n^m(t_0) + (tt_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

magnetic_elements
¶ 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: Union[datetime.date, int, float] = datetime.date(2020, 11, 23)) → NoReturn¶ 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 = tt_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 thanarctan
, 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: Union[datetime.date, int, float] = None) → NoReturn¶ 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: Union[datetime.date, int, float]) → NoReturn¶ 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 predefined flatness and eccentricity, we can better approximate the values of the conversion.
In this function the Earth’s major and minor semiaxis are considered. However, we can convert the coordinates of different ellipsoidal bodies, by giving the dimensions of its semiaxes.
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 semiaxis, in kilometers. Defaults to Earth’s equatorial radius
 b (float, default: 6356.752314245) – Minor semiaxis, 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.