Add Earth deflection dataset

Author: yvonnefroehlich

doc/api/index.rst

@@ -235,6 +235,7 @@ and store them in GMT's user data directory.
     datasets.load_blue_marble
     datasets.load_earth_age
     datasets.load_earth_dist
+    datasets.load_earth_east_west_deflection
     datasets.load_earth_free_air_anomaly
     datasets.load_earth_geoid
     datasets.load_earth_magnetic_anomaly
@@ -248,6 +249,7 @@ and store them in GMT's user data directory.
     datasets.load_moon_relief
     datasets.load_pluto_relief
     datasets.load_venus_relief
+    datasets.load_earth_south_north_deflection
     datasets.load_sample_data
 
 In addition, there is also a special function to load XYZ tile maps via

pygmt/datasets/__init__.py

@@ -7,6 +7,7 @@
 from pygmt.datasets.earth_age import load_earth_age
 from pygmt.datasets.earth_day import load_blue_marble
 from pygmt.datasets.earth_dist import load_earth_dist
+from pygmt.datasets.earth_east_west_deflection import load_earth_east_west_deflection
 from pygmt.datasets.earth_free_air_anomaly import load_earth_free_air_anomaly
 from pygmt.datasets.earth_geoid import load_earth_geoid
 from pygmt.datasets.earth_magnetic_anomaly import load_earth_magnetic_anomaly
@@ -17,6 +18,9 @@
 from pygmt.datasets.earth_mean_sea_surface import load_earth_mean_sea_surface
 from pygmt.datasets.earth_night import load_black_marble
 from pygmt.datasets.earth_relief import load_earth_relief
+from pygmt.datasets.earth_south_north_deflection import (
+    load_earth_south_north_deflection,
+)
 from pygmt.datasets.earth_vertical_gravity_gradient import (
     load_earth_vertical_gravity_gradient,
 )

pygmt/datasets/earth_deflection.py

