Remote datasets: Add "load_earth_mean_dynamic_topography" to load "CN…

…ES Earth Mean Dynamic Topography" dataset (#3718)

Co-authored-by: Dongdong Tian <[email protected]>

doc/api/index.rst

@@ -238,6 +238,7 @@ and store them in GMT's user data directory.
     datasets.load_earth_geoid
     datasets.load_earth_magnetic_anomaly
     datasets.load_earth_mask
+    datasets.load_earth_mean_dynamic_topography
     datasets.load_earth_relief
     datasets.load_earth_vertical_gravity_gradient
     datasets.load_mars_relief

pygmt/datasets/__init__.py

@@ -11,6 +11,9 @@
 from pygmt.datasets.earth_geoid import load_earth_geoid
 from pygmt.datasets.earth_magnetic_anomaly import load_earth_magnetic_anomaly
 from pygmt.datasets.earth_mask import load_earth_mask
+from pygmt.datasets.earth_mean_dynamic_topography import (
+    load_earth_mean_dynamic_topography,
+)
 from pygmt.datasets.earth_night import load_black_marble
 from pygmt.datasets.earth_relief import load_earth_relief
 from pygmt.datasets.earth_vertical_gravity_gradient import (

pygmt/datasets/earth_mean_dynamic_topography.py

@@ -0,0 +1,103 @@
+"""
+Function to download the CNES Earth mean dynamic topography dataset from the GMT data
+server, and load as :class:`xarray.DataArray`.
+
+The grids are available in various resolutions.
+"""
+
+from collections.abc import Sequence
+from typing import Literal
+
+import xarray as xr
+from pygmt.datasets.load_remote_dataset import _load_remote_dataset
+
+__doctest_skip__ = ["load_earth_mean_dynamic_topography"]
+
+
+def load_earth_mean_dynamic_topography(
+    resolution: Literal["01d", "30m", "20m", "15m", "10m", "07m"] = "01d",
+    region: Sequence[float] | str | None = None,
+    registration: Literal["gridline", "pixel"] = "gridline",
+) -> xr.DataArray:
+    r"""
+    Load the CNES Earth mean dynamic topography dataset in various resolutions.
+
+    .. figure:: https://www.generic-mapping-tools.org/remote-datasets/_images/GMT_earth_mdt.jpg
+       :width: 80 %
+       :align: center
+
+       CNES Earth mean dynamic topography dataset.
+
+    The grids are downloaded to a user data directory (usually
+    ``~/.gmt/server/earth/earth_mdt/``) 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_mdt**\_\ *res*\[_\ *reg*] to any grid processing function or plotting
+    method. *res* is the grid resolution (see below), and *reg* is the grid registration
+    type (**p** for pixel registration or **g** for gridline registration).
+
+    The default color palette table (CPT) for this dataset is *@earth_mdt.cpt*. It's
+    implicitly used when passing in the file name of the dataset to any grid plotting
+    method if no CPT is explicitly specified. When the dataset is loaded and plotted
+    as an :class:`xarray.DataArray` object, the default CPT is ignored, and GMT's
+    default CPT (*turbo*) is used. To use the dataset-specific CPT, you need to
+    explicitly set ``cmap="@earth_mdt.cpt"``.
+
+    Refer to :gmt-datasets:`earth-mdt.html` for more details about available datasets,
+    including version information and references.
+
+    Parameters
+    ----------
+    resolution
+        The grid resolution. The suffix ``d`` and ``m`` stand for arc-degrees and
+        arc-minutes. Note that ``"07m"`` refers to a resolution of 7.5 arc-minutes.
+    region
+        The subregion of the grid to load, in the form of a sequence [*xmin*, *xmax*,
+        *ymin*, *ymax*] or an ISO country code.
+    registration
+        Grid registration type. Either ``"pixel"`` for pixel registration or
+        ``"gridline"`` for gridline registration.
+
+    Returns
+    -------
+    grid
+        The CNES Earth mean dynamic topography grid. Coordinates are latitude and
+        longitude in degrees. Values are in meters.
+
+    Note
+    ----
+    The registration and coordinate system type of the returned
+    :class:`xarray.DataArray` grid can be accessed via the GMT accessors (i.e.,
+    ``grid.gmt.registration`` and ``grid.gmt.gtype`` respectively). However, these
+    properties may be lost after specific grid operations (such as slicing) and will
+    need to be manually set before passing the grid to any PyGMT data processing or
+    plotting functions. Refer to :class:`pygmt.GMTDataArrayAccessor` for detailed
+    explanations and workarounds.
+
+    Examples
+    --------
+
+    >>> from pygmt.datasets import load_earth_mean_dynamic_topography
+    >>> # load the default grid (gridline-registered 1 arc-degree grid)
+    >>> grid = load_earth_mean_dynamic_topography()
+    >>> # load the 30 arc-minutes grid with "gridline" registration
+    >>> grid = load_earth_mean_dynamic_topography(
+            resolution="30m", registration="gridline"
+    ...)
+    >>> # load high-resolution (5 arc-minutes) grid for a specific region
+    >>> grid = load_earth_mean_dynamic_topography(
+    ...     resolution="05m",
+    ...     region=[120, 160, 30, 60],
+    ...     registration="gridline",
+    ... )
+    """
+    grid = _load_remote_dataset(
+        name="earth_mdt",
+        prefix="earth_mdt",
+        resolution=resolution,
+        region=region,
+        registration=registration,
+    )
+    return grid

pygmt/datasets/load_remote_dataset.py

@@ -246,6 +246,19 @@ class GMTRemoteDataset(NamedTuple):
             "30s": Resolution("30s", registrations=["pixel"]),
         },
     ),
