Add some test functions to ppigrf

README.md

@@ -55,6 +55,29 @@ Be, Bn, Bu = ppigrf.igrf(lon, lat, h, dates)
 ```
 The output will have shape `(4, 2, 3)`.
 
+## How to test?
+
+It is possible to run some tests that corroborates that `ppigrf` works as
+expected. To do so you need to install [`pytest`](https://docs.pytest.org). You
+can do so through `pip`:
+
+```
+pip install pytest
+```
+
+Then you need to clone this repository and install `ppigrf`:
+
+```
+python setup.py install
+```
+
+Finally you can test the installed version of `ppigrf` with:
+
+```
+pytest -v
+```
+
+
 ## Why?
 There are lots of Python modules that can calculate IGRF values. Most are wrappers of Fortran code, which can be tricky to compile. This version is pure Python. For most applications it is still quite fast. I also prefer the ppigrf interface over the alternatives. 
 

test/data/README.md

@@ -0,0 +1,17 @@
+# Data files used for testing ppigrf
+
+## Precomputed IGRF
+
+Each one of these folders contain a set of `b_e.csv`, `b_n.csv` and `b_z.csv`
+files that host precomputed values of the IGRF field obtained through the [NCEI
+Geomagnetic
+Calculator](https://www.ngdc.noaa.gov/geomag/calculators/magcalc.shtml)
+provided by [NOAA](https://www.ngdc.noaa.gov).
+
+These values were computed in a regular grid with a spacing of 1 degree in both
+longitudinal and latitudinal directions, at a height of 5km above the WGS84
+ellipsoid. The date of each IGRF field files is specified in the folder name
+following the `YYYY-MM-DD` format.
+
+The `b_z.csv` contains the **downward** components of the magnetic vector on
+each location.

test/test_igrf.py

@@ -0,0 +1,99 @@
+"""
+Test IGRF field
+"""
+import os
+import pytest
+import numpy as np
+import numpy.testing as npt
+import pandas as pd
+from pathlib import Path
+from datetime import datetime
+
+from ppigrf import igrf
+from ppigrf.ppigrf import yearfrac_to_datetime
+
+# Define paths to test directory and test data directory
+TEST_DIR = Path(os.path.dirname(__file__))
+TEST_DATA_DIR = TEST_DIR / "data"
+
+
+def load_precomputed_igrf(date):
+    """
+    Loads the precomputed IGRF files from a given date
+
+    Available dates:
+        * 2020-01-01
+        * 2022-10-05
+
+    Parameters
+    ----------
+    date : :class:`datetime.datetime` object
+        Date of the precomputed IGRF field files that will be loaded.
+
+    Returns
+    -------
+    igrf_precomputed : :class:`pandas.Dataframe`
+        Dataframe containing the precomputed values of the IGRF on the given
+        date.
+    """
+    # Read the csv files
+    date_dir = TEST_DATA_DIR / date.strftime("%Y-%m-%d")
+    first_columns = ["date", "latitude", "longitude", "altitude_km"]
+    components = ("b_e", "b_n", "b_z")
+    dataframes = []
+    for component in components:
+        columns = first_columns + [component, component + "_sv"]
+        fname = date_dir / f"{component}.csv"
+        df = pd.read_csv(fname, skiprows=13, names=columns)
+        dataframes.append(df)
+    # Merge the dataframes
+    igrf_precomputed = pd.merge(dataframes[0], dataframes[1])
+    igrf_precomputed = pd.merge(igrf_precomputed, dataframes[-1])
+    # Convert the data in the dataframe into a datetime object
+    decimal_date = igrf_precomputed.date.values[0]
+    (date,) = yearfrac_to_datetime([decimal_date])
+    igrf_precomputed = igrf_precomputed.assign(date=date)
+    return igrf_precomputed
+
+
+class TestIGRFKnownValues:
+    """
+    Test the IGRF field against precomputed values
+    """
+
+    rtol = 1e-2  # 1% of error
+
+    @pytest.mark.parametrize(
+        "date, atol",
+        [[datetime(2020, 1, 1), 1], [datetime(2022, 10, 5), 4]],
+        ids=["2020-01-01", "2022-10-05"],
+    )
+    def test_igrf(self, date, atol):
+        """
+        Test IGRF against the precomputed values
+
+        The test on 2020-01-01 doesn't involve any interpolation on the
+        dates. The atol (in nT) has been chosen for each case to account points
+        where the component is close to zero. For the second date that involves
+        an interpolation in time the atol has been increased to account for
+        differences due to different types of dates interpolations.
+        """
+        # Get precomputed IGRF field
+        precomputed_igrf = load_precomputed_igrf(date)
+        # Overwrite the date with the one in the data file
+        # date = precomputed_igrf.date.values[0]
+        # Compute igrf using ppigrf
+        b_e, b_n, b_u = igrf(
+            precomputed_igrf.longitude,
+            precomputed_igrf.latitude,
+            precomputed_igrf.altitude_km,
+            date,
+        )
+        # Ravel the arrays
+        b_e, b_n, b_u = tuple(np.ravel(component) for component in (b_e, b_n, b_u))
+        # Check if the values are equal to the expected ones
+        npt.assert_allclose(b_e, precomputed_igrf.b_e, rtol=self.rtol, atol=atol)
+        npt.assert_allclose(b_n, precomputed_igrf.b_n, rtol=self.rtol, atol=atol)
+        npt.assert_allclose(
+            b_u, -precomputed_igrf.b_z, rtol=self.rtol, atol=atol
+        )  # invert the direction of b_z