Fix interpolated coefficient shape

The interpolation was always going to the maximum degree instead of
stopping at the desired maximum.

Author: leouieda

harmonica/_spherical_harmonics/igrf.py

@@ -76,7 +76,7 @@ def load_igrf(path):
     return years, coeffs["g"], coeffs["h"], coeffs["g_sv"], coeffs["h_sv"]
 
 
-def interpolate_coefficients(date, years, g, h, g_sv, h_sv):
+def interpolate_coefficients(date, max_degree, years, g, h, g_sv, h_sv):
     """
     Interpolate the Gauss coefficients to the given date.
 
@@ -87,6 +87,8 @@ def interpolate_coefficients(date, years, g, h, g_sv, h_sv):
     ----------
     date : :class:`datetime.datetime`
         The date and time at which to interpolate the coefficients.
+    max_degree : int
+        The maximum degree that will be interpolated.
     years : 1d-array
         The years for the knot points (epochs) of the time interpolation.
     g, h : 3d-arrays
@@ -106,7 +108,8 @@ def interpolate_coefficients(date, years, g, h, g_sv, h_sv):
     g_date, h_date : 2d-arrays
         The interpolated Gauss coefficients. The first dimension is the degree
         n and the second is the order m. At m > n, the coefficients are
-        assigned zero. Units are nT.
+        assigned zero. The shape of the arrays will be ``(max_degree + 1,
+        max_degree + 1)``. Units are nT.
     """
     if date.year < years[0] or date.year >= years[-1] + 5:
         message = (
@@ -124,9 +127,9 @@ def interpolate_coefficients(date, years, g, h, g_sv, h_sv):
         - datetime.datetime(year=years[index], month=1, day=1)
     ).total_seconds()
     year_in_seconds = 365.25 * 24 * 60 * 60
-    g_date = np.zeros(g.shape[1:])
-    h_date = np.zeros(h.shape[1:])
-    for n in range(1, g.shape[1]):
+    g_date = np.zeros((max_degree + 1, max_degree + 1))
+    h_date = np.zeros((max_degree + 1, max_degree + 1))
+    for n in range(1, max_degree + 1):
         for m in range(n + 1):
             if index >= years.size - 1:
                 variation_g = g_sv[n, m] / year_in_seconds
@@ -294,7 +297,7 @@ def __init__(
     def coefficients(self):
         if self._coefficients is None:
             self._coefficients = interpolate_coefficients(
-                self.date, *load_igrf(self._fetch_coefficient_file())
+                self.date, self.max_degree, *load_igrf(self._fetch_coefficient_file())
             )
         return self._coefficients
 

harmonica/tests/test_igrf.py

@@ -163,6 +163,7 @@ def test_interpolate_coefficients():
     for i, year in enumerate(years):
         g_date, h_date = interpolate_coefficients(
             datetime.datetime(year, month=1, day=1, hour=0, minute=1, second=0),
+            g.shape[1] - 1,
             years,
             g,
             h,
@@ -175,6 +176,23 @@ def test_interpolate_coefficients():
                 npt.assert_allclose(h_date[n, m], h[i, n, m], atol=0.001)
 
 
+def test_interpolate_coefficients_max_degree():
+    "Check that the returned coefficients stop at the max degree"
+    max_degree = 13
+    years, g, h, g_sv, h_sv = load_igrf(IGRF14("2020-02-10")._fetch_coefficient_file())
+    g_date, h_date = interpolate_coefficients(
+        datetime.datetime(year=2023, month=1, day=20, hour=0, minute=1, second=0),
+        max_degree,
+        years,
+        g,
+        h,
+        g_sv,
+        h_sv,
+    )
+    assert g_date.shape == h_date.shape
+    assert g_date.shape == (max_degree + 1, max_degree + 1)
+
+
 @pytest.mark.parametrize(
     "date",
     [datetime.datetime(1800, month=1, day=1), datetime.datetime(2030, month=1, day=1)],
@@ -184,6 +202,7 @@ def test_interpolate_coefficients_invalid_date(date):
     with pytest.raises(ValueError, match="Invalid date"):
         g_date, h_date = interpolate_coefficients(
             date,
+            None,
             list(range(1900, 2030, 5)),
             None,
             None,