‹‹‹ Blog Overview

Where to start?

This discussion on Github beautifully summarizes the problem: "Too many IGRFs"

In all seriousness, it is a lot of implementations, at least 10 as of writing this text. This does not need to be a bad thing - competition is usually good and I like when science gets reproduced. However, what I found when I entered this arena was that few of the many implementations had been properly tested or verified. Besides, almost everyone of them had its smaller or bigger issues, anything from structure, design, speed, packaging, documentation etc. pyIGRF by zzyztyy was by far not the best, but it was (a) the one closest to IAGA's original Fortran implementation and, given the previous constraint, (b) the most "pure Python" one. So I based my work on it. I needed something I could trust, fast and verified.

What are the options?

There are other, independent implementations named pyIGRF, not to be confused with the previously mentioned one:

• another recent independent pyIGRF, Python
• another older independent pyIGRF, Python & Fortran
• another small independent pyigrf, Python

The official IGRF implementations by IAGA are available via the NOAA website:

• IAGA IGRF, Fortran
• IAGA IGRF, Python

A good, modern and pure Python implementation:

• ppigrf, Python

Another Python implementation is part of the navtools package:

• navtools IGRF, Python

Another Python wrapper around a cleaned-up version of IAGA's Fortran code is maintained as part of the space physics project:

• space-physics IGRF, Python & Fortran

A more or less undocumented, pure Python implementation can be found here:

• ADCS_IGRF, Python

All of the above models rely on the IGRF-13 coefficients¹ as provided by NOAA.

The IGRF is not to be confused with the World Magnetic Model (WMM). Recent implementations of the WMM accessible via Python are maintained as part of the space physics project:

• space-physics WMM2015, Python & C
• space-physics WMM2020, Python & C

A new implementation

Just to make a small difference, I named my implementation pyCRGI, using the French translation of IGRF's name: Champ de Référence Géomagnétique International.

pyCRGI is a massively cleaned-up and modernized fork of pyIGRF by zzyztyy. Be aware that there are a number of small function and module name differences to the original pyIGRF package. The fork's main goals are verified results, tests, speed and ease of maintainability. This is work in progress. The package offers different implementations, pure Python, Python JIT-compiled via numba and array-oriented via numpy and numba.

The current project status can be summarized as follows:

• [x] package structure cleanup, allowing proper testing and packaging
• [x] type annotations
• [x] doc strings completed and prepared for Sphinx autodoc
• [x] debug mode via environment variable PYCGIR_DEBUG=1
• [x] pure Python 3 implementation without dependency to numpy
• [x] JIT-compiled implementation depending on numba and numpy, installation target jited
• [x] array implementation depending on numba and numpy, installation target array
• [x] unit tests via test makefile target
• [x] verification of field's values against original Fortran implementation
• [ ] verification of field's annual/seasonal variation -> huge differences to original Fortran implementation
• [x] benchmark
• [ ] clean up and add missing relevant coordinate system conversions

How does it compare?

How to Install?

Use pip to install the latest development version from Github:

pip install git+https://github.com/pleiszenburg/pyCRGI.git@master

How to Use it?

First import the package, either as pure Python 3 or JIT-compiled via numba and numpy:

from pyCRGI.pure import get_value, get_variation # pure Python 3
# or
from pyCRGI.jited import get_value, get_variation # JIT-compiled via `numba` and `numpy`
# or
from pyCRGI.array import get_value, get_variation # array implementation via `numba` and `numpy`

You can calculate the magnetic field's intensity:

get_value(lat, lon, alt, year)

You can calculate the annual variation of the magnetic field's intensity:

get_variation(lat, lon, alt, year)

The return value is a tuple (or array) of seven floating point numbers representing the local magnetic field:

• D: declination (+ve east) [degree]
• I: inclination (+ve down) [degree]
• H: horizontal intensity [nT]
• X: north component [nT]
• Y: east component [nT]
• Z: vertical component (+ve down) [nT]
• F: total intensity [nT]

If you want to use the IGRF-13 model in a more flexible manner, you can use the functions geodetic2geocentric and get_syn. They are somewhat closer to the original Fortran implementation.

Another function, get_coeffs, can be used to get g[m][n] or h[m][n] corresponding to the IGRF's formula.

‹‹‹ Blog Overview