diff --git a/pygmt/datasets/earth_relief.py b/pygmt/datasets/earth_relief.py index 9a1bd911e3e..980f405337a 100644 --- a/pygmt/datasets/earth_relief.py +++ b/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) diff --git a/pygmt/helpers/testing.py b/pygmt/helpers/testing.py index fbb17339258..13fc1766abb 100644 --- a/pygmt/helpers/testing.py +++ b/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", diff --git a/pygmt/tests/test_datasets_earth_relief.py b/pygmt/tests/test_datasets_earth_relief.py index e6f1ed40406..3529e555df4 100644 --- a/pygmt/tests/test_datasets_earth_relief.py +++ b/pygmt/tests/test_datasets_earth_relief.py @@ -8,7 +8,8 @@ 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. """ @@ -16,15 +17,17 @@ def test_earth_relief_fails(): 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)) @@ -32,12 +35,29 @@ def test_earth_relief_01d(): 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", + )