Add IGRF forward calculation

[leouieda]

Add a class to calculate the IGRF14 field. Uses Pooch to download IGRF14 from Zenodo. The class takes a date to instantiate and then calculates the field for that data through predict and grid methods. Inputs and outputs are in geocentric spherical coordinates by default but can be changed to geodetic if a Boule ellipsoid is given to the constructor.

Relevant issues/PRs: Fixes #504

[leouieda]

Running a benchmark against ppigrf:

Grid size: (594, 1201) = 713394
Harmonica: 0.156 +- 0.447 s 
ppigrf: 9.998 +- 0.163 s 
Harmonica is 63.89 times faster.

The Harmonica results include the numba compilation, hence the large std. Results are median times. This is the code for the benchmark:

import numpy as np
import time
import datetime
import ppigrf
import harmonica as hm
import verde as vd


date = datetime.datetime.fromisoformat("2024-01-20")
coords = vd.grid_coordinates((0, 360, -89, 89), spacing=0.3, extra_coords=0)
print(f"Grid size: {coords[0].shape} = {coords[0].size}")

times_harmonica = []
for i in range(10):
    start = time.time()
    igrf = hm.IGRF14(date)
    b_harmonica = igrf.predict(coords)
    times_harmonica.append(time.time() - start)
median_harmonica = np.median(times_harmonica)
print(f"Harmonica: {median_harmonica:.3f} +- {np.std(times_harmonica):.3f} s ")

times_ppigrf = []
for i in range(5):
    start = time.time()
    b_ppigrf = ppigrf.igrf(*coords, date)
    times_ppigrf.append(time.time() - start)
median_ppigrf = np.median(times_ppigrf)
print(f"ppigrf: {median_ppigrf:.3f} +- {np.std(times_ppigrf):.3f} s ")

print(f"Harmonica is {median_ppigrf / median_harmonica:.2f} times faster.")

# Check that all components are close for both software
np.testing.assert_allclose(b_harmonica, [c[0] for c in b_ppigrf], rtol=0.001)

[leouieda]

Would be good to benchmark against SHTools as well.

[leouieda]

Here's an updated benchmark against ppigrf and SHTools:

Grid size: (601, 1201) = 721801
Harmonica       : 0.158 +- 0.451 s 
Harmonica (grid): 0.075 +- 0.315 s 
SHTools         : 0.015 +- 0.125 s 
ppigrf          : 10.253 +- 0.130 s 

Again, what's reported is the median time. So we're way faster than ppigrf but SHTools is still king (as expected). I'm actually surprised that we're not that much slower than SHTools but I think that'll change drastically if we do higher degree models. Still, the grid implementation is able to but down the time in half by avoid repeat computations.

This is the benchmark script:

import numpy as np
import warnings
import time
import datetime
import ppigrf
import harmonica as hm
import verde as vd
import pyshtools as psh
import boule as bl


# ppigrf doesn't handle the poles very well
warnings.filterwarnings("ignore")

date = datetime.datetime.fromisoformat("1900-01-01")
n = 600
region = (0, 360, -90, 90)
shape = (n + 1, 2 * n + 1)
height = 0
coords = vd.grid_coordinates(region, shape=shape, extra_coords=height)
print(f"Grid size: {coords[0].shape} = {coords[0].size}")

times_harmonica = []
for _ in range(10):
    start = time.perf_counter()
    igrf = hm.IGRF14(date)
    b_harmonica = igrf.predict(coords)
    times_harmonica.append(time.perf_counter() - start)
median_harmonica = np.median(times_harmonica)
print(f"Harmonica       : {median_harmonica:.3f} +- {np.std(times_harmonica):.3f} s ")

times_harmonica_g = []
for _ in range(10):
    start = time.perf_counter()
    igrf = hm.IGRF14(date)
    b_harmonica_g = igrf.grid(region, height, shape=shape)
    times_harmonica_g.append(time.perf_counter() - start)
median_harmonica_g = np.median(times_harmonica_g)
print(
    f"Harmonica (grid): {median_harmonica_g:.3f} +- {np.std(times_harmonica_g):.3f} s "
)

times_shtools = []
for _ in range(10):
    start = time.perf_counter()
    igrf = psh.datasets.Earth.IGRF_13(year=date.year)
    b_shtools = igrf.expand(
        a=bl.WGS84.semimajor_axis,
        f=bl.WGS84.flattening,
        lmax=(n - 2) / 2,
        lmax_calc=13,
    )
    times_shtools.append(time.perf_counter() - start)
median_shtools = np.median(times_shtools)
print(f"SHTools         : {median_shtools:.3f} +- {np.std(times_shtools):.3f} s ")

times_ppigrf = []
for _ in range(5):
    start = time.perf_counter()
    b_ppigrf = ppigrf.igrf(*coords, date)
    times_ppigrf.append(time.perf_counter() - start)
median_ppigrf = np.median(times_ppigrf)
print(f"ppigrf          : {median_ppigrf:.3f} +- {np.std(times_ppigrf):.3f} s ")

This was as much as my computer could handle in terms of RAM from ppigrf. It used almost 8 Gb to calculate on the ~700k points because of all the numpy broadcasting and pandas operations they do. Harmonica and SHTools used virtually nothing.

[leouieda]

Need to do:

• Implement more tests, particularly the grid method.
• Add metadata to the generated grid.
• Automatically calculate the ideal spacing based on degree.
• Finish docstrings.
• Add the class to the API docs.
• Create an example.

[leouieda]

@santisoler @mdtanker this is ready for a review. Since you both commented on #504, it would be great to have your thoughts on the API in particular.

[mdtanker]

This looks great! I added a few comments to the User Guide for some small and totally optional changes, which are included in PR #582. I also noticed a small mistake in the docs (or confusion on my part) about the shape of the coefficients attribute.

I tested this with yet another Python implementation (https://github.com/zzyztyy/pyIGRF), which had good agreement.

Thanks for implementing this, I will definitely be using it!

[leouieda]

Thanks @mdtanker! I added all your suggestions.