Wrap grdclip function (GenericMappingTools#1261)

Wrap the gmt function grdclip and add tests for it.

Co-authored-by: Wei Ji <[email protected]>
Co-authored-by: Dongdong Tian <[email protected]>
Co-authored-by: Wei Ji <[email protected]>

doc/api/index.rst

@@ -88,6 +88,7 @@ Operations on grids:
 .. autosummary::
     :toctree: generated
 
+    grdclip
     grdcut
     grdfilter
     grdtrack

pygmt/__init__.py

@@ -33,6 +33,7 @@
     blockmedian,
     config,
     grd2cpt,
+    grdclip,
     grdcut,
     grdfilter,
     grdinfo,

pygmt/src/__init__.py

@@ -10,6 +10,7 @@
 from pygmt.src.config import config
 from pygmt.src.contour import contour
 from pygmt.src.grd2cpt import grd2cpt
+from pygmt.src.grdclip import grdclip
 from pygmt.src.grdcontour import grdcontour
 from pygmt.src.grdcut import grdcut
 from pygmt.src.grdfilter import grdfilter

pygmt/src/grdclip.py

@@ -0,0 +1,89 @@
+"""
+grdclip - Change the range and extremes of grid values.
+"""
+
+import xarray as xr
+from pygmt.clib import Session
+from pygmt.helpers import (
+    GMTTempFile,
+    build_arg_string,
+    fmt_docstring,
+    kwargs_to_strings,
+    use_alias,
+)
+
+
+@fmt_docstring
+@use_alias(
+    G="outgrid",
+    R="region",
+    Sa="above",
+    Sb="below",
+    Si="between",
+    Sr="new",
+    V="verbose",
+)
+@kwargs_to_strings(
+    R="sequence",
+    Sa="sequence",
+    Sb="sequence",
+    Si="sequence",
+    Sr="sequence",
+)
+def grdclip(grid, **kwargs):
+    r"""
+    Sets values in a grid that meet certain criteria to a new value.
+
+    Produce a clipped ``outgrid`` or :class:`xarray.DataArray` version of the
+    input ``grid`` file.
+
+    The parameters ``above`` and ``below`` allow for a given value to be set
+    for values above or below a set amount, respectively. This allows for
+    extreme values in a grid, such as points below a certain depth when
+    plotting Earth relief, to all be set to the same value.
+
+    Full option list at :gmt-docs:`grdclip.html`
+
+    {aliases}
+
+    Parameters
+    ----------
+    grid : str or xarray.DataArray
+        The file name of the input grid or the grid loaded as a DataArray.
+    outgrid : str or None
+        The name of the output netCDF file with extension .nc to store the grid
+        in.
+    {R}
+    above : str or list or tuple
+        [*high*, *above*].
+        Set all data[i] > *high* to *above*.
+    below : str or list or tuple
+        [*low*, *below*].
+        Set all data[i] < *low* to *below*.
+    between : str or list or tuple
+        [*low*, *high*, *between*].
+        Set all data[i] >= *low* and <= *high* to *between*.
+    new : str or list or tuple
+        [*old*, *new*].
+        Set all data[i] == *old* to *new*. This is mostly useful when
+        your data are known to be integer values.
+    {V}
+    """
+    with GMTTempFile(suffix=".nc") as tmpfile:
+        with Session() as lib:
+            file_context = lib.virtualfile_from_data(check_kind="raster", data=grid)
+            with file_context as infile:
+                if "G" not in kwargs.keys():  # if outgrid is unset, output to tempfile
+                    kwargs.update({"G": tmpfile.name})
+                outgrid = kwargs["G"]
+                arg_str = " ".join([infile, build_arg_string(kwargs)])
+                lib.call_module("grdclip", arg_str)
+
+        if outgrid == tmpfile.name:  # if user did not set outgrid, return DataArray
+            with xr.open_dataarray(outgrid) as dataarray:
+                result = dataarray.load()
+                _ = result.gmt  # load GMTDataArray accessor information
+        else:
+            result = None  # if user sets an outgrid, return None
+
+        return result

pygmt/tests/test_grdclip.py

@@ -0,0 +1,47 @@
+"""
+Tests for grdclip.
+"""
+import os
+
+import numpy.testing as npt
+import pytest
+from pygmt import grdclip, grdinfo
+from pygmt.datasets import load_earth_relief
+from pygmt.helpers import GMTTempFile
+
+
+@pytest.fixture(scope="module", name="grid")
+def fixture_grid():
+    """
+    Load the grid data from the sample earth_relief file.
+    """
+    return load_earth_relief(resolution="01d", region=[-5, 5, -5, 5])
+
+
+def test_grdclip_outgrid(grid):
+    """
+    Test the below and above parameters for grdclip and creates a test outgrid.
+    """
+    with GMTTempFile(suffix=".nc") as tmpfile:
+        result = grdclip(
+            grid=grid, outgrid=tmpfile.name, below=[-1500, -1800], above=[-200, 40]
+        )
+        assert result is None  # return value is None
+        assert os.path.exists(path=tmpfile.name)  # check that outgrid exists
+        result = (
+            grdinfo(grid=tmpfile.name, force_scan=0, per_column="n").strip().split()
+        )
+    assert int(result[4]) == -1800  # minimum value
+    assert int(result[5]) == 40  # maximum value
+
+
+def test_grdclip_no_outgrid(grid):
+    """
+    Test the below and above parameters for grdclip with no set outgrid.
+    """
+    temp_grid = grdclip(grid=grid, below=[-1500, -1800], above=[-200, 40])
+    assert temp_grid.dims == ("lat", "lon")
+    assert temp_grid.gmt.gtype == 1  # Geographic grid
+    assert temp_grid.gmt.registration == 1  # Pixel registration
+    npt.assert_allclose(temp_grid.min(), -1800)
+    npt.assert_allclose(temp_grid.max(), 40)