IGRF calculation

[leouieda]

Description of the desired feature:

I'd like to implement calculation of the IGRF field. This would require the following components:

1. A way to calculate the associated Legendre functions $P_n^m$ without the Condon-Shortly phase and with Schmidt quasi-normalization. Equations for this can be found in https://www.gnu.org/software/gsl/tr/tr001.pdf
2. Reading the model coefficients from a file.
3. Interpolating the model coefficients to the specified date. Not hard to do with Python's datetime module.

Since the model is only degree 13, we don't have to worry too much about the Legendre functions blowing up at larger orders. It would be nice to have the Holmes and Featherstone (2002) method implemented later but it's not required for this job.

Design-wise, I have a few questions:

1. Should this be a function or a class? An igrf function would be very straight forward: calculate the 3-components at the given time and location. The location can be a float or array. But this won't produce an xarray grid. A class could have IGRF.predict and IGRF.grid like our equivalent-sources. But it won't have a fit since the model is already fit. In that case, IGRF could be an instance of a SphericalHarmonics class but I'm not sure I like that.
2. The calculations are done in geocentric spherical coordinates. But usually we'll have geodetic coordinates as inputs and the vector components output should be rotated to a local geodetic frame as well. At least, this is what web calculators and the NOAA and BGS programs calculate. Doing the conversion requires knowing the ellipsoid and so Boule comes in. So: should the function/class require an bould.Ellipsoid to be passed so it can do the conversions? Should it be optional (with ellipsoid=None assuming values in geocentric coordinates)? Or do we do everything in spherical and delegate the conversions entirely to Boule outside the IGRF function (this would require adding a vector rotation function to boule.Ellipsoid)?
3. Do we package the coefficient file or download it from NOAA? I'd be hesitant of NOAA changing their download location or simply changing the file when a new IGRF comes out. So maybe packaging it would be best.

I might break this down into multiple PRs adding different part as possibly private functions (Legendre and coefficient reading ones).

Are you willing to help implement and maintain this feature?

Yes!

[santisoler]

Excited to see this.

A few shallow comments after you ideas:

1. I have no personal preference between having a igrf function or a IGRF class. I would consider how often we would compute the IGRF on a grid and if we want it to be an xarray object. In case we choose it to be a class, I would avoid creating premature abstractions and just implement a single class for it without inheritance. If in the future we find use for a SphericalHarmoincs class, we can refactor it and introduce the inheritance later on.
2. I think we should consider geodetic coordinates as the default, and to simplify the process, I would make this new class/function to handle the rotation by taking an ellipsoid as argument. By default we could use WGS84. In case we choose to implement the IGRF class, should ellipsoid be an argument of the constructor of the predict/grid method?
3. Those files aren't too large, I would package them. Also, packaging them allows us to use/test the feature even without internet access.

Bonus: until we implement the Holmes and Featherstone (2002) method, I think the functions for computing the Legendre polynomials should include a warning (maybe just in the docs) informing users to avoid using the function for higher orders.

[MarkWieczorek]

A couple comments:

1. Should this implementation be specific to IGRF, or more general for an arbitrary magnetic field model? For Earth, there are higher order models that include the lithospheric field, and we can also include planetary magnetic field models of Mercury, the Moon, and Mars. Nevertheless, I understand that these other models don't have a time dependence (Mercury might after BepiColombo). A class would probably be easier to deal with IGRF specifically, because otherwise, you would need to (a) compute the spherical harmonic coefficients at a given time, (b) evaluate the function at a position or grid, and (c) convert from spherical to geodetic coordinates.
2. For Boule, we should at some point consider implementing functions that give the normal and tangent vectors to the ellipsoid. With these, it would be simple to convert vectors from one reference frame to another.
3. I would definitely include the Holmes and Featherstone scaling, as it is really only a line or two of code. If you want inspiration, here is my version in fortran: https://github.com/SHTOOLS/SHTOOLS/blob/master/src/PlmSchmidt.f95. Once you cut out all the boiler plate, its pretty short.
4. If its any help, I wrote this routine to read IGRF 11-13 and output the coefficients at a given time: https://github.com/SHTOOLS/SHTOOLS/blob/master/pyshtools/shio/read_igrf.py
5. If you don't like relying on the NOAA url, I would just upload the file to zenodo and then use pooch to download it. Perhaps someone already did this.

[LL-Geo]

Looks great! 👍

Just noticed a Python package on the NOAA website https://www.ngdc.noaa.gov/IAGA/vmod/pyIGRF.zip

I think pack the coefficient will be very reasonable, which makes life easier without the internet. It is 40 kb, so it won't be an issue.

[abbyazari]

Commenting on Mark W.'s thoughts re: IGRF vs SH generally "Should this implementation be specific to IGRF, or more general for an arbitrary magnetic field model?"

