From 230ff31b4a6dd9770a226e1d23d7e3cbf3a3a7ec Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yvonne=20Fr=C3=B6hlich?= <94163266+yvonnefroehlich@users.noreply.github.com> Date: Sat, 28 Dec 2024 15:15:40 +0100 Subject: [PATCH] Remote datasets: Add "uncertainty" parameter to "load_earth_free_air_anomaly" to load the "free-air anomaly uncertainty" dataset (#3727) --- pygmt/datasets/earth_free_air_anomaly.py | 93 ++++++++++--------- pygmt/datasets/load_remote_dataset.py | 18 ++++ pygmt/helpers/caching.py | 2 + .../test_datasets_earth_free_air_anomaly.py | 51 ++++++++++ 4 files changed, 122 insertions(+), 42 deletions(-) diff --git a/pygmt/datasets/earth_free_air_anomaly.py b/pygmt/datasets/earth_free_air_anomaly.py index da48977d688..df1154f2aaa 100644 --- a/pygmt/datasets/earth_free_air_anomaly.py +++ b/pygmt/datasets/earth_free_air_anomaly.py @@ -1,6 +1,6 @@ """ -Function to download the IGPP Earth free-air anomaly dataset from the GMT data server, -and load as :class:`xarray.DataArray`. +Function to download the IGPP Earth free-air anomaly and uncertainty datasets from +the GMT data server, and load as :class:`xarray.DataArray`. The grids are available in various resolutions. """ @@ -20,36 +20,43 @@ def load_earth_free_air_anomaly( ] = "01d", region: Sequence[float] | str | None = None, registration: Literal["gridline", "pixel", None] = None, + uncertainty: bool = False, ) -> xr.DataArray: r""" - Load the IGPP Earth free-air anomaly dataset in various resolutions. + Load the IGPP Earth free-air anomaly and uncertainty datasets in various + resolutions. - .. figure:: https://www.generic-mapping-tools.org/remote-datasets/_images/GMT_earth_faa.jpg - :width: 80 % - :align: center + .. list-table:: + :widths: 50 50 + :header-rows: 1 - IGPP Earth free-air anomaly dataset. + * - IGPP Earth free-Air anomaly + - IGPP Earth free-Air anomaly uncertainty + * - .. figure:: https://www.generic-mapping-tools.org/remote-datasets/_images/GMT_earth_faa.jpg + - .. figure:: https://www.generic-mapping-tools.org/remote-datasets/_images/GMT_earth_faaerror.jpg - The grids are downloaded to a user data directory - (usually ``~/.gmt/server/earth/earth_faa/``) 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. + The grids are downloaded to a user data directory (usually + ``~/.gmt/server/earth/earth_faa/`` or ``~/.gmt/server/earth/earth_faaerror/``) the + first time you invoke this function. Afterwards, it will load the grid from 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_faa**\_\ *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_faa.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_faa.cpt"``. - - Refer to :gmt-datasets:`earth-faa.html` for more details about available - datasets, including version information and references. + **@earth_faa_type**\_\ *res*\[_\ *reg*] to any grid processing function or + plotting method. *earth_faa_type* is the GMT name for the dataset. The available + options are **earth_faa** and **earth_faaerror**. *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 tables (CPTs) for these datasets are *@earth_faa.cpt* and + *@earth_faaerror.cpt*. The dataset-specific CPT is 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_faa.cpt"`` or + ``cmap="@earth_faaerror.cpt"``. + + Refer to :gmt-datasets:`earth-faa.html` and :gmt-datasets:`earth-faaerror.html` for + more details about available datasets, including version information and references. Parameters ---------- @@ -62,26 +69,27 @@ def load_earth_free_air_anomaly( 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. + ``"gridline"`` for gridline registration. Default is ``None`` which means + ``"gridline"`` for all resolutions except ``"01m"`` which is ``"pixel"`` only. + uncertainty + By default, the Earth free-air anomaly values are returned. Set to ``True`` to + return the related uncertainties instead. Returns ------- grid - The Earth free-air anomaly grid. Coordinates are latitude and - longitude in degrees. Units are in mGal. + The Earth free-air anomaly (uncertainty) grid. Coordinates are latitude and + longitude in degrees. Values and uncertainties are in mGal. 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. + :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 -------- @@ -89,18 +97,19 @@ def load_earth_free_air_anomaly( >>> from pygmt.datasets import load_earth_free_air_anomaly >>> # load the default grid (gridline-registered 1 arc-degree grid) >>> grid = load_earth_free_air_anomaly() + >>> # load the uncertainties related to the default grid + >>> grid = load_earth_free_air_anomaly(uncertainty=True) >>> # load the 30 arc-minutes grid with "gridline" registration >>> grid = load_earth_free_air_anomaly(resolution="30m", registration="gridline") >>> # load high-resolution (5 arc-minutes) grid for a specific region >>> grid = load_earth_free_air_anomaly( - ... resolution="05m", - ... region=[120, 160, 30, 60], - ... registration="gridline", + ... resolution="05m", region=[120, 160, 30, 60], registration="gridline" ... ) """ + prefix = "earth_faaerror" if uncertainty is True else "earth_faa" grid = _load_remote_dataset( - name="earth_faa", - prefix="earth_faa", + name=prefix, + prefix=prefix, resolution=resolution, region=region, registration=registration, diff --git a/pygmt/datasets/load_remote_dataset.py b/pygmt/datasets/load_remote_dataset.py index 3949a2e6f95..da21cc7ae94 100644 --- a/pygmt/datasets/load_remote_dataset.py +++ b/pygmt/datasets/load_remote_dataset.py @@ -128,6 +128,24 @@ class GMTRemoteDataset(NamedTuple): "01m": Resolution("01m", registrations=["pixel"], tiled=True), }, ), + "earth_faaerror": GMTRemoteDataset( + description="IGPP Earth free-air anomaly errors", + units="mGal", + 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_gebco": GMTRemoteDataset( description="GEBCO Earth relief", units="meters", diff --git a/pygmt/helpers/caching.py b/pygmt/helpers/caching.py index d0bfd8fb46e..0175839e97d 100644 --- a/pygmt/helpers/caching.py +++ b/pygmt/helpers/caching.py @@ -16,6 +16,7 @@ def cache_data(): "@earth_day_01d", "@earth_dist_01d", "@earth_faa_01d_g", + "@earth_faaerror_01d_g", "@earth_gebco_01d_g", "@earth_gebcosi_01d_g", "@earth_gebcosi_15m_p", @@ -51,6 +52,7 @@ def cache_data(): "@N30E090.earth_age_01m_g.nc", "@N00W030.earth_dist_01m_g.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", diff --git a/pygmt/tests/test_datasets_earth_free_air_anomaly.py b/pygmt/tests/test_datasets_earth_free_air_anomaly.py index 517a1bb9b89..7ce5fc2662c 100644 --- a/pygmt/tests/test_datasets_earth_free_air_anomaly.py +++ b/pygmt/tests/test_datasets_earth_free_air_anomaly.py @@ -52,3 +52,54 @@ def test_earth_faa_01m_default_registration(): npt.assert_allclose(data.coords["lon"].data.max(), -9.00833333) npt.assert_allclose(data.min(), -49.225, atol=0.025) npt.assert_allclose(data.max(), 115.0, atol=0.025) + + +def test_earth_faaerror_01d(): + """ + Test some properties of the free air anomaly error 01d data. + """ + data = load_earth_free_air_anomaly(resolution="01d", uncertainty=True) + assert data.name == "z" + assert data.attrs["long_name"] == "faaerror (mGal)" + assert data.attrs["description"] == "IGPP Earth free-air anomaly errors" + assert data.attrs["units"] == "mGal" + 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(), 0.0, atol=0.04) + npt.assert_allclose(data.max(), 49.16, atol=0.04) + + +def test_earth_faaerror_01d_with_region(): + """ + Test loading low-resolution earth free air anomaly error with 'region'. + """ + data = load_earth_free_air_anomaly( + resolution="01d", region=[-10, 10, -5, 5], uncertainty=True + ) + 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.72, atol=0.04) + npt.assert_allclose(data.max(), 21.04, atol=0.04) + + +def test_earth_faaerror_01m_default_registration(): + """ + Test that the grid returned by default for the 1 arc-minute resolution has a "pixel" + registration. + """ + data = load_earth_free_air_anomaly( + resolution="01m", region=[-10, -9, 3, 5], uncertainty=True + ) + 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(), 0.40, atol=0.04) + npt.assert_allclose(data.max(), 13.36, atol=0.04)