diff --git a/pygmt/datasets/earth_magnetic_anomaly.py b/pygmt/datasets/earth_magnetic_anomaly.py
index 80a66358161..1ad9bfae223 100644
--- a/pygmt/datasets/earth_magnetic_anomaly.py
+++ b/pygmt/datasets/earth_magnetic_anomaly.py
@@ -7,6 +7,8 @@
 from pygmt.datasets.load_remote_dataset import _load_remote_dataset
 from pygmt.helpers import kwargs_to_strings
 
+__doctest_skip__ = ["load_earth_magnetic_anomaly"]
+
 
 @kwargs_to_strings(region="sequence")
 def load_earth_magnetic_anomaly(
@@ -69,6 +71,27 @@ def load_earth_magnetic_anomaly(
     The :class:`xarray.DataArray` grid doesn't support slice operation, for
     Earth magnetic anomaly with resolutions of 5 arc-minutes or higher,
     which are stored as smaller tiles.
+
+    Examples
+    --------
+
+    >>> from pygmt.datasets import load_earth_magnetic_anomaly
+    >>> # load the default grid (gridline-registered 1 arc-degree grid)
+    >>> grid = load_earth_magnetic_anomaly()
+    >>> # load the 30 arc-minutes grid with "gridline" registration
+    >>> grid = load_earth_magnetic_anomaly(
+    ...     resolution="30m", registration="gridline"
+    ... )
+    >>> # load high-resolution (5 arc-minutes) grid for a specific region
+    >>> grid = load_earth_magnetic_anomaly(
+    ...     resolution="05m",
+    ...     region=[120, 160, 30, 60],
+    ...     registration="gridline",
+    ... )
+    >>> # load the 20 arc-minutes grid of the mag4km dataset
+    >>> grid = load_earth_magnetic_anomaly(
+    ...     resolution="20m", registration="gridline", mag4km=True
+    ... )
     """
     dataset_prefix = "earth_mag4km_" if mag4km is True else "earth_mag_"
     grid = _load_remote_dataset(