-> It would greatly enable general use if the class was written more arbitrarily. I think as a first step though obviously IGRF functionality given the user base but the more that can be done to enable general use for arbitrary spherical harmonics in the future would be useful and let folks use different versions of SH models. It is worth noting that a fair amount of the tools in SHTools (like Mark pointed out) provide some of this base functionality already.

[leouieda]

Hi all, thank you for the great inputs! To answer some of them:

A lot of this is in SHTools already

Yes, absolutely. And we have no intention of trying to replicate more than necessary. The niche we aim to cover is for evaluating SH models at arbitrary points (which SH tools can do but not as conveniently as we'd like) and doing inversions (synthesis but with non-gridded input data).

Be just IGRF or more broadly used

We're planning on adding IGRF calculation first and then generalizing this to gravity and magnetic spherical harmonics. But we'll do it in stages to avoid overengineering things (which we tend to do). I'm personally going to need this for research soon so there is a good chance that it will happen relatively quickly (famous last words). I don't think it's likely we'll implement really general spherical harmonics but only the gravity and magnetics specific part.

The NOAA "python" code

It's not really Python code. It's a wrapper around a fortran or C code that does all the actual calculation.

Packaging coefficients vs downloading from NOAA or somewhere else

I like @MarkWieczorek's suggestion of duplicating the file somewhere we control and downloading with Pooch. I think I'd commit it to the repository instead of using Zenodo to avoid spamming their API any more than we already do.

Including scaling in the Plm calculation

Alright, I'll definitely have a look at that. I need to read the paper more carefully. The SHTools implementation is a huge help, thanks @MarkWieczorek! In #505 I'm calculating without it and can only go up to degree 600. That's plenty for what I want out of this right now but if it's not too hard to add I'll make the effort.

SHTools routine to read the IGRF

That does help, thanks. I wrote one for a class I'm teaching and I use Python's datetime, which is very good. The file itself isn't too hard to parse.

Rotating vectors to be tangential to the ellipsoid in Boule

Yes, definitely! It's a single calculation so I'll do it here for now but it would be great to have in Boule.

Class versus function

We were brainstorming this at the Fatiando call earlier today and came with this usage example for a class (which is our preferred option):

# Start with this class
class IGRF():
    def __init__(self, date, ellipsoid=boule.WGS84):
        # Fetch and read the coeffs from the file
        # Interpolate to the given date
        # Define g and h, normalization, reference radius

    def predict(coordinates, field="b"):
        # Predict on the given geodetic coordinates
        # Choose what to predict: V, bn, be, bu

    def grid(coordinates, field="b", ...):
        # Generate an xarray grid at the specified coordinates
        # made with verde.grid_coordinates

# Usage
be, bn, bu = IGRF(datetime(1990, 3, 14)).predict((lon, lat, height), field="b")
grid = IGRF(datetime(1990, 3, 14)).grid(grid_coords, field="b")
# Plot an intensity grid
np.sqrt(grid.bu**2 + grid.be**2 + grid.bn**2).plot()

# Later, could be generalized to this (and then IGRF inherits from it)
class MagneticSphericalHarmonics():
    def __init__(
        self, 
        max_degree, 
        reference_radius, 
        normalization="schmidt", 
        damping=None,    # Regularization for the inversion
        ellipsoid=None,  # None means spherical coordinates
        g=None, h=None,  # Can initialize the model with coefficients
)

# Usage for fitting
sh = MagneticSphericalHarmonics(max_degree=100, reference_radius=bl.WGS84.mean_radius)
coordinates = (longitude, latitude_sph, radius)
sh.fit(b_e=(coordinates, b_e_data, b_e_weights), b_n=(coordinates, b_n_data, b_n_weights))
sh.predict(..., field="b")
sh.grid(..., field="b")

What do you all think? Is there anything we're missing?

[santisoler]

Thanks @leouieda for capturing all these ideas.

To that draft, I'd just add that the fit method for the MagneticSphericalHarmonics, that would look something like:

    def fit(self, be=None, bn=None, bu=None):
        # Fit the sph coefficients with these data. At least one of the arguments
        # should not be None. They would be in the form of `(coordinates, data, weights)`.
        ...
        return self

This design for the fit method was originally thought for the equivalent sources classes.

[ckohnke]

Not sure if this is useful at this point, but I recently found myself going down the IGRF rabbit hole and found this repo with a decent overview of a lot of the current projects to to calculate the IGRF in the readme. Unfortunately there seems to be a lack of modern ways to efficiently compute the reference field over a grid-like structure, which this would hopefully solve.

[leouieda]

@ckohnke thanks for the link! I've seen several of those and many of them are wrappers of the same Fortran code from NOAA. But it's nice to see what's available.