Let load_earth_relief() support 'region' argument for all resolutions (#873)

* Let load_earth_relief() support 'region' argument for all resolutions

* Improve docstrings of load_earth_relief()

Author: seisman

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
     --------
@@ -83,26 +87,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"

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():