Trying to make test data

Using the BGS API to make some grids of the IGRF to test against their
results.

Author: leouieda

bx.npy

@@ -1,6 +1,13 @@
+# Copyright (c) 2018 The Harmonica Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
 """
 Calculation of the IGRF magnetic field.
 """
+
 import pathlib
 
 import boule
@@ -14,8 +21,7 @@
 
 
 def fetch_igrf13():
-    """
-    """
+    """ """
     path = pooch.retrieve(
         url="doi:10.5281/zenodo.11269410/igrf13coeffs.txt",
         path=pooch.os_cache("harmonica"),
@@ -40,9 +46,9 @@ def interpolate_coefficients(date, g, h, years):
     return g_date, h_date
 
 
-class IGRF13():
+class IGRF14:
     """
-    13th generation of the International Geomagnetic Reference Field
+    14th generation of the International Geomagnetic Reference Field
     """
 
     def __init__(self, date, ellipsoid=boule.WGS84):
@@ -82,12 +88,16 @@ def predict(self, coordinates, field="b"):
         b_north = np.zeros(n_data)
         b_up = np.zeros(n_data)
         g, h = self.coefficients
-        _evaluate_igrf(longitude, colatitude, radius, g, h, self.max_degree, b_east, b_north, b_up)
+        _evaluate_igrf(
+            longitude, colatitude, radius, g, h, self.max_degree, b_east, b_north, b_up
+        )
         return b_east, b_north, b_up
 
 
 @numba.jit(parallel=True, nopython=True)
-def _evaluate_igrf(longitude, colatitude, radius, g, h, max_degree, b_east, b_north, b_up):
+def _evaluate_igrf(
+    longitude, colatitude, radius, g, h, max_degree, b_east, b_north, b_up
+):
     n_data = longitude.size
     for i in numba.prange(n_data):
-        plm =
+        plm = 1

by.npy

@@ -0,0 +1,72 @@
+# Copyright (c) 2018 The Harmonica Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
+"""
+Script to use the BGS API to get values of the IGRF for testing purposes.
+
+"""
+
+import datetime
+import json
+
+import numpy as np
+import requests
+import verde as vd
+import xarray as xr
+
+coordinates = vd.grid_coordinates((0, 359, -90, 90), spacing=10, extra_coords=0)
+
+dates = ["1930-04-20", "1986-07-29", "2024-01-10", "2029-10-30"]
+
+shape = (len(dates), *coordinates[0].shape)
+bx = np.empty(shape)
+by = np.empty(shape)
+bz = np.empty(shape)
+
+url = (
+    "https://geomag.bgs.ac.uk/web_service/GMModels/igrf/14/?"
+    "latitude={latitude}&longitude={longitude}&altitude={altitude}&date={date}&"
+    "format=json"
+)
+
+for k in range(shape[0]):
+    for i in range(shape[1]):
+        for j in range(shape[2]):
+            response = requests.get(
+                url.format(
+                    longitude=coordinates[0][i, j],
+                    latitude=coordinates[1][i, j],
+                    altitude=coordinates[2][i, j],
+                    date=dates[k],
+                )
+            )
+            response.raise_for_status()
+            results = json.loads(response.text)
+            field = results["geomagnetic-field-model-result"]["field-value"]
+            bx[k, i, j] = field["east-intensity"]["value"]
+            by[k, i, j] = field["north-intensity"]["value"]
+            # The BGS gives vertical downward
+            bz[k, i, j] = -field["vertical-intensity"]["value"]
+
+np.save("bx.npy", bx)
+np.save("by.npy", by)
+np.save("bz.npy", bz)
+
+bx = np.load("bx.npy")
+by = np.load("by.npy")
+bz = np.load("bz.npy")
+
+time = [datetime.datetime.fromisoformat(d) for d in dates]
+dims = ("time", "latitude", "longitude")
+grid = xr.Dataset(
+    {"bx": (dims, bx), "by": (dims, by), "bz": (dims, bz)},
+    coords={
+        "time": time,
+        "longitude": coordinates[0][0, :],
+        "latitude": coordinates[1][:, 0],
+    },
+)
+grid.to_netcdf("igrf14-bgs.nc")

bz.npy

@@ -0,0 +1,14 @@
+# Copyright (c) 2018 The Harmonica Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
+"""
+Test the IGRF results against the ones calculated by the BGS.
+"""
+import pytest
+
+
+from .._spherical_harmonics.igrf import load_igrf, fetch_igrf,
+interpolate_coefficients, IGRF14