Add load_earth_free_air_anomaly function for Earth free-air anomaly d…

…ataset (#2238)

* Add earth_free_air_anomaly.py
* Add load_earth_free_air_anomaly to API index
* Add tests and cache files for load_earth_free_air_anomaly

Co-authored-by: Dongdong Tian <[email protected]>
Co-authored-by: Michael Grund <[email protected]>

doc/api/index.rst

@@ -222,6 +222,7 @@ and store them in GMT's user data directory.
 
     datasets.list_sample_data
     datasets.load_earth_age
+    datasets.load_earth_free_air_anomaly
     datasets.load_earth_geoid
     datasets.load_earth_magnetic_anomaly
     datasets.load_earth_relief

pygmt/datasets/__init__.py

@@ -3,6 +3,7 @@
 # Load sample data included with GMT (downloaded from the GMT cache server).
 
 from pygmt.datasets.earth_age import load_earth_age
+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
 from pygmt.datasets.earth_relief import load_earth_relief

pygmt/datasets/earth_free_air_anomaly.py

@@ -0,0 +1,68 @@
+"""
+Function to download the IGPP Global Earth Free-Air Anomaly datasets from the
+GMT data server, and load as :class:`xarray.DataArray`.
+
+The grids are available in various resolutions.
+"""
+from pygmt.datasets.load_remote_dataset import _load_remote_dataset
+from pygmt.helpers import kwargs_to_strings
+
+
+@kwargs_to_strings(region="sequence")
+def load_earth_free_air_anomaly(resolution="01d", region=None, registration=None):
+    r"""
+    Load an Earth Free-Air Anomaly grid in various resolutions.
+
+    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.
+
+    These grids can also be accessed by passing in the file name
+    **@earth_faa**\_\ *res*\[_\ *reg*] to any grid plotting/processing
+    function. *res* is the grid resolution (see below), and *reg* is grid
+    registration type (**p** for pixel registration or **g** for gridline
+    registration).
+
+    Refer to :gmt-datasets:`earth-faa.html` for more details.
+
+    Parameters
+    ----------
+    resolution : str
+        The grid resolution. The suffix ``d`` and ``m`` stand for
+        arc-degree and arc-minute. It can be ``"01d"``, ``"30m"``,
+        ``"20m"``, ``"15m"``, ``"10m"``, ``"06m"``, ``"05m"``, ``"04m"``,
+        ``"03m"``, ``"02m"``, or ``"01m"``.
+
+    region : str or list
+        The subregion of the grid to load, in the forms of a list
+        [*xmin*, *xmax*, *ymin*, *ymax*] or a string *xmin/xmax/ymin/ymax*.
+        Required for grids with resolutions higher than 5
+        arc-minutes (i.e., ``"05m"``).
+
+    registration : str
+        Grid registration type. Either ``"pixel"`` for pixel registration or
+        ``"gridline"`` for gridline registration. Default is ``None``, where
+        a pixel-registered grid is returned unless only the
+        gridline-registered grid is available.
+
+    Returns
+    -------
+    grid : :class:`xarray.DataArray`
+        The Earth free-air anomaly grid. Coordinates are latitude and
+        longitude in degrees. Units are in mGal.
+
+    Note
+    ----
+    The :class:`xarray.DataArray` grid doesn't support slice operation, for
+    Earth free-air anomaly with resolutions of 5 arc-minutes or higher,
+    which are stored as smaller tiles.
+    """
+    grid = _load_remote_dataset(
+        dataset_name="earth_free_air_anomaly",
+        dataset_prefix="earth_faa_",
+        resolution=resolution,
+        region=region,
+        registration=registration,
+    )
+    return grid

pygmt/datasets/load_remote_dataset.py

@@ -147,6 +147,26 @@ class GMTRemoteDataset(NamedTuple):
             "02m": Resolution(["pixel"], True),
         },
     ),
+    "earth_free_air_anomaly": GMTRemoteDataset(
+        title="free air anomaly",
+        name="free_air_anomaly",
+        long_name="IGPP Global Earth Free-Air Anomaly",
+        units="mGal",
+        extra_attributes={"horizontal_datum": "WGS84"},
+        resolutions={
+            "01d": Resolution(["pixel", "gridline"], False),
+            "30m": Resolution(["pixel", "gridline"], False),
+            "20m": Resolution(["pixel", "gridline"], False),
+            "15m": Resolution(["pixel", "gridline"], False),
+            "10m": Resolution(["pixel", "gridline"], False),
+            "06m": Resolution(["pixel", "gridline"], False),
+            "05m": Resolution(["pixel", "gridline"], True),
+            "04m": Resolution(["pixel", "gridline"], True),
+            "03m": Resolution(["pixel", "gridline"], True),
+            "02m": Resolution(["pixel", "gridline"], True),
+            "01m": Resolution(["pixel"], True),
+        },
+    ),
 }
 
 

