diff --git a/doc/api/index.rst b/doc/api/index.rst
index 40ab49ac430..536c9ced9d3 100644
--- a/doc/api/index.rst
+++ b/doc/api/index.rst
@@ -95,6 +95,7 @@ Operations on grids:
     grdlandmask
     grdgradient
     grdtrack
+    xyz2grd
 
 Crossover analysis with x2sys:
 
diff --git a/pygmt/__init__.py b/pygmt/__init__.py
index 972f6474286..f25e8a7e853 100644
--- a/pygmt/__init__.py
+++ b/pygmt/__init__.py
@@ -47,6 +47,7 @@
     which,
     x2sys_cross,
     x2sys_init,
+    xyz2grd,
 )
 
 # Get semantic version through setuptools-scm
diff --git a/pygmt/src/__init__.py b/pygmt/src/__init__.py
index f8bae667f02..a2ced5dbada 100644
--- a/pygmt/src/__init__.py
+++ b/pygmt/src/__init__.py
@@ -41,3 +41,4 @@
 from pygmt.src.wiggle import wiggle
 from pygmt.src.x2sys_cross import x2sys_cross
 from pygmt.src.x2sys_init import x2sys_init
+from pygmt.src.xyz2grd import xyz2grd
diff --git a/pygmt/src/xyz2grd.py b/pygmt/src/xyz2grd.py
new file mode 100644
index 00000000000..ef40abc889d
--- /dev/null
+++ b/pygmt/src/xyz2grd.py
@@ -0,0 +1,74 @@
+"""
+xyz2grd - Convert data table to a grid.
+"""
+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",
+    I="spacing",
+    R="region",
+    V="verbose",
+)
+@kwargs_to_strings(R="sequence")
+def xyz2grd(table, **kwargs):
+    """
+    Create a grid file from table data.
+
+    xyz2grd reads one or more z or xyz tables and creates a binary grid file.
+    xyz2grd will report if some of the nodes are not filled in with data. Such
+    unconstrained nodes are set to a value specified by the user [Default is
+    NaN]. Nodes with more than one value will be set to the mean value.
+
+    Full option list at :gmt-docs:`xyz2grd.html`
+
+    Parameters
+    ----------
+    table : str or {table-like}
+        Pass in either a file name to an ASCII data table, a 1D/2D
+        {table-classes}.
+
+    outgrid : str or None
+        Optional. The name of the output netCDF file with extension .nc to
+        store the grid in.
+    {I}
+    {R}
+    {V}
+
+    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="vector", data=table)
+            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 = build_arg_string(kwargs)
+                arg_str = " ".join([infile, arg_str])
+                lib.call_module("xyz2grd", 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_xyz2grd.py b/pygmt/tests/test_xyz2grd.py
new file mode 100644
index 00000000000..530c6b0905a
--- /dev/null
+++ b/pygmt/tests/test_xyz2grd.py
@@ -0,0 +1,69 @@
+"""
+Tests for xyz2grd.
+"""
+import os
+
+import numpy as np
+import pytest
+import xarray as xr
+from pygmt import grdinfo, xyz2grd
+from pygmt.datasets import load_sample_bathymetry
+from pygmt.helpers import GMTTempFile
+
+
+@pytest.fixture(scope="module", name="ship_data")
+def fixture_ship_data():
+    """
+    Load the data from the sample bathymetry dataset.
+    """
+    return load_sample_bathymetry()
+
+
+def test_xyz2grd_input_file():
+    """
+    Run xyz2grd by passing in a filename.
+    """
+    output = xyz2grd(table="@tut_ship.xyz", spacing=5, region=[245, 255, 20, 30])
+    assert isinstance(output, xr.DataArray)
+    assert output.gmt.registration == 0  # Gridline registration
+    assert output.gmt.gtype == 0  # Cartesian type
+    return output
+
+
+def test_xyz2grd_input_array(ship_data):
+    """
+    Run xyz2grd by passing in a numpy array.
+    """
+    output = xyz2grd(table=np.array(ship_data), spacing=5, region=[245, 255, 20, 30])
+    assert isinstance(output, xr.DataArray)
+    assert output.gmt.registration == 0  # Gridline registration
+    assert output.gmt.gtype == 0  # Cartesian type
+    return output
+
+
+def test_xyz2grd_input_df(ship_data):
+    """
+    Run xyz2grd by passing in a data frame.
+    """
+    output = xyz2grd(table=ship_data, spacing=5, region=[245, 255, 20, 30])
+    assert isinstance(output, xr.DataArray)
+    assert output.gmt.registration == 0  # Gridline registration
+    assert output.gmt.gtype == 0  # Cartesian type
+    return output
+
+
+def test_xyz2grd_input_array_file_out(ship_data):
+    """
+    Run xyz2grd by passing in a numpy array and set an outgrid file.
+    """
+    with GMTTempFile(suffix=".nc") as tmpfile:
+        result = xyz2grd(
+            table=np.array(ship_data),
+            spacing=5,
+            region=[245, 255, 20, 30],
+            outgrid=tmpfile.name,
+        )
+        assert result is None  # return value is None
+        assert os.path.exists(path=tmpfile.name)
+        result = grdinfo(tmpfile.name, per_column=True).strip()
+        assert result == "245 255 20 30 -3651.06079102 -352.379486084 5 5 3 3 0 0"