From d07ff6bc079457a2527f76dcc1c88ca240e58f19 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Fri, 12 Feb 2021 13:03:09 -0500 Subject: [PATCH 1/2] Let load_earth_relief() support 'region' argument for all resolutions --- pygmt/datasets/earth_relief.py | 30 ++++++++++++++---------------- pygmt/tests/test_datasets.py | 10 ++++++++-- 2 files changed, 22 insertions(+), 18 deletions(-) diff --git a/pygmt/datasets/earth_relief.py b/pygmt/datasets/earth_relief.py index d3d57114173..2a3bf5f52a5 100644 --- a/pygmt/datasets/earth_relief.py +++ b/pygmt/datasets/earth_relief.py @@ -83,26 +83,24 @@ def load_earth_relief(resolution="01d", region=None, registration=None): "gridline-registered grid is available." ) + if resolution not in non_tiled_resolutions + tiled_resolutions: + raise GMTInvalidInput(f"Invalid Earth relief resolution '{resolution}'.") + # different ways to load tiled and non-tiled earth relief data - if resolution in non_tiled_resolutions: - if region is not None: - raise NotImplementedError( - f"'region' is not supported for Earth relief resolution '{resolution}'" - ) - fname = which(f"@earth_relief_{resolution}{reg}", download="a") - with xr.open_dataarray(fname) as dataarray: - grid = dataarray.load() - _ = grid.gmt # load GMTDataArray accessor information - elif resolution in tiled_resolutions: - # Titled grid can't be sliced. - # See https://github.com/GenericMappingTools/pygmt/issues/524 - if region is None: + # Known issue: tiled grids don't support slice operation + # See https://github.com/GenericMappingTools/pygmt/issues/524 + if region is None: + if resolution in non_tiled_resolutions: + fname = which(f"@earth_relief_{resolution}{reg}", download="a") + with xr.open_dataarray(fname) as dataarray: + grid = dataarray.load() + _ = grid.gmt # load GMTDataArray accessor information + else: raise GMTInvalidInput( - f"'region' is required for Earth relief resolution '{resolution}'" + f"'region' is required for Earth relief resolution '{resolution}'." ) - grid = grdcut(f"@earth_relief_{resolution}{reg}", region=region) else: - raise GMTInvalidInput(f'Invalid Earth relief resolution "{resolution}"') + grid = grdcut(f"@earth_relief_{resolution}{reg}", region=region) # Add some metadata to the grid grid.name = "elevation" diff --git a/pygmt/tests/test_datasets.py b/pygmt/tests/test_datasets.py index 72bdcca3d20..6bdde8da40b 100644 --- a/pygmt/tests/test_datasets.py +++ b/pygmt/tests/test_datasets.py @@ -93,8 +93,14 @@ def test_earth_relief_01d_with_region(): """ Test loading low-resolution earth relief with 'region'. """ - with pytest.raises(NotImplementedError): - load_earth_relief("01d", region=[0, 180, 0, 90]) + data = load_earth_relief( + 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(), -5145) + npt.assert_allclose(data.max(), 805.5) def test_earth_relief_30m(): From 43e14ded930c20e2c77a110c8131d4768a65d809 Mon Sep 17 00:00:00 2001 From: Dongdong Tian Date: Fri, 12 Feb 2021 13:30:51 -0500 Subject: [PATCH 2/2] Improve docstrings of load_earth_relief() --- pygmt/datasets/earth_relief.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/pygmt/datasets/earth_relief.py b/pygmt/datasets/earth_relief.py index 2a3bf5f52a5..c906d064f70 100644 --- a/pygmt/datasets/earth_relief.py +++ b/pygmt/datasets/earth_relief.py @@ -1,6 +1,6 @@ """ Function to download the Earth relief datasets from the GMT data server, and -load as DataArray. +load as :class:`xarray.DataArray`. The grids are available in various resolutions. """ @@ -12,7 +12,7 @@ @kwargs_to_strings(region="sequence") def load_earth_relief(resolution="01d", region=None, registration=None): - """ + r""" Load Earth relief grids (topography and bathymetry) in various resolutions. The grids are downloaded to a user data directory @@ -21,8 +21,11 @@ def load_earth_relief(resolution="01d", region=None, registration=None): 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_rru[_reg]'`` to any grid plotting/processing function. - Refer to :gmt-docs:`datasets/remote-data.html` for more details. + **@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). Refer to :gmt-docs:`datasets/remote-data.html` for more + details. Parameters ---------- @@ -35,7 +38,7 @@ def load_earth_relief(resolution="01d", region=None, registration=None): region : str or list The subregion of the grid to load. Required for Earth relief grids with - resolutions <= 05m. + resolutions higher than 5 arc-minute (i.e., ``05m``). registration : str Grid registration type. Either ``pixel`` for pixel registration or @@ -45,14 +48,15 @@ def load_earth_relief(resolution="01d", region=None, registration=None): Returns ------- - grid : xarray.DataArray + grid : :class:`xarray.DataArray` The Earth relief grid. Coordinates are latitude and longitude in degrees. Relief is in meters. Notes ----- - The DataArray doesn's support slice operation, for Earth relief data with - resolutions higher than "05m", which are stored as smaller tiles. + The :class:`xarray.DataArray` grid doesn't support slice operation, for + Earth relief data with resolutions higher than "05m", which are stored as + smaller tiles. Examples --------