pygmt/helpers/testing.py

@@ -181,13 +181,16 @@ def download_test_data():
         "@earth_age_01d_g",
         "@S90W180.earth_age_05m_g.nc",  # Specific grid for 05m test
         # Earth geoid grids
-        "@earth_geoid_01d_g.grd",
+        "@earth_geoid_01d_g",
         "@S90W180.earth_geoid_05m_g.nc",  # Specific grid for 05m test
         # Earth magnetic anomaly grids
         "@earth_mag_01d_g",
         "@S90W180.earth_mag_05m_g.nc",  # Specific grid for 05m test
         "@earth_mag4km_01d_g",
         "@S90W180.earth_mag4km_05m_g.nc",  # Specific grid for 05m test
+        # Earth free-air anomaly grids
+        "@earth_faa_01d_g",
+        "@S90W180.earth_faa_05m_g.nc",  # Specific grid for 05m test
         # Other cache files
         "@capitals.gmt",
         "@earth_relief_20m_holes.grd",

pygmt/tests/test_datasets_earth_free_air_anomaly.py

@@ -0,0 +1,82 @@
+"""
+Test basic functionality for loading Earth free air anomaly datasets.
+"""
+import numpy as np
+import numpy.testing as npt
+import pytest
+from pygmt.datasets import load_earth_free_air_anomaly
+from pygmt.exceptions import GMTInvalidInput
+
+
+def test_earth_faa_fails():
+    """
+    Make sure earth_free_air_anomaly fails for invalid resolutions.
+    """
+    resolutions = "1m 1d bla 60d 001m 03".split()
+    resolutions.append(60)
+    for resolution in resolutions:
+        with pytest.raises(GMTInvalidInput):
+            load_earth_free_air_anomaly(resolution=resolution)
+
+
+def test_earth_faa_incorrect_registration():
+    """
+    Test loading earth_free_air_anomaly with incorrect registration type.
+    """
+    with pytest.raises(GMTInvalidInput):
+        load_earth_free_air_anomaly(registration="improper_type")
+
+
+def test_earth_faa_01d():
+    """
+    Test some properties of the free air anomaly 01d data.
+    """
+    data = load_earth_free_air_anomaly(resolution="01d", registration="gridline")
+    assert data.name == "free_air_anomaly"
+    assert data.attrs["long_name"] == "IGPP Global Earth Free-Air Anomaly"
+    assert data.attrs["units"] == "mGal"
+    assert data.attrs["horizontal_datum"] == "WGS84"
+    assert data.shape == (181, 361)
+    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(), -275.75)
+    npt.assert_allclose(data.max(), 308.375)
+
+
+def test_earth_faa_01d_with_region():
+    """
+    Test loading low-resolution earth free air anomaly with 'region'.
+    """
+    data = load_earth_free_air_anomaly(
+        resolution="01d", region=[-10, 10, -5, 5], registration="gridline"
+    )
+    assert data.shape == (11, 21)
+    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(), -58.75)
+    npt.assert_allclose(data.max(), 69.524994)
+
+
+def test_earth_faa_05m_with_region():
+    """
+    Test loading a subregion of high-resolution earth free air anomaly data.
+    """
+    data = load_earth_free_air_anomaly(
+        resolution="05m", region=[-115, -112, 4, 6], registration="gridline"
+    )
+    assert data.shape == (25, 37)
+    assert data.lat.min() == 4
+    assert data.lat.max() == 6
+    assert data.lon.min() == -115
+    assert data.lon.max() == -112
+    npt.assert_allclose(data.min(), -20.5)
+    npt.assert_allclose(data.max(), -3.9500122)
+
+
+def test_earth_faa_05m_without_region():
+    """
+    Test loading high-resolution earth free air anomaly without passing
+    'region'.
+    """
+    with pytest.raises(GMTInvalidInput):
+        load_earth_free_air_anomaly("05m")