diff --git a/doc/api/index.rst b/doc/api/index.rst
index 591784e82c9..dff153d1aad 100644
--- a/doc/api/index.rst
+++ b/doc/api/index.rst
@@ -129,6 +129,7 @@ Operations on tabular data
     blockmedian
     blockmode
     filter1d
+    fitcircle
     nearneighbor
     project
     select
diff --git a/pygmt/__init__.py b/pygmt/__init__.py
index 8adecdf429a..a827baa277d 100644
--- a/pygmt/__init__.py
+++ b/pygmt/__init__.py
@@ -35,6 +35,7 @@
     config,
     dimfilter,
     filter1d,
+    fitcircle,
     grd2cpt,
     grd2xyz,
     grdclip,
diff --git a/pygmt/helpers/caching.py b/pygmt/helpers/caching.py
index e315a0eddcc..51c886f4fee 100644
--- a/pygmt/helpers/caching.py
+++ b/pygmt/helpers/caching.py
@@ -116,6 +116,7 @@ def cache_data() -> None:
         "@RidgeTest.prj",
         "@RidgeTest.shp",
         "@RidgeTest.shx",
+        "@sat_03.txt",
         "@SOEST_block4.png",
         "@Table_5_11.txt",
         "@Table_5_11_mean.xyz",
diff --git a/pygmt/src/__init__.py b/pygmt/src/__init__.py
index 7f77088edb9..4e500b5b7be 100644
--- a/pygmt/src/__init__.py
+++ b/pygmt/src/__init__.py
@@ -9,6 +9,7 @@
 from pygmt.src.config import config
 from pygmt.src.dimfilter import dimfilter
 from pygmt.src.filter1d import filter1d
+from pygmt.src.fitcircle import fitcircle
 from pygmt.src.grd2cpt import grd2cpt
 from pygmt.src.grd2xyz import grd2xyz
 from pygmt.src.grdclip import grdclip
