Repository: GenericMappingTools/pygmt

Issue:
Allow load_earth_relief() to load the original land-only SRTM tiles
**Description of the desired feature**

Refer to https://docs.generic-mapping-tools.org/latest/datasets/remote-data.html#technical-information for details.

> The 3 and 1 arc second data are the SRTM 1x1 degree tiles from NASA. Note: The 3 and 1 arc second grids only extend to latitudes ±60˚ and are only available over land. When these grids are accessed as @earth_relief_01s or @earth_relief_03s we will automatically up-sample the relevant @earth_relief_15s tiles to fill in the missing ocean values. If you just want the original land-only SRTM tiles you may use the special names @srtm_relief_03s or @srtm_relief_01s instead. 

Currently, the `load_earth_relief()` function can load `01s` and `03s` Earth relief grid, but can't load the original land-only SRTM grid. 

It seems easy to add the feature. We just need to add a new parameter `use_srtm` (default to `False`) to the `load_earth_relief` function.

If `resolution="01s"` or `resolution="03s"`, we check the value of `use_srtm` to determine the remote file prefix (either `earth_relief_` or `srtm_relief_`).

**Are you willing to help implement and maintain this feature?** Yes, but volunteers are welcomed.

<!-- Every feature we add is code that we will have to maintain and keep updated. This takes a lot of effort. If you are willing to be involved in the project and help maintain your feature, it will make it easier for us to accept it. -->


Discussion:
seisman: @core-man I just converted your PR to draft to save CI resources.
core-man: > We also need to add a test to make sure that the change work as expected.
> 
> Please see `pygmt/tests/test_datasets.py` for similar tests.

`test_earth_relief_03s_landonly_srtm` was added in 849df3b.
seisman: @core-man Sorry, I didn't realize that you already added a test to `tests/test_datasets.py`.

