diff --git a/doc/api/index.rst b/doc/api/index.rst
index 536c9ced9d3..3536e83f640 100644
--- a/doc/api/index.rst
+++ b/doc/api/index.rst
@@ -94,6 +94,7 @@ Operations on grids:
     grdfilter
     grdlandmask
     grdgradient
+    grdsample
     grdtrack
     xyz2grd
 
diff --git a/pygmt/__init__.py b/pygmt/__init__.py
index f25e8a7e853..25ad96a2127 100644
--- a/pygmt/__init__.py
+++ b/pygmt/__init__.py
@@ -40,6 +40,7 @@
     grdgradient,
     grdinfo,
     grdlandmask,
+    grdsample,
     grdtrack,
     info,
     makecpt,
diff --git a/pygmt/src/__init__.py b/pygmt/src/__init__.py
index a2ced5dbada..2e1f8b56186 100644
--- a/pygmt/src/__init__.py
+++ b/pygmt/src/__init__.py
@@ -19,6 +19,7 @@
 from pygmt.src.grdimage import grdimage
 from pygmt.src.grdinfo import grdinfo
 from pygmt.src.grdlandmask import grdlandmask
+from pygmt.src.grdsample import grdsample
 from pygmt.src.grdtrack import grdtrack
 from pygmt.src.grdview import grdview
 from pygmt.src.histogram import histogram
diff --git a/pygmt/src/grdsample.py b/pygmt/src/grdsample.py
new file mode 100644
index 00000000000..b1c807ccf94
--- /dev/null
+++ b/pygmt/src/grdsample.py
@@ -0,0 +1,95 @@
+"""
+grdsample - Resample a grid onto a new lattice
+"""
+
+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",
+    J="projection",
+    I="spacing",
+    R="region",
+    T="translate",
+    V="verbose",
+    f="coltypes",
+    n="interpolation",
+    r="registration",
+    x="cores",
+)
+@kwargs_to_strings(I="sequence", R="sequence")
+def grdsample(grid, **kwargs):
+    r"""
+    Change the registration, spacing, or nodes in a grid file.
+
+    This reads a grid file and interpolates it to create a new grid
+    file. It can change the registration with ``translate`` or
+    ``registration``, change the grid-spacing or number of nodes with
+    ``spacing``, and set a new sub-region using ``region``. A bicubic
+    [Default], bilinear, B-spline or nearest-neighbor interpolation is set
+    with ``interpolation``.
+
+    When ``region`` is omitted, the output grid will cover the same region as
+    the input grid. When ``spacing`` is omitted, the grid spacing of the
+    output grid will be the same as the input grid. Either ``registration`` or
+    ``translate`` can be used to change the grid registration. When omitted,
+    the output grid will have the same registration as the input grid.
+
+    {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.
+    {I}
+    {R}
+    translate : bool
+        Translate between grid and pixel registration; if the input is
+        grid-registered, the output will be pixel-registered and vice-versa.
+    registration : str or bool
+        [**g**\ |\ **p**\ ].
+        Set registration to **g**\ ridline or **p**\ ixel.
+    {V}
+    {f}
+    {n}
+    {x}
+
+    Returns
+    -------
+    ret: xarray.DataArray or None
+        Return type depends on whether the ``outgrid`` parameter is set:
+
+        - :class:`xarray.DataArray` if ``outgrid`` is not set
+        - None if ``outgrid`` is set (grid output will be stored in file set by
+          ``outgrid``)
+    """
+    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("grdsample", 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
diff --git a/pygmt/tests/test_grdsample.py b/pygmt/tests/test_grdsample.py
new file mode 100644
index 00000000000..5222e7cf987
--- /dev/null
+++ b/pygmt/tests/test_grdsample.py
@@ -0,0 +1,43 @@
+"""
+Tests for grdsample.
+"""
+import os
+
+import pytest
+from pygmt import grdinfo, grdsample
+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], registration="pixel"
+    )
+
+
+def test_grdsample_file_out(grid):
+    """
+    grdsample with an outgrid set and the spacing is changed.
+    """
+    with GMTTempFile(suffix=".nc") as tmpfile:
+        result = grdsample(grid=grid, outgrid=tmpfile.name, spacing=[1, 0.5])
+        assert result is None  # return value is None
+        assert os.path.exists(path=tmpfile.name)  # check that outgrid exists
+        result = grdinfo(tmpfile.name, per_column=True).strip().split()
+        assert float(result[6]) == 1  # x-increment
+        assert float(result[7]) == 0.5  # y-increment
+
+
+def test_grdsample_no_outgrid(grid):
+    """
+    Test grdsample with no set outgrid and applying registration changes.
+    """
+    assert grid.gmt.registration == 1  # Pixel registration
+    translated_grid = grdsample(grid=grid, translate=True)
+    assert translated_grid.gmt.registration == 0  # Gridline registration
+    registration_grid = grdsample(grid=translated_grid, registration="p")
+    assert registration_grid.gmt.registration == 1  # Pixel registration