diff --git a/pygmt/src/fitcircle.py b/pygmt/src/fitcircle.py
new file mode 100644
index 00000000000..c451d15d4c6
--- /dev/null
+++ b/pygmt/src/fitcircle.py
@@ -0,0 +1,141 @@
+"""
+fitcircle - Find mean position and great [or small] circle fit to points on
+sphere.
+"""
+
+from typing import Literal
+
+import numpy as np
+import pandas as pd
+from pygmt._typing import PathLike, TableLike
+from pygmt.alias import Alias, AliasSystem
+from pygmt.clib import Session
+from pygmt.exceptions import GMTParameterError, GMTValueError
+from pygmt.helpers import build_arg_list, fmt_docstring, validate_output_table_type
+
+
+@fmt_docstring
+def fitcircle(
+    data: PathLike | TableLike | None = None,
+    x=None,
+    y=None,
+    output_type: Literal["pandas", "numpy", "file"] = "pandas",
+    outfile: PathLike | None = None,
+    norm: Literal["absolutes", "squares", "both"] | None = None,
+    small_circle: bool | float = False,
+    verbose: Literal["quiet", "error", "warning", "timing", "info", "compat", "debug"]
+    | bool = False,
+    **kwargs,
+) -> pd.DataFrame | np.ndarray | None:
+    r"""
+    Find mean position and great [or small] circle fit to points on sphere.
+
+    **fitcircle** reads (longitude, latitude) or (latitude, longitude) values from the
+    first two columns of the input data. These are converted to Cartesian
+    three-vectors on the unit sphere. Then two locations are found: the mean
+    of the input positions, and the pole to the great circle which best fits
+    the input positions. The user may choose one or both of two possible
+    solutions to this problem. When the data are closely grouped along a
+    great circle both solutions are similar. If the data have large
+    dispersion, the pole to the great circle will be less well determined
+    than the mean. Compare both solutions as a qualitative check.
+
+    Setting ``norm`` to ``"absolutes"`` approximates the minimization of the
+    sum of absolute values of cosines of angular distances. This solution
+    finds the mean position as the Fisher average of the data, and the pole
+    position as the Fisher average of the cross-products between the mean
+    and the data. Averaging cross-products gives weight to points in
+    proportion to their distance from the mean, analogous to the "leverage"
+    of distant points in linear regression in the plane.
+
+    Setting ``norm`` to ``"squares"`` approximates the minimization of the
+    sum of squares of cosines of angular distances. It creates a 3 by 3
+    matrix of sums of squares of components of the data vectors. The
+    eigenvectors of this matrix give the mean and pole locations. This
+    method may be more subject to roundoff errors when there are thousands
+    of data. The pole is given by the eigenvector corresponding to the
+    smallest eigenvalue; it is the least-well represented factor in the data
+    and is not easily estimated by either method.
+
+    Takes a matrix, (x, y) pairs, or a file name as input.
+
+    Must provide either ``data`` or ``x`` and ``y``.
+
+    Full GMT docs at :gmt-docs:`fitcircle.html`.
+
+    $aliases
+       - V = verbose
+
+    Parameters
+    ----------
+    data
+        Pass in (longitude, latitude) or (latitude, longitude) values by
+        providing a file name to an ASCII data table, a 2-D
+        $table_classes.
+    x/y : 1-D arrays
+        Arrays of x and y coordinates of the data points.
+    $output_type
+    $outfile
+    norm
+        Specify the desired norm. Use ``"absolutes"`` or ``"squares"`` to
+        select a single solution, or ``"both"`` to see both solutions. Note
+        that ``output_type="pandas"`` is not supported when ``norm`` is
+        ``"both"``; use ``output_type="numpy"`` or ``output_type="file"``
+        instead.
+    small_circle : bool or float
+        Attempt to fit a small circle instead of a great circle. The pole
+        will be constrained to lie on the great circle connecting the pole
+        of the best-fit great circle and the mean location of the data.
+        Optionally append the desired fixed latitude of the small circle
+        [Default will determine the optimal latitude].
+    $verbose
+
+    Returns
+    -------
+    ret
+        Return type depends on ``outfile`` and ``output_type``:
+
+        - ``None`` if ``outfile`` is set (output will be stored in the file set by
+          ``outfile``)
+        - :class:`pandas.DataFrame` or :class:`numpy.ndarray` if ``outfile`` is not set
+          (depends on ``output_type``)
+    """
+    if norm is None:
+        raise GMTParameterError(required="norm")
+
+    output_type = validate_output_table_type(output_type, outfile=outfile)
+    if output_type == "pandas" and norm == "both":
+        raise GMTValueError(
+            norm,
+            description="value for parameter 'norm'",
+            reason=(
+                "Pandas output is not supported when 'norm' is set to 'both' "
+                "since both solutions are stacked in the same rows. Use "
+                "output_type='numpy' or output_type='file' instead."
+            ),
+        )
+
+    aliasdict = AliasSystem(
+        L=Alias(norm, name="norm", mapping={"absolutes": 1, "squares": 2, "both": 3}),
+        S=Alias(small_circle, name="small_circle"),
+    ).add_common(
+        V=verbose,
+    )
+    aliasdict.merge(kwargs)
+
+    with Session() as lib:
+        with (
+            lib.virtualfile_in(
+                check_kind="vector", data=data, x=x, y=y, mincols=2
+            ) as vintbl,
+            lib.virtualfile_out(kind="dataset", fname=outfile) as vouttbl,
+        ):
+            lib.call_module(
+                module="fitcircle",
+                args=build_arg_list(aliasdict, infile=vintbl, outfile=vouttbl),
+            )
+        return lib.virtualfile_to_dataset(
+            vfname=vouttbl,
+            output_type=output_type,
+            column_names=["longitude", "latitude", "method"],
+        )
diff --git a/pygmt/tests/test_fitcircle.py b/pygmt/tests/test_fitcircle.py
new file mode 100644
index 00000000000..b4cbd41807b
--- /dev/null
+++ b/pygmt/tests/test_fitcircle.py
@@ -0,0 +1,138 @@
+"""
+Test pygmt.fitcircle.
+"""
+
+from pathlib import Path
+
+import numpy as np
+import numpy.testing as npt
+import pandas as pd
+import pytest
+from pygmt import fitcircle
+from pygmt.exceptions import GMTParameterError, GMTValueError
+from pygmt.helpers import GMTTempFile
+from pygmt.src import which
+
+
+@pytest.fixture(scope="module", name="data")
+def fixture_data():
+    """
+    Load the sample data from the @sat_03 remote file.
+    """
+    fname = which("@sat_03.txt", download="c")
+    return pd.read_csv(
+        fname, header=None, skiprows=1, sep="\t", names=["longitude", "latitude", "z"]
+    )
+
+
+@pytest.mark.benchmark
+def test_fitcircle_no_outfile(data):
+    """
+    Test fitcircle with no set outfile.
+    """
+    result = fitcircle(data=data, norm="squares")
+    assert isinstance(result, pd.DataFrame)
+    assert result.shape == (4, 3)
+    # Test longitude results
+    npt.assert_allclose(result.longitude.min(), 52.7449849947)
+    npt.assert_allclose(result.longitude.max(), 330.243649573)
+    # Test latitude results
+    npt.assert_allclose(result.latitude.min(), -21.2046833116)
+    npt.assert_allclose(result.latitude.max(), 21.2046833116)
+
+
+def test_fitcircle_file_output(data):
+    """
+    Test that fitcircle returns a file output when it is specified.
+    """
+    with GMTTempFile(suffix=".txt") as tmpfile:
+        result = fitcircle(
+            data=data, norm="both", outfile=tmpfile.name, output_type="file"
+        )
+        assert result is None  # return value is None
+        assert Path(tmpfile.name).stat().st_size > 0  # check that outfile exists
+
+
+def test_fitcircle_invalid_format(data):
+    """
+    Test that fitcircle fails with an incorrect format for output_type.
+    """
+    with pytest.raises(GMTValueError):
+        fitcircle(data=data, norm="both", output_type="a")
+
+
+def test_fitcircle_no_norm(data):
+    """
+    Test that fitcircle fails when the required "norm" parameter is missing.
+    """
+    with pytest.raises(GMTParameterError):
+        fitcircle(data=data)
+
+
+def test_fitcircle_no_outfile_specified(data):
+    """
+    Test that fitcircle fails when output_type is set to "file" but no outfile
+    is specified.
+    """
+    with pytest.raises(GMTParameterError):
+        fitcircle(data=data, norm="both", output_type="file")
+
+
+def test_fitcircle_outfile_incorrect_output_type(data):
+    """
+    Test that fitcircle raises a warning when an outfile filename is set but the
+    output_type is not set to "file".
+    """
+    with GMTTempFile(suffix=".txt") as tmpfile:
+        with pytest.warns(RuntimeWarning) as record:
+            result = fitcircle(
+                data=data, norm="both", outfile=tmpfile.name, output_type="numpy"
+            )
+        assert len(record) == 1  # check that only one warning was raised
+        assert result is None  # return value is None
+        assert Path(tmpfile.name).stat().st_size > 0  # check that outfile exists
+
+
+def test_fitcircle_format(data):
+    """
+    Test that correct formats are returned.
+    """
+    circle_default = fitcircle(data=data, norm="squares")
+    assert isinstance(circle_default, pd.DataFrame)
+    circle_array = fitcircle(data=data, norm="squares", output_type="numpy")
+    assert isinstance(circle_array, np.ndarray)
+    circle_df = fitcircle(data=data, norm="squares", output_type="pandas")
+    assert isinstance(circle_df, pd.DataFrame)
+
+
+def test_fitcircle_pandas_unsupported_for_both_norms(data):
+    """
+    Test that fitcircle raises an exception when output_type is "pandas" (the
+    default) and norm is "both", since the two solutions are stacked in the
+    same rows and can't be represented as a single pandas.DataFrame.
+    """
+    with pytest.raises(GMTValueError):
+        fitcircle(data=data, norm="both")
+    with pytest.raises(GMTValueError):
+        fitcircle(data=data, norm="both", output_type="pandas")
+    result = fitcircle(data=data, norm="both", output_type="numpy")
+    assert isinstance(result, np.ndarray)
+
+
+def test_fitcircle_small_circle(data):
+    """
+    Test that fitcircle can fit a small circle instead of a great circle.
+    """
+    result = fitcircle(data=data, norm="squares", small_circle=True)
+    assert isinstance(result, pd.DataFrame)
+    assert result.shape == (5, 3)
+    assert "Small Circle Pole" in result.method.iloc[-1]
+
+
+def test_fitcircle_input_xy(data):
+    """
+    Run fitcircle by passing in x/y as input.
+    """
+    result = fitcircle(x=data.longitude, y=data.latitude, norm="squares")
+    assert isinstance(result, pd.DataFrame)
+    assert result.shape == (4, 3)