diff --git a/doc/api/index.rst b/doc/api/index.rst index 09d20fa1415..1f573dcb8eb 100644 --- a/doc/api/index.rst +++ b/doc/api/index.rst @@ -235,6 +235,7 @@ and store them in GMT's user data directory. datasets.load_black_marble datasets.load_blue_marble datasets.load_earth_age + datasets.load_earth_deflection datasets.load_earth_dist datasets.load_earth_free_air_anomaly datasets.load_earth_geoid diff --git a/pygmt/datasets/__init__.py b/pygmt/datasets/__init__.py index 9e163b33e3c..3d44a8fe676 100644 --- a/pygmt/datasets/__init__.py +++ b/pygmt/datasets/__init__.py @@ -6,6 +6,7 @@ from pygmt.datasets.earth_age import load_earth_age from pygmt.datasets.earth_day import load_blue_marble +from pygmt.datasets.earth_deflection import load_earth_deflection from pygmt.datasets.earth_dist import load_earth_dist from pygmt.datasets.earth_free_air_anomaly import load_earth_free_air_anomaly from pygmt.datasets.earth_geoid import load_earth_geoid diff --git a/pygmt/datasets/earth_deflection.py b/pygmt/datasets/earth_deflection.py new file mode 100644 index 00000000000..c0a9cbf406f --- /dev/null +++ b/pygmt/datasets/earth_deflection.py @@ -0,0 +1,119 @@ +""" +Function to download the IGPP Earth east-west and north-south 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, + component: Literal["east", "north"] = "east", +) -> xr.DataArray: + r""" + Load the IGPP Earth east-west and north-south deflection datasets in various + resolutions. + + .. list-table:: + :widths: 50 50 + :header-rows: 1 + + * - IGPP Earth east-west deflection + - IGPP Earth north-south 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``, which means + ``"gridline"`` for all resolutions except ``"01m"`` which is ``"pixel"`` only. + component + By default, the east-west deflection (``component="east"``) is returned, + set ``component="north"`` to return the north-south deflection. + + Returns + ------- + grid + The Earth east-west or north-south deflection grid. Coordinates are latitude + and longitude in degrees. Deflection values are in micro-radians, where + positive (negative) values indicate a deflection to the east or north (west + or south). + + 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 north-south deflection + >>> grid = load_earth_deflection(component="north") + >>> # 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 component == "north" else "earth_edefl" + grid = _load_remote_dataset( + name=prefix, + prefix=prefix, + resolution=resolution, + region=region, + registration=registration, + ) + return grid diff --git a/pygmt/datasets/load_remote_dataset.py b/pygmt/datasets/load_remote_dataset.py index da21cc7ae94..b19c055425c 100644 --- a/pygmt/datasets/load_remote_dataset.py +++ b/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 deflection", + 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", @@ -295,6 +313,24 @@ class GMTRemoteDataset(NamedTuple): "07m": Resolution("07m", registrations=["gridline"]), }, ), + "earth_ndefl": GMTRemoteDataset( + description="IGPP Earth north-south deflection", + 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", diff --git a/pygmt/helpers/caching.py b/pygmt/helpers/caching.py index 0175839e97d..19d3eae0559 100644 --- a/pygmt/helpers/caching.py +++ b/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_faaerror_01d_g", "@earth_gebco_01d_g", @@ -27,6 +28,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", @@ -51,12 +53,14 @@ 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_faaerror_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", diff --git a/pygmt/tests/test_datasets_earth_deflection.py b/pygmt/tests/test_datasets_earth_deflection.py new file mode 100644 index 00000000000..9118779a379 --- /dev/null +++ b/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 deflection" + 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(), -142.64, atol=0.04) + npt.assert_allclose(data.max(), 178.32, 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(), -28.92, atol=0.04) + npt.assert_allclose(data.max(), 24.72, 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(), -62.24, atol=0.04) + npt.assert_allclose(data.max(), 15.52, atol=0.04) + + +def test_earth_ndefl_01d(): + """ + Test some properties of the Earth north-south deflection 01d data. + """ + data = load_earth_deflection(resolution="01d", component="north") + assert data.name == "z" + assert data.attrs["long_name"] == "ndefl (microradians)" + assert data.attrs["description"] == "IGPP Earth north-south deflection" + 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(), -214.8, atol=0.04) + npt.assert_allclose(data.max(), 163.04, atol=0.04) + + +def test_earth_ndefl_01d_with_region(): + """ + Test loading low-resolution Earth north-south deflection with "region". + """ + data = load_earth_deflection( + resolution="01d", region=[-10, 10, -5, 5], component="north" + ) + 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(), -48.08, atol=0.04) + npt.assert_allclose(data.max(), 18.92, 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], component="north" + ) + 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(), -107.04, atol=0.04) + npt.assert_allclose(data.max(), 20.28, atol=0.04)