+    "earth_mdt": GMTRemoteDataset(
+        description="CNES Earth mean dynamic topography",
+        units="meters",
+        extra_attributes={"horizontal_datum": "WGS84"},
+        resolutions={
+            "01d": Resolution("01d"),
+            "30m": Resolution("30m"),
+            "20m": Resolution("20m"),
+            "15m": Resolution("15m"),
+            "10m": Resolution("10m"),
+            "07m": Resolution("07m", registrations=["gridline"]),
+        },
+    ),
     "earth_vgg": GMTRemoteDataset(
         description="IGPP Earth vertical gravity gradient",
         units="Eotvos",

pygmt/helpers/caching.py

@@ -23,6 +23,8 @@ def cache_data():
         "@earth_mag_01d_g",
         "@earth_mag4km_01d_g",
         "@earth_mask_01d_g",
+        "@earth_mdt_01d_g",
+        "@earth_mdt_07m_g",
         "@earth_night_01d",
         "@earth_relief_01d_g",
         "@earth_relief_01d_p",

pygmt/tests/test_datasets_mean_dynamic_topography.py

@@ -0,0 +1,53 @@
+"""
+Test basic functionality for loading Earth mean dynamic topography datasets.
+"""
+
+import numpy as np
+import numpy.testing as npt
+from pygmt.datasets import load_earth_mean_dynamic_topography
+
+
+def test_earth_mdt_01d():
+    """
+    Test some properties of the Earth mean dynamic topography 01d data.
+    """
+    data = load_earth_mean_dynamic_topography(resolution="01d")
+    assert data.name == "z"
+    assert data.attrs["description"] == "CNES Earth mean dynamic topography"
+    assert data.attrs["units"] == "meters"
+    assert data.attrs["horizontal_datum"] == "WGS84"
+    assert data.shape == (181, 361)
+    assert data.gmt.registration == 0
+    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(), -1.4668, atol=0.0001)
+    npt.assert_allclose(data.max(), 1.7151, atol=0.0001)
+
+
+def test_earth_mdt_01d_with_region():
+    """
+    Test loading low-resolution Earth mean dynamic topography with "region".
+    """
+    data = load_earth_mean_dynamic_topography(resolution="01d", region=[-10, 10, -5, 5])
+    assert data.shape == (11, 21)
+    assert data.gmt.registration == 0
+    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(), 0.346, atol=0.0001)
+    npt.assert_allclose(data.max(), 0.4839, atol=0.0001)
+
+
+def test_earth_mdt_07m_default_registration():
+    """
+    Test that the grid returned by default for the 7 arc-minutes resolution has a
+    "gridline" registration.
+    """
+    data = load_earth_mean_dynamic_topography(resolution="07m", region=[-10, -9, 3, 5])
+    assert data.shape == (17, 9)
+    assert data.gmt.registration == 0
+    assert data.coords["lat"].data.min() == 3.0
+    assert data.coords["lat"].data.max() == 5.0
+    assert data.coords["lon"].data.min() == -10.0
+    assert data.coords["lon"].data.max() == -9.0
+    npt.assert_allclose(data.min(), 0.4138, atol=0.0001)
+    npt.assert_allclose(data.max(), 0.4302, atol=0.0001)