@@ -0,0 +1,117 @@
+"""
+Function to download the IGPP Earth east-west and south-north deflection datasets 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_deflection"]
+
+
+def load_earth_deflection(
+    resolution: Literal[
+        "01d", "30m", "20m", "15m", "10m", "06m", "05m", "04m", "03m", "02m", "01m"
+    ] = "01d",
+    region: Sequence[float] | str | None = None,
+    registration: Literal["gridline", "pixel", None] = None,
+    direction: Literal["edefl", "ndefl"] = "edefl",
+) -> xr.DataArray:
+    r"""
+    Load the IGPP Earth east-west and south-north deflections datasets in various
+    resolutions.
+
+    .. list-table::
+       :widths: 50 50
+       :header-rows: 1
+
+       * - IGPP Earth east-west deflection
+         - IGPP Earth south-north deflection
+       * - .. figure:: https://www.generic-mapping-tools.org/remote-datasets/_images/GMT_earth_edefl.jpg
+         - .. figure:: https://www.generic-mapping-tools.org/remote-datasets/_images/GMT_earth_ndefl.jpg
+
+    The grids are downloaded to a user data directory (usually
+    ``~/.gmt/server/earth/earth_edefl/`` and``~/.gmt/server/earth/earth_ndefl/`` 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_defl_type*\_\ *res*\[_\ *reg*] to any grid processing function or
+    plotting method. *earth_defl_type* is the GMT name for the dataset. The available
+    options are **earth_edefl** and **earth_ndefl**. *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 (CPTs) for this dataset is *@earth_defl.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_defl.cpt"``.
+
+    Refer to :gmt-datasets:`earth-edefl.html` and :gmt-datasets:`earth-ndefl.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.
+    region
+        The subregion of the grid to load, in the form of a sequence [*xmin*, *xmax*,
+        *ymin*, *ymax*] or an ISO country code. 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 ``"01m"`` which is ``"pixel"`` only.
+    direction
+        By default, the east-west deflection (``direction="edefl"``) is returned, set
+        ``direction="ndefl"`` to return the south-north deflection.
+
+    Returns
+    -------
+    grid
+        The Earth east-west or south-north deflection grid. Coordinates are latitude
+        and longitude in degrees. Units are in micro-radians.
+
+    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_deflection
+    >>> # load the default grid for east-west deflection (gridline-registered 1
+    >>> # arc-degree grid)
+    >>> grid = load_earth_deflection()
+    >>> # load the default grid for south-north deflection
+    >>> grid = load_earth_deflection(direction="ndefl")
+    >>> # load the 30 arc-minutes grid with "gridline" registration
+    >>> grid = load_earth_deflection(resolution="30m", registration="gridline")
+    >>> # load high-resolution (5 arc-minutes) grid for a specific region
+    >>> grid = load_earth_deflection(
+    ...     resolution="05m", region=[120, 160, 30, 60], registration="gridline"
+    ... )
+    """
+    prefix = "earth_ndefl" if direction == "ndefl" else "earth_edefl"
+    grid = _load_remote_dataset(
+        name=prefix,
+        prefix=prefix,
+        resolution=resolution,
+        region=region,
+        registration=registration,
+    )
+    return grid

pygmt/datasets/load_remote_dataset.py

@@ -110,6 +110,24 @@ class GMTRemoteDataset(NamedTuple):
             "01m": Resolution("01m", registrations=["gridline"], tiled=True),
         },
     ),
+    "earth_edefl": GMTRemoteDataset(
+        description="IGPP Earth east-west deflections",
+        units="micro-radians",
+        extra_attributes={"horizontal_datum": "WGS84"},
+        resolutions={
+            "01d": Resolution("01d"),
+            "30m": Resolution("30m"),
+            "20m": Resolution("20m"),
+            "15m": Resolution("15m"),
+            "10m": Resolution("10m"),
+            "06m": Resolution("06m"),
+            "05m": Resolution("05m", tiled=True),
+            "04m": Resolution("04m", tiled=True),
+            "03m": Resolution("03m", tiled=True),
+            "02m": Resolution("02m", tiled=True),
+            "01m": Resolution("01m", registrations=["pixel"], tiled=True),
+        },
+    ),
     "earth_faa": GMTRemoteDataset(
         description="IGPP Earth free-air anomaly",
         units="mGal",
@@ -277,6 +295,24 @@ class GMTRemoteDataset(NamedTuple):
             "07m": Resolution("07m", registrations=["gridline"]),
         },
     ),
+    "earth_ndefl": GMTRemoteDataset(
+        description="IGPP Earth south-north deflections",
+        units="micro-radians",
+        extra_attributes={"horizontal_datum": "WGS84"},
+        resolutions={
+            "01d": Resolution("01d"),
+            "30m": Resolution("30m"),
+            "20m": Resolution("20m"),
+            "15m": Resolution("15m"),
+            "10m": Resolution("10m"),
+            "06m": Resolution("06m"),
+            "05m": Resolution("05m", tiled=True),
+            "04m": Resolution("04m", tiled=True),
+            "03m": Resolution("03m", tiled=True),
+            "02m": Resolution("02m", tiled=True),
+            "01m": Resolution("01m", registrations=["pixel"], tiled=True),
+        },
+    ),
     "earth_vgg": GMTRemoteDataset(
         description="IGPP Earth vertical gravity gradient",
         units="Eotvos",

pygmt/helpers/caching.py

@@ -15,6 +15,7 @@ def cache_data():
         "@earth_age_01d_g",
         "@earth_day_01d",
         "@earth_dist_01d",
+        "@earth_edefl_01d",
         "@earth_faa_01d_g",
         "@earth_gebco_01d_g",
         "@earth_gebcosi_01d_g",
@@ -26,6 +27,7 @@ def cache_data():
         "@earth_mdt_01d_g",
         "@earth_mdt_07m_g",
         "@earth_mss_01d_g",
+        "@earth_ndefl_01d",
         "@earth_night_01d",
         "@earth_relief_01d_g",
         "@earth_relief_01d_p",
@@ -50,11 +52,13 @@ def cache_data():
         "@N30E060.earth_age_01m_g.nc",
         "@N30E090.earth_age_01m_g.nc",
         "@N00W030.earth_dist_01m_g.nc",
+        "@N00W030.earth_edefl_01m_p.nc",
         "@N00W030.earth_faa_01m_p.nc",
         "@N00W030.earth_geoid_01m_g.nc",
         "@S30W060.earth_mag_02m_p.nc",
         "@S30W120.earth_mag4km_02m_p.nc",
         "@N30E090.earth_mss_01m_g.nc",
+        "@N30E090.earth_ndefl_01m_p.nc",
         "@N00W090.earth_relief_03m_p.nc",
         "@N00E135.earth_relief_30s_g.nc",
         "@N00W010.earth_relief_15s_p.nc",

pygmt/tests/test_datasets_earth_deflection.py

@@ -0,0 +1,106 @@
+"""
+Test basic functionality for loading IGPP Earth east-west and south-north deflection
+datasets.
+"""
+
+import numpy as np
+import numpy.testing as npt
+from pygmt.datasets import load_earth_deflection
+
+
+def test_earth_edefl_01d():
+    """
+    Test some properties of the Earth east-west deflection 01d data.
+    """
+    data = load_earth_deflection(resolution="01d")
+    assert data.name == "z"
+    assert data.attrs["long_name"] == "edefl (microradians)"
+    assert data.attrs["description"] == "IGPP Earth east-west deflections"
+    assert data.attrs["units"] == "micro-radians"
+    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(), -188.85, atol=0.04)
+    npt.assert_allclose(data.max(), 161.25, atol=0.04)
+
+
+def test_earth_edefl_01d_with_region():
+    """
+    Test loading low-resolution Earth east-west deflection with "region".
+    """
+    data = load_earth_deflection(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(), -36.125, atol=0.04)
+    npt.assert_allclose(data.max(), 45.3, atol=0.04)
+
+
+def test_earth_edefl_01m_default_registration():
+    """
+    Test that the grid returned by default for the 1 arc-minute resolution has a "pixel"
+    registration.
+    """
+    data = load_earth_deflection(resolution="01m", region=[-10, -9, 3, 5])
+    assert data.shape == (120, 60)
+    assert data.gmt.registration == 1
+    npt.assert_allclose(data.coords["lat"].data.min(), 3.008333333)
+    npt.assert_allclose(data.coords["lat"].data.max(), 4.991666666)
+    npt.assert_allclose(data.coords["lon"].data.min(), -9.99166666)
+    npt.assert_allclose(data.coords["lon"].data.max(), -9.00833333)
+    npt.assert_allclose(data.min(), -49.225, atol=0.04)
+    npt.assert_allclose(data.max(), 115.0, atol=0.04)
+
+
+def test_earth_ndefl_01d():
+    """
+    Test some properties of the Earth south-north deflection 01d data.
+    """
+    data = load_earth_deflection(resolution="01d", direction="ndefl")
+    assert data.name == "z"
+    assert data.attrs["long_name"] == "ndefl (microradians)"
+    assert data.attrs["description"] == "IGPP Earth south-north deflections"
+    assert data.attrs["units"] == "micro-radians"
+    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(), -188.85, atol=0.04)
+    npt.assert_allclose(data.max(), 161.25, atol=0.04)
+
+
+def test_earth_ndefl_01d_with_region():
+    """
+    Test loading low-resolution Earth south-north deflection with "region".
+    """
+    data = load_earth_deflection(
+        resolution="01d", region=[-10, 10, -5, 5], direction="ndefl"
+    )
+    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(), -36.125, atol=0.04)
+    npt.assert_allclose(data.max(), 45.3, atol=0.04)
+
+
+def test_earth_ndefl_01m_default_registration():
+    """
+    Test that the grid returned by default for the 1 arc-minute resolution has a "pixel"
+    registration.
+    """
+    data = load_earth_deflection(
+        resolution="01m", region=[-10, -9, 3, 5], direction="ndefl"
+    )
+    assert data.shape == (120, 60)
+    assert data.gmt.registration == 1
+    npt.assert_allclose(data.coords["lat"].data.min(), 3.008333333)
+    npt.assert_allclose(data.coords["lat"].data.max(), 4.991666666)
+    npt.assert_allclose(data.coords["lon"].data.min(), -9.99166666)
+    npt.assert_allclose(data.coords["lon"].data.max(), -9.00833333)
+    npt.assert_allclose(data.min(), -49.225, atol=0.04)
+    npt.assert_allclose(data.max(), 115.0, atol=0.04)