load_moon_relief: Load moon relief dataset in various resolutions and…

… registrations

doc/api/index.rst

@@ -226,6 +226,7 @@ and store them in GMT's user data directory.
     datasets.load_earth_mask
     datasets.load_earth_relief
     datasets.load_earth_vertical_gravity_gradient
+    datasets.load_moon_relief
     datasets.load_sample_data
 
 In addition, there is also a special function to load XYZ tile maps via

pygmt/datasets/__init__.py

@@ -13,5 +13,6 @@
 from pygmt.datasets.earth_vertical_gravity_gradient import (
     load_earth_vertical_gravity_gradient,
 )
+from pygmt.datasets.moon_relief import load_moon_relief
 from pygmt.datasets.samples import list_sample_data, load_sample_data
 from pygmt.datasets.tile_map import load_tile_map

pygmt/datasets/load_remote_dataset.py

@@ -229,6 +229,29 @@ class GMTRemoteDataset(NamedTuple):
             "03m": Resolution(["gridline"], True),
         },
     ),
+    "moon_relief": GMTRemoteDataset(
+        title="Moon relief",
+        name="moon relief",
+        long_name="Moon relief",
+        units="meters",
+        extra_attributes={},
+        resolutions={
+            "01d": Resolution(["gridline", "pixel"], False),
+            "30m": Resolution(["gridline", "pixel"], False),
+            "20m": Resolution(["gridline", "pixel"], False),
+            "15m": Resolution(["gridline", "pixel"], False),
+            "10m": Resolution(["gridline", "pixel"], False),
+            "06m": Resolution(["gridline", "pixel"], False),
+            "05m": Resolution(["gridline", "pixel"], True),
+            "04m": Resolution(["gridline", "pixel"], True),
+            "03m": Resolution(["gridline", "pixel"], True),
+            "02m": Resolution(["gridline", "pixel"], True),
+            "01m": Resolution(["gridline", "pixel"], True),
+            "30s": Resolution(["gridline", "pixel"], True),
+            "15s": Resolution(["gridline", "pixel"], True),
+            "14s": Resolution(["pixel"], True),
+        },
+    ),
 }
 
 

pygmt/datasets/moon_relief.py

@@ -0,0 +1,110 @@
+"""
+Function to download the Moon relief dataset from the GMT data server, and load
+as :class:`xarray.DataArray`.
+
+The grids are available in various resolutions.
+"""
+from typing import Literal
+
+from pygmt.datasets.load_remote_dataset import _load_remote_dataset
+from pygmt.helpers import kwargs_to_strings
+
+__doctest_skip__ = ["load_moon_relief"]
+
+
+@kwargs_to_strings(region="sequence")
+def load_moon_relief(
+    resolution="01d",
+    region=None,
+    registration: Literal["gridline", "pixel", None] = None,
+):
+    r"""
+    Load the Moon relief dataset in various resolutions.
+
+    .. figure:: https://www.generic-mapping-tools.org/remote-datasets/_images/GMT_moon_relief.jpg
+       :width: 80%
+       :align: center
+
+       Moon relief dataset.
+
+    The grids are downloaded to a user data directory
+    (usually ``~/.gmt/server/moon/moon_relief/``) 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
+    **@moon_relief**\_\ *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
+    *@moon_relief.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="@moon_relief.cpt"``.
+
+    Refer to :gmt-datasets:`moon-relief.html` for more details about available
+    datasets, including version information and references.
+
+    Parameters
+    ----------
+    resolution : str
+        The grid resolution. The suffix ``d``, ``m`` and ``s`` stand for
+        arc-degrees, arc-minutes and arc-seconds. It can be ``"01d"``,
+        ``"30m"``, ``"20m"``, ``"15m"``, ``"10m"``, ``"06m"``, ``"05m"``,
+        ``"04m"``, ``"03m"``, ``"02m"``, ``"01m"``, ``"30s"``, ``"15s"``, and
+        ``"14s"``.
+    region : str or list
+        The subregion of the grid to load, in the form of a list
+        [*xmin*, *xmax*, *ymin*, *ymax*] or a string *xmin/xmax/ymin/ymax*.
+        Required for grids with resolutions higher than 5 arc-minutes
+        (i.e., ``"05m"``).
+    registration
+        Grid registration type. Either ``"pixel"`` for pixel registration or
+        ``"gridline"`` for gridline registration. Default is ``None`` , means
+        ``"gridline"`` for all resolutions except for ``"14s"`` which is
+        ``"pixel"`` only.
+
+    Returns
+    -------
+    grid : :class:`xarray.DataArray`
+        The Moon relief grid. Coordinates are latitude and longitude in
+        degrees. Relief is 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_moon_relief
+    >>> # load the default grid (gridline-registered 1 arc-degree grid)
+    >>> grid = load_moon_relief()
+    >>> # load the 30 arc-minutes grid with "gridline" registration
+    >>> grid = load_moon_relief(resolution="30m", registration="gridline")
+    >>> # load high-resolution (5 arc-minutes) grid for a specific region
+    >>> grid = load_moon_relief(
+    ...     resolution="05m",
+    ...     region=[120, 160, 30, 60],
+    ...     registration="gridline",
+    ... )
+    """
+    grid = _load_remote_dataset(
+        dataset_name="moon_relief",
+        dataset_prefix="moon_relief_",
+        resolution=resolution,
+        region=region,
+        registration=registration,
+    )
+    return grid

pygmt/helpers/testing.py

@@ -200,6 +200,9 @@ def download_test_data():
         # Earth WDMAM grids
         "@earth_wdmam_01d_g",
         "@S90E000.earth_wdmam_03m_g.nc",  # Specific grid for 03m test
+        # Moon relief grids
+        "@moon_relief_01d_g",
+        "@N00W030.moon_relief_01m_g.nc",  # Specific grid for 01m test
         # Other cache files
         "@capitals.gmt",
         "@earth_relief_20m_holes.grd",

pygmt/tests/test_datasets_moon_relief.py

@@ -0,0 +1,51 @@
+"""
+Test basic functionality for loading Moon relief datasets.
+"""
+import numpy as np
+import numpy.testing as npt
+from pygmt.datasets import load_moon_relief
+
+
+def test_moon_relief_01d():
+    """
+    Test some properties of the Moon relief 01d data.
+    """
+    data = load_moon_relief(resolution="01d")
+    assert data.name == "moon relief"
+    assert data.attrs["units"] == "meters"
+    assert data.attrs["long_name"] == "Moon relief"
+    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(), -7669.0, atol=0.5)
+    npt.assert_allclose(data.max(), 8507.5, atol=0.5)
+
+
+def test_moon_relief_01d_with_region():
+    """
+    Test loading low-resolution Moon relief with 'region'.
+    """
+    data = load_moon_relief(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(), -1122.0, atol=0.5)
+    npt.assert_allclose(data.max(), 943.0, atol=0.5)
+
+
+def test_moon_relief_01m_default_registration():
+    """
+    Test that the grid returned by default for the 1 arc-minute resolution has
+    a "gridline" registration.
+    """
+    data = load_moon_relief(resolution="01m", region=[-10, -9, 3, 5])
+    assert data.shape == (121, 61)
+    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(), -1840.5, atol=0.5)
+    npt.assert_allclose(data.max(), -235.0, atol=0.5)