load_earth_magnetic_anomaly: Add mag4km parameter to support earth_mag4km dataset

pygmt/datasets/earth_magnetic_anomaly.py

@@ -9,19 +9,24 @@
 
 
 @kwargs_to_strings(region="sequence")
-def load_earth_magnetic_anomaly(resolution="01d", region=None, registration=None):
+def load_earth_magnetic_anomaly(
+    resolution="01d", region=None, registration=None, mag4km=False
+):
     r"""
     Load an Earth magnetic anomaly grid in various resolutions.
 
     The grids are downloaded to a user data directory
-    (usually ``~/.gmt/server/earth/earth_mag/``) the first time you invoke
+    (usually ``~/.gmt/server/earth/earth_mag/`` or
+    ``~/.gmt/server/earth/earth_mag4km/``) the first time you invoke
     this function. Afterwards, it will load the grid from the data directory.
     So you'll need an internet connection the first time around.
 
     These grids can also be accessed by passing in the file name
-    **@earth_mag**\_\ *res*\[_\ *reg*] to any grid plotting/processing
-    function. *res* is the grid resolution (see below), and *reg* is grid
-    registration type (**p** for pixel registration or **g** for gridline
+    **@**\ *earth_mag_type*\_\ *res*\[_\ *reg*] to any grid plotting/processing
+    function. *earth_mag_type* is the GMT name
+    for the dataset. The available options are **earth_mag** and
+    **earth_mag4km**. *res* is the grid resolution (see below), and *reg* is
+    grid registration type (**p** for pixel registration or **g** for gridline
     registration).
 
     Refer to :gmt-datasets:`earth-mag.html` for more details.
@@ -46,6 +51,13 @@ def load_earth_magnetic_anomaly(resolution="01d", region=None, registration=None
         a pixel-registered grid is returned unless only the
         gridline-registered grid is available.
 
+    mag4km : bool
+        Choose the data version to use. The default is ``False``, which is
+        observed at sea level over oceanic regions and has no data over land.
+        Set ``mag4km`` to ``True`` to use a version where all observations
+        are relative to an altitude of 4 km above the geoid and include data
+        over land.
+
     Returns
     -------
     grid : :class:`xarray.DataArray`
@@ -58,7 +70,7 @@ def load_earth_magnetic_anomaly(resolution="01d", region=None, registration=None
     Earth magnetic anomaly with resolutions of 5 arc-minutes or higher,
     which are stored as smaller tiles.
     """
-    dataset_prefix = "earth_mag_"
+    dataset_prefix = "earth_mag4km_" if mag4km is True else "earth_mag_"
     grid = _load_remote_dataset(
         dataset_name="earth_magnetic_anomaly",
         dataset_prefix=dataset_prefix,

pygmt/helpers/testing.py

@@ -184,8 +184,10 @@ def download_test_data():
         "@earth_geoid_01d_g.grd",
         "@S90W180.earth_geoid_05m_g.nc",  # Specific grid for 05m test
         # Earth magnetic anomaly grids
-        "@earth_mag_01d_g.grd",
+        "@earth_mag_01d_g",
         "@S90W180.earth_mag_05m_g.nc",  # Specific grid for 05m test
+        "@earth_mag4km_01d_g",
+        "@S90W180.earth_mag4km_05m_g.nc",  # Specific grid for 05m test
         # Other cache files
         "@capitals.gmt",
         "@earth_relief_20m_holes.grd",

pygmt/tests/test_datasets_earth_magnetic_anomaly.py

@@ -80,3 +80,69 @@ def test_earth_mag_05m_without_region():
     """
     with pytest.raises(GMTInvalidInput):
         load_earth_magnetic_anomaly("05m")
+
+
+def test_earth_mag_incorrect_resolution_registration():
+    """
+    Test that an error is raised when trying to load a grid registration with
+    an unavailable resolution.
+    """
+    with pytest.raises(GMTInvalidInput):
+        load_earth_magnetic_anomaly(
+            resolution="02m", region=[0, 1, 3, 5], registration="gridline", mag4km=False
+        )
+
+
+def test_earth_mag4km_01d():
+    """
+    Test some properties of the magnetic anomaly 4km 01d data.
+    """
+    data = load_earth_magnetic_anomaly(
+        resolution="01d", registration="gridline", mag4km=True
+    )
+    assert data.name == "magnetic_anomaly"
+    assert data.attrs["long_name"] == "Earth magnetic anomaly"
+    assert data.attrs["units"] == "nT"
+    assert data.attrs["horizontal_datum"] == "WGS84"
+    assert data.shape == (181, 361)
+    npt.assert_allclose(data.lat, np.arange(-90, 91, 1))
+    npt.assert_allclose(data.lon, np.arange(-180, 181, 1))
+    npt.assert_allclose(data.min(), -799.19995)
+    npt.assert_allclose(data.max(), 3226.4)
+
+
+def test_earth_mag4km_01d_with_region():
+    """
+    Test loading low-resolution earth magnetic anomaly 4km 01d with 'region'.
+    """
+    data = load_earth_magnetic_anomaly(
+        resolution="01d",
+        region=[-10, 10, -5, 5],
+        registration="gridline",
+        mag4km=True,
+    )
+    assert data.shape == (11, 21)
+    npt.assert_allclose(data.lat, np.arange(-5, 6, 1))
+    npt.assert_allclose(data.lon, np.arange(-10, 11, 1))
+    npt.assert_allclose(data.min(), -153.19995)
+    npt.assert_allclose(data.max(), 113.59985)
+
+
+def test_earth_mag4km_05m_with_region():
+    """
+    Test loading a subregion of high-resolution earth magnetic anomaly 4km
+    data.
+    """
+    data = load_earth_magnetic_anomaly(
+        resolution="05m",
+        region=[-115, -112, 4, 6],
+        registration="gridline",
+        mag4km=True,
+    )
+    assert data.shape == (25, 37)
+    assert data.lat.min() == 4
+    assert data.lat.max() == 6
+    assert data.lon.min() == -115
+    assert data.lon.max() == -112
+    npt.assert_allclose(data.min(), -128.40015)
+    npt.assert_allclose(data.max(), 76.80005)