Add GEBCO data source to load_earth_relief (#1818)

Co-authored-by: Dongdong Tian <seisman.info@gmail.com>
Co-authored-by: Wei Ji <23487320+weiji14@users.noreply.github.com>

pygmt/datasets/earth_relief.py

@@ -4,29 +4,38 @@
 
 The grids are available in various resolutions.
 """
-from pygmt.exceptions import GMTInvalidInput
+from packaging.version import Version
+from pygmt.clib import Session
+from pygmt.exceptions import GMTInvalidInput, GMTVersionError
 from pygmt.helpers import kwargs_to_strings
 from pygmt.io import load_dataarray
 from pygmt.src import grdcut, which
 
 
 @kwargs_to_strings(region="sequence")
-def load_earth_relief(resolution="01d", region=None, registration=None, use_srtm=False):
+def load_earth_relief(
+    resolution="01d",
+    region=None,
+    registration=None,
+    use_srtm=False,
+    data_source="igpp",
+):
     r"""
     Load Earth relief grids (topography and bathymetry) in various resolutions.
 
     The grids are downloaded to a user data directory
-    (usually ``~/.gmt/server/earth/earth_relief/``) 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.
+    (usually a subdirectory under ~/.gmt/server/earth/) 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_relief**\_\ *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).
+    This module downloads the grids that can also be accessed by
+    passing in the file name **@**\ *earth_relief_type*\_\ *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-relief.html` for more details.
+    Refer to :gmt-datasets:`earth-relief.html` for more details about available
+    datasets, including version information and references.
 
     Parameters
     ----------
@@ -55,9 +64,21 @@ def load_earth_relief(resolution="01d", region=None, registration=None, use_srtm
     use_srtm : bool
         By default, the land-only SRTM tiles from NASA are used to generate the
         ``"03s"`` and ``"01s"`` grids, and the missing ocean values are filled
-        by up-sampling the SRTM15+ tiles which have a resolution of 15
+        by up-sampling the SRTM15 tiles which have a resolution of 15
         arc-second (i.e., ``"15s"``). If True, will only load the original
-        land-only SRTM tiles.
+        land-only SRTM tiles. Only works when ``data_source="igpp"``.
+
+    data_source : str
+        Select the source for the Earth relief data.
+
+        Available options:
+
+        - **igpp** : IGPP Global Earth Relief [Default option]. See
+          :gmt-datasets:`earth-relief.html`.
+
+        - **gebco** : GEBCO Global Earth Relief with only observed relief and
+          inferred relief via altimetric gravity. See
+          :gmt-datasets:`earth-gebco.html`.
 
     Returns
     -------
@@ -90,7 +111,7 @@ def load_earth_relief(resolution="01d", region=None, registration=None, use_srtm
     ...     use_srtm=True,
     ... )
     """
-
+    # pylint: disable=too-many-branches
     # earth relief data stored as single grids for low resolutions
     non_tiled_resolutions = ["01d", "30m", "20m", "15m", "10m", "06m"]
     # earth relief data stored as tiles for high resolutions
@@ -120,11 +141,30 @@ def load_earth_relief(resolution="01d", region=None, registration=None, use_srtm
             f"{registration}-registered Earth relief data for "
             f"resolution '{resolution}' is not supported."
         )
-
+    earth_relief_sources = {"igpp": "earth_relief_", "gebco": "earth_gebco_"}
+    if data_source not in earth_relief_sources:
+        raise GMTInvalidInput(
+            f"Invalid earth relief 'data_source' {data_source}, "
+            "valid values are 'igpp' and 'gebco'."
+        )
+    if data_source != "igpp":
+        with Session() as lib:
+            if Version(lib.info["version"]) < Version("6.4.0"):
+                raise GMTVersionError(
+                    f"The {data_source} option is not available for GMT"
+                    " versions before 6.4.0."
+                )
     # Choose earth relief data prefix
-    earth_relief_prefix = "earth_relief_"
     if use_srtm and resolution in land_only_srtm_resolutions:
-        earth_relief_prefix = "srtm_relief_"
+        if data_source == "igpp":
+            earth_relief_prefix = "srtm_relief_"
+        else:
+            raise GMTInvalidInput(
+                f"The {data_source} option is not available if 'use_srtm=True'."
+                " Set data_source to 'igpp'."
+            )
+    else:
+        earth_relief_prefix = earth_relief_sources.get(data_source)
 
     # different ways to load tiled and non-tiled earth relief data
     # Known issue: tiled grids don't support slice operation
@@ -134,7 +174,7 @@ def load_earth_relief(resolution="01d", region=None, registration=None, use_srtm
             raise GMTInvalidInput(
                 f"'region' is required for Earth relief resolution '{resolution}'."
             )
-        fname = which(f"@earth_relief_{resolution}{reg}", download="a")
+        fname = which(f"@{earth_relief_prefix}{resolution}{reg}", download="a")
         grid = load_dataarray(fname, engine="netcdf4")
     else:
         grid = grdcut(f"@{earth_relief_prefix}{resolution}{reg}", region=region)

pygmt/helpers/testing.py

@@ -153,6 +153,7 @@ def download_test_data():
     # List of datasets to download
     datasets = [
         # Earth relief grids
+        "@earth_gebco_01d_g",
         "@earth_relief_01d_p",
         "@earth_relief_01d_g",
         "@earth_relief_30m_p",

pygmt/tests/test_datasets_earth_relief.py

@@ -8,36 +8,56 @@
 from pygmt.exceptions import GMTInvalidInput
 
 
-def test_earth_relief_fails():
+@pytest.mark.parametrize("data_source", ["igpp", "gebco"])
+def test_earth_relief_fails(data_source):
     """
     Make sure earth relief 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_relief(resolution=resolution)
+            load_earth_relief(resolution=resolution, data_source=data_source)
 
 
 # Only test 01d and 30m to avoid downloading large datasets in CI
-def test_earth_relief_01d():
+def test_earth_relief_01d_igpp():
     """
-    Test some properties of the earth relief 01d data.
+    Test some properties of the earth relief 01d data with IGPP data.
     """
-    data = load_earth_relief(resolution="01d", registration="gridline")
+    data = load_earth_relief(
+        resolution="01d", registration="gridline", data_source="igpp"
+    )
     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(), -8600.5)
     npt.assert_allclose(data.max(), 5559.0)
 
 
-def test_earth_relief_01d_with_region():
+def test_earth_relief_01d_gebco():
+    """
+    Test some properties of the earth relief 01d data with GEBCO data.
+    """
+    data = load_earth_relief(
+        resolution="01d", registration="gridline", data_source="gebco"
+    )
+    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(), -8598)
+    npt.assert_allclose(data.max(), 5559.0)
+
+
+def test_earth_relief_01d_with_region_srtm():
     """
-    Test loading low-resolution earth relief with 'region'.
+    Test loading low-resolution earth relief with 'region' with IGPP data.
     """
     data = load_earth_relief(
-        resolution="01d", region=[-10, 10, -5, 5], registration="gridline"
+        resolution="01d",
+        region=[-10, 10, -5, 5],
+        registration="gridline",
+        data_source="igpp",
     )
     assert data.shape == (11, 21)
     npt.assert_allclose(data.lat, np.arange(-5, 6, 1))
@@ -46,6 +66,23 @@ def test_earth_relief_01d_with_region():
     npt.assert_allclose(data.max(), 805.5)
 
 
+def test_earth_relief_01d_with_region_gebco():
+    """
+    Test loading low-resolution earth relief with 'region' with GEBCO data.
+    """
+    data = load_earth_relief(
+        resolution="01d",
+        region=[-10, 10, -5, 5],
+        registration="gridline",
+        data_source="gebco",
+    )
+    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(), -5146)
+    npt.assert_allclose(data.max(), 806)
+
+
 def test_earth_relief_30m():
     """
     Test some properties of the earth relief 30m data.
@@ -122,3 +159,28 @@ def test_earth_relief_invalid_resolution_registration_combination():
     ]:
         with pytest.raises(GMTInvalidInput):
             load_earth_relief(resolution=resolution, registration=registration)
+
+
+def test_earth_relief_invalid_data_source():
+    """
+    Test loading earth relief with invalid data_source argument.
+    """
+    with pytest.raises(GMTInvalidInput):
+        load_earth_relief(
+            resolution="01d", registration="gridline", data_source="invalid_source"
+        )
+
+
+def test_earth_relief_invalid_data_source_with_use_srtm():
+    """
+    Test loading earth relief with use_srtm=True and an incompatible
+    data_source argument.
+    """
+    with pytest.raises(GMTInvalidInput):
+        load_earth_relief(
+            resolution="03s",
+            region=[135, 136, 35, 36],
+            registration="gridline",
+            use_srtm=True,
+            data_source="gebco",
+        )