From 16919adbfeb4bfe36bb0c03be4270270d265878f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yvonne=20Fr=C3=B6hlich?= <94163266+yvonnefroehlich@users.noreply.github.com> Date: Thu, 26 Dec 2024 12:43:07 +0100 Subject: [PATCH] Remote datasets: Add "load_earth_mean_sea_surface" to load "CNES Earth Mean Sea Surface" dataset (#3717) --- doc/api/index.rst | 1 + pygmt/datasets/__init__.py | 1 + pygmt/datasets/earth_mean_sea_surface.py | 104 ++++++++++++++++++ pygmt/datasets/load_remote_dataset.py | 18 +++ pygmt/helpers/caching.py | 2 + .../test_datasets_earth_mean_sea_surface.py | 53 +++++++++ 6 files changed, 179 insertions(+) create mode 100644 pygmt/datasets/earth_mean_sea_surface.py create mode 100644 pygmt/tests/test_datasets_earth_mean_sea_surface.py diff --git a/doc/api/index.rst b/doc/api/index.rst index 6f6d87b8b7f..7f46afdf413 100644 --- a/doc/api/index.rst +++ b/doc/api/index.rst @@ -239,6 +239,7 @@ and store them in GMT's user data directory. datasets.load_earth_magnetic_anomaly datasets.load_earth_mask datasets.load_earth_mean_dynamic_topography + datasets.load_earth_mean_sea_surface datasets.load_earth_relief datasets.load_earth_vertical_gravity_gradient datasets.load_mars_relief diff --git a/pygmt/datasets/__init__.py b/pygmt/datasets/__init__.py index f60739a6e1b..9e163b33e3c 100644 --- a/pygmt/datasets/__init__.py +++ b/pygmt/datasets/__init__.py @@ -14,6 +14,7 @@ from pygmt.datasets.earth_mean_dynamic_topography import ( load_earth_mean_dynamic_topography, ) +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_vertical_gravity_gradient import ( diff --git a/pygmt/datasets/earth_mean_sea_surface.py b/pygmt/datasets/earth_mean_sea_surface.py new file mode 100644 index 00000000000..f4856d98626 --- /dev/null +++ b/pygmt/datasets/earth_mean_sea_surface.py @@ -0,0 +1,104 @@ +""" +Function to download the CNES Earth mean sea surface 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_sea_surface"] + + +def load_earth_mean_sea_surface( + resolution: Literal[ + "01d", "30m", "20m", "15m", "10m", "06m", "05m", "04m", "03m", "02m", "01m" + ] = "01d", + region: Sequence[float] | str | None = None, + registration: Literal["gridline", "pixel"] = "gridline", +) -> xr.DataArray: + r""" + Load the CNES Earth mean sea surface dataset in various resolutions. + + .. figure:: https://www.generic-mapping-tools.org/remote-datasets/_images/GMT_earth_mss.jpg + :width: 80 % + :align: center + + CNES Earth mean sea surface dataset. + + The grids are downloaded to a user data directory (usually + ``~/.gmt/server/earth/earth_mss/``) 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_mss**\_\ *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_mss.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_mss.cpt"``. + + Refer to :gmt-datasets:`earth-mss.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. + + Returns + ------- + grid + The CNES Earth mean sea surface 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_sea_surface + >>> # load the default grid (gridline-registered 1 arc-degree grid) + >>> grid = load_earth_mean_sea_surface() + >>> # load the 30 arc-minutes grid with "gridline" registration + >>> grid = load_earth_mean_sea_surface(resolution="30m", registration="gridline") + >>> # load high-resolution (5 arc-minutes) grid for a specific region + >>> grid = load_earth_mean_sea_surface( + ... resolution="05m", + ... region=[120, 160, 30, 60], + ... registration="gridline", + ... ) + """ + grid = _load_remote_dataset( + name="earth_mss", + prefix="earth_mss", + 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 9275d694707..0cd68796a7c 100644 --- a/pygmt/datasets/load_remote_dataset.py +++ b/pygmt/datasets/load_remote_dataset.py @@ -227,6 +227,24 @@ class GMTRemoteDataset(NamedTuple): "15s": Resolution("15s"), }, ), + "earth_mss": GMTRemoteDataset( + description="CNES Earth mean sea surface", + units="meters", + 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", tiled=True, registrations=["gridline"]), + }, + ), "earth_night": GMTRemoteDataset( description="NASA Night Images", units=None, diff --git a/pygmt/helpers/caching.py b/pygmt/helpers/caching.py index 054de7a4f85..d0bfd8fb46e 100644 --- a/pygmt/helpers/caching.py +++ b/pygmt/helpers/caching.py @@ -25,6 +25,7 @@ def cache_data(): "@earth_mask_01d_g", "@earth_mdt_01d_g", "@earth_mdt_07m_g", + "@earth_mss_01d_g", "@earth_night_01d", "@earth_relief_01d_g", "@earth_relief_01d_p", @@ -53,6 +54,7 @@ def cache_data(): "@N00W030.earth_geoid_01m_g.nc", "@S30W060.earth_mag_02m_p.nc", "@S30W120.earth_mag4km_02m_p.nc", + "@N30E090.earth_mss_01m_g.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_mean_sea_surface.py b/pygmt/tests/test_datasets_earth_mean_sea_surface.py new file mode 100644 index 00000000000..84b2b7123cc --- /dev/null +++ b/pygmt/tests/test_datasets_earth_mean_sea_surface.py @@ -0,0 +1,53 @@ +""" +Test basic functionality for loading Earth mean sea surface datasets. +""" + +import numpy as np +import numpy.testing as npt +from pygmt.datasets import load_earth_mean_sea_surface + + +def test_earth_mss_01d(): + """ + Test some properties of the Earth mean sea surface 01d data. + """ + data = load_earth_mean_sea_surface(resolution="01d") + assert data.name == "z" + assert data.attrs["description"] == "CNES Earth mean sea surface" + 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(), -104.71, atol=0.01) + npt.assert_allclose(data.max(), 82.38, atol=0.01) + + +def test_earth_mss_01d_with_region(): + """ + Test loading low-resolution Earth mean sea surface with "region". + """ + data = load_earth_mean_sea_surface(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(), 6.53, atol=0.01) + npt.assert_allclose(data.max(), 29.31, atol=0.01) + + +def test_earth_mss_01m_default_registration(): + """ + Test that the grid returned by default for the 1 arc-minute resolution has a + "gridline" registration. + """ + data = load_earth_mean_sea_surface(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(), 21.27, atol=0.01) + npt.assert_allclose(data.max(), 31.11, atol=0.01)