diff --git a/doc/api/index.rst b/doc/api/index.rst
index 69685060086..c6f58cb9468 100644
--- a/doc/api/index.rst
+++ b/doc/api/index.rst
@@ -83,6 +83,7 @@ Operations on tabular data:
     blockmedian
     blockmode
     nearneighbor
+    sph2grd
     surface
 
 Operations on grids:
diff --git a/pygmt/__init__.py b/pygmt/__init__.py
index 7a1a1079f37..14e4bad6c4c 100644
--- a/pygmt/__init__.py
+++ b/pygmt/__init__.py
@@ -49,6 +49,7 @@
     info,
     makecpt,
     nearneighbor,
+    sph2grd,
     sphdistance,
     surface,
     which,
diff --git a/pygmt/helpers/testing.py b/pygmt/helpers/testing.py
index 1d4e5ec4159..98b866ae09a 100644
--- a/pygmt/helpers/testing.py
+++ b/pygmt/helpers/testing.py
@@ -165,6 +165,7 @@ def download_test_data():
         "@N37W120.earth_relief_03s_g.nc",
         "@N00W090.earth_relief_03m_p.nc",
         # Other cache files
+        "@EGM96_to_36.txt",
         "@fractures_06.txt",
         "@hotspots.txt",
         "@ridge.txt",
diff --git a/pygmt/src/__init__.py b/pygmt/src/__init__.py
index a3dded681eb..95b82387192 100644
--- a/pygmt/src/__init__.py
+++ b/pygmt/src/__init__.py
@@ -37,6 +37,7 @@
 from pygmt.src.plot3d import plot3d
 from pygmt.src.rose import rose
 from pygmt.src.solar import solar
+from pygmt.src.sph2grd import sph2grd
 from pygmt.src.sphdistance import sphdistance
 from pygmt.src.subplot import set_panel, subplot
 from pygmt.src.surface import surface
diff --git a/pygmt/src/sph2grd.py b/pygmt/src/sph2grd.py
new file mode 100644
index 00000000000..7070b390a93
--- /dev/null
+++ b/pygmt/src/sph2grd.py
@@ -0,0 +1,77 @@
+"""
+sph2grd - Compute grid from spherical harmonic coefficients
+"""
+from pygmt.clib import Session
+from pygmt.helpers import (
+    GMTTempFile,
+    build_arg_string,
+    fmt_docstring,
+    kwargs_to_strings,
+    use_alias,
+)
+from pygmt.io import load_dataarray
+
+
+@fmt_docstring
+@use_alias(
+    G="outgrid",
+    I="spacing",
+    R="region",
+    V="verbose",
+    b="binary",
+    h="header",
+    i="incols",
+    r="registration",
+    x="cores",
+)
+@kwargs_to_strings(I="sequence", R="sequence", i="sequence_comma")
+def sph2grd(data, **kwargs):
+    r"""
+    Create spherical grid files in tension of data.
+
+    Reads a spherical harmonics coefficient table with records of L, M,
+    C[L,M], S[L,M] and evaluates the spherical harmonic model on the
+    specified grid.
+
+    Full option list at :gmt-docs:`sph2grd.html`
+
+    {aliases}
+
+    Parameters
+    ----------
+    data : str or {table-like}
+        Pass in data with L, M, C[L,M], S[L,M] values by
+        providing a file name to an ASCII data table, a 2D
+        {table-classes}.
+    outgrid : str or None
+        The name of the output netCDF file with extension .nc to store the grid
+        in.
+    {I}
+    {R}
+    {V}
+    {b}
+    {h}
+    {i}
+    {r}
+    {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="vector", data=data)
+            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("sph2grd", arg_str)
+
+        return load_dataarray(outgrid) if outgrid == tmpfile.name else None
diff --git a/pygmt/tests/test_sph2grd.py b/pygmt/tests/test_sph2grd.py
new file mode 100644
index 00000000000..316e6d52408
--- /dev/null
+++ b/pygmt/tests/test_sph2grd.py
@@ -0,0 +1,34 @@
+"""
+Tests for sph2grd.
+"""
+import os
+
+import numpy.testing as npt
+from pygmt import sph2grd
+from pygmt.helpers import GMTTempFile
+
+
+def test_sph2grd_outgrid():
+    """
+    Test sph2grd with a set outgrid.
+    """
+    with GMTTempFile(suffix=".nc") as tmpfile:
+        result = sph2grd(
+            data="@EGM96_to_36.txt", outgrid=tmpfile.name, spacing=1, region="g"
+        )
+        assert result is None  # return value is None
+        assert os.path.exists(path=tmpfile.name)  # check that outgrid exists
+
+
+def test_sph2grd_no_outgrid():
+    """
+    Test sph2grd with no set outgrid.
+    """
+    temp_grid = sph2grd(data="@EGM96_to_36.txt", spacing=1, region="g")
+    assert temp_grid.dims == ("y", "x")
+    assert temp_grid.gmt.gtype == 0  # Cartesian grid
+    assert temp_grid.gmt.registration == 0  # Gridline registration
+    npt.assert_allclose(temp_grid.max(), 0.00021961, rtol=1e-4)
+    npt.assert_allclose(temp_grid.min(), -0.0004326, rtol=1e-4)
+    npt.assert_allclose(temp_grid.median(), -0.00010894, rtol=1e-4)
+    npt.assert_allclose(temp_grid.mean(), -0.00010968, rtol=1e-4)