Remote datasets: Add "uncertainty" parameter to "load_earth_free_air_anomaly" to load the "free-air anomaly uncertainty" dataset (#3727)

Author: yvonnefroehlich

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,45 +69,47 @@ 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
     --------
 
     >>> 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,

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",

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",

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)