We just split the tests in `tests/test_dataset.py` into two files `tests/test_datasets_tutorial.py` and `tests/test_datasets_earth_relief.py` (https://github.com/GenericMappingTools/pygmt/pull/978). So you have to deal with the conflict and move your new test into the file `tests/test_datasets_earth_relief.py`.

Sorry for the trouble. Please let me know if you need any help resolving the conflicts.
core-man: > @core-man Sorry, I didn't realize that you already added a test to `tests/test_datasets.py`.
> 
> We just split the tests in `tests/test_dataset.py` into two files `tests/test_datasets_tutorial.py` and `tests/test_datasets_earth_relief.py` (#978). So you have to deal with the conflict and move your new test into the file `tests/test_datasets_earth_relief.py`.
> 
> Sorry for the trouble. Please let me know if you need any help resolving the conflicts.

@seisman done in a3a6a05.
seisman: This PR is ready to merge. 

@weiji14 @michaelgrund @willschlitzer @meghanrjones Please see if you have further comments.

Patch:
diff --git a/.github/ISSUE_TEMPLATE/release_checklist.md b/.github/ISSUE_TEMPLATE/release_checklist.md
index f3479658..60c16647 100644
--- a/.github/ISSUE_TEMPLATE/release_checklist.md
+++ b/.github/ISSUE_TEMPLATE/release_checklist.md
@@ -38,7 +38,7 @@ assignees: ''
   - [ ] GMT [forum](https://forum.generic-mapping-tools.org/c/news/)
   - [ ] [Major/Minor releases only] GMT [website](https://github.com/GenericMappingTools/website) (News)
   - [ ] [ResearchGate](https://www.researchgate.net/project/PyGMT-A-Python-interface-for-the-Generic-Mapping-Tools)
-
+  - [ ] [Twitter](https://twitter.com/gmt_dev)
 ---
 
 - [ ] Party :tada: (don't tick before all other checkboxes are ticked!)
diff --git a/pygmt/datasets/earth_relief.py b/pygmt/datasets/earth_relief.py
index 1818b8a8..1a949cc3 100644
--- a/pygmt/datasets/earth_relief.py
+++ b/pygmt/datasets/earth_relief.py
@@ -10,8 +10,8 @@ from pygmt.helpers import kwargs_to_strings
 from pygmt.src import grdcut, which
 
 
-@kwargs_to_strings(region="sequence")
-def load_earth_relief(resolution="01d", region=None, registration=None):
+@kwargs_to_strings(convert_bools=False, region="sequence")
+def load_earth_relief(resolution="01d", region=None, registration=None, use_srtm=False):
     r"""
     Load Earth relief grids (topography and bathymetry) in various resolutions.
 
@@ -50,6 +50,13 @@ def load_earth_relief(resolution="01d", region=None, registration=None):
         a pixel-registered grid is returned unless only the
         gridline-registered grid is available.
 
+    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+V2.1 tiles which have a resolution of 15
+        arc-second (i.e., ``'15s'``). If True, will only load the original
+        land-only SRTM tiles.
+
     Returns
     -------
     grid : :class:`xarray.DataArray`
@@ -73,12 +80,21 @@ def load_earth_relief(resolution="01d", region=None, registration=None):
     >>> grid = load_earth_relief(
     ...     "05m", region=[120, 160, 30, 60], registration="gridline"
     ... )
+    >>> # load the original 3 arc-second land-only SRTM tiles from NASA
+    >>> grid = load_earth_relief(
+    ...     "03s",
+    ...     region=[135, 136, 35, 36],
+    ...     registration="gridline",
+    ...     use_srtm=True,
+    ... )
     """
 
     # 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
     tiled_resolutions = ["05m", "04m", "03m", "02m", "01m", "30s", "15s", "03s", "01s"]
+    # resolutions of original land-only SRTM tiles from NASA
+    land_only_srtm_resolutions = ["03s", "01s"]
 
     if registration in ("pixel", "gridline", None):
         # If None, let GMT decide on Pixel/Gridline type
@@ -103,21 +119,25 @@ def load_earth_relief(resolution="01d", region=None, registration=None):
             f"resolution '{resolution}' is not supported."
         )
 
+    # 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_"
+
     # different ways to load tiled and non-tiled earth relief data
     # 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:
+        if resolution not in non_tiled_resolutions:
             raise GMTInvalidInput(
                 f"'region' is required 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
     else:
-        grid = grdcut(f"@earth_relief_{resolution}{reg}", region=region)
+        grid = grdcut(f"@{earth_relief_prefix}{resolution}{reg}", region=region)
 
     # Add some metadata to the grid
     grid.name = "elevation"


Test patch:
diff --git a/pygmt/tests/test_datasets_earth_relief.py b/pygmt/tests/test_datasets_earth_relief.py
index 297cd3bb..f735a267 100644
--- a/pygmt/tests/test_datasets_earth_relief.py
+++ b/pygmt/tests/test_datasets_earth_relief.py
@@ -83,6 +83,25 @@ def test_earth_relief_05m_without_region():
         load_earth_relief("05m")
 
 
+def test_earth_relief_03s_landonly_srtm():
+    """
+    Test loading original 3 arc-second land-only SRTM tiles.
+    """
+    data = load_earth_relief(
+        "03s", region=[135, 136, 35, 36], registration="gridline", use_srtm=True
+    )
+
+    assert data.coords["lat"].data.min() == 35.0
+    assert data.coords["lat"].data.max() == 36.0
+    assert data.coords["lon"].data.min() == 135.0
+    assert data.coords["lon"].data.max() == 136.0
+    # data.data.min() == -305.51846 if use_srtm is False.
+    assert data.data.min() == -6.0
+    assert data.data.max() == 1191.0
+    assert data.sizes["lat"] == 1201
+    assert data.sizes["lon"] == 1201
+
+
 def test_earth_relief_incorrect_registration():
     """
     Test loading earth relief with incorrect registration type.