diff --git a/doc/api/index.rst b/doc/api/index.rst
index 75dedbf737d..197362503c7 100644
--- a/doc/api/index.rst
+++ b/doc/api/index.rst
@@ -55,6 +55,7 @@ Plotting tabular data
     Figure.plot
     Figure.plot3d
     Figure.rose
+    Figure.sac
     Figure.ternary
     Figure.velo
     Figure.wiggle
diff --git a/pygmt/figure.py b/pygmt/figure.py
index 34c345bb707..598157cd019 100644
--- a/pygmt/figure.py
+++ b/pygmt/figure.py
@@ -33,6 +33,7 @@
 from pygmt.src.psconvert import psconvert as _psconvert
 from pygmt.src.pygmtlogo import pygmtlogo as _pygmtlogo
 from pygmt.src.rose import rose as _rose
+from pygmt.src.sac import sac as _sac
 from pygmt.src.scalebar import scalebar as _scalebar
 from pygmt.src.shift_origin import shift_origin as _shift_origin
 from pygmt.src.solar import solar as _solar
@@ -469,6 +470,7 @@ def _repr_html_(self) -> str:
     psconvert = _psconvert
     pygmtlogo = _pygmtlogo
     rose = _rose
+    sac = _sac
     scalebar = _scalebar
     set_panel = _set_panel
     shift_origin = _shift_origin
diff --git a/pygmt/src/sac.py b/pygmt/src/sac.py
new file mode 100644
index 00000000000..12f7e313ea8
--- /dev/null
+++ b/pygmt/src/sac.py
@@ -0,0 +1,208 @@
+"""
+sac - Plot seismograms in SAC format.
+"""
+
+from collections.abc import Sequence
+from typing import Literal
+
+from pygmt._typing import PathLike
+from pygmt.alias import Alias, AliasSystem
+from pygmt.clib import Session
+from pygmt.helpers import build_arg_list, fmt_docstring
+from pygmt.params import Axis, Frame
+
+
+@fmt_docstring
+def sac(
+    self,
+    data: PathLike | Sequence[PathLike],
+    pen: str | None = None,
+    time_window: Sequence[float] | bool = False,
+    offset: float | Sequence[float] | None = None,
+    profile: str | None = None,
+    preprocess: str | None = None,
+    fill: str | Sequence[str] | bool = False,
+    amplitude_scale: float | str | Sequence[float | str] | None = None,
+    vertical: bool = False,
+    time_scale: float | str | None = None,
+    reduction_velocity: float | None = None,
+    time_shift: float | None = None,
+    time_reference: Literal[
+        "b",
+        "e",
+        "o",
+        "a",
+        "t0",
+        "t1",
+        "t2",
+        "t3",
+        "t4",
+        "t5",
+        "t6",
+        "t7",
+        "t8",
+        "t9",
+    ]
+    | None = None,
+    projection: str | None = None,
+    region: Sequence[float | str] | str | None = None,
+    frame: Frame | Axis | Literal["none"] | str | Sequence[str] | bool = False,
+    verbose: Literal["quiet", "error", "warning", "timing", "info", "compat", "debug"]
+    | bool = False,
+    panel: int | Sequence[int] | bool = False,
+    perspective: float | Sequence[float] | str | bool = False,
+    transparency: float | None = None,
+    **kwargs,
+):
+    """
+    Plot seismograms in SAC format.
+
+    Reads SAC waveform files and plots seismic traces. Only evenly spaced SAC
+    data is supported.
+
+    Full GMT docs at :gmt-docs:`supplements/seis/sac.html`.
+
+    **Aliases:**
+
+    .. hlist::
+       :columns: 3
+
+       - B = frame
+       - C = time_window
+       - D = offset
+       - E = profile
+       - F = preprocess
+       - G = fill
+       - J = projection
+       - M = amplitude_scale
+       - Q = vertical
+       - R = region
+       - S = time_scale
+       - T = **+r**: reduction_velocity, **+s**: time_shift, **+t**: time_reference
+       - V = verbose
+       - W = pen
+       - c = panel
+       - p = perspective
+       - t = transparency
+
+    Parameters
+    ----------
+    data
+        The SAC waveform file(s) to plot.
+    pen
+        Set pen attributes for all traces [Default is ``"0.25p,black,solid"``].
+    time_window
+        Read and plot seismograms in the time window between *t0* and *t1* only,
+        where *t0* and *t1* are relative to the reference time (see
+        ``time_reference``). If no reference time is set, the reference time in
+        the SAC header is used. Set to ``True`` to determine *t0*/*t1* from the
+        *xmin*/*xmax* of the ``region`` parameter.
+    offset
+        Offset the seismogram positions by the given amounts *dx*[/ *dy*]
+        [Default is no offset]. If *dy* is not given, it is set equal to *dx*.
+    profile
+        Choose the profile type, i.e., the type of the y axis. Use ``a`` for
+        azimuth, ``b`` for back-azimuth, ``k`` for epicentral distance in km,
+        ``d`` for epicentral distance in degrees, ``n`` for trace number (the
+        first trace is numbered *n*, e.g., ``n0``), and ``u`` for user-defined
+        profile (the y positions are determined by the SAC header variable
+        ``usern``, e.g., ``u0``).
+    preprocess
+        Preprocess the data before plotting. Use ``i`` for integral, ``q`` for
+        square, and ``r`` for removing the mean value. The letters can repeat
+        multiple times, and the order controls the processing order, e.g.,
+        ``"rii"`` converts acceleration to displacement.
+    fill
+        Paint the positive or negative portion of the traces. Use ``p``/``n`` to
+        paint the positive/negative portion [Default paints the positive
+        portion], ``+g`` *fill* to set the fill color [Default is ``"black"``],
+        ``+z`` *zero* to define the zero line, and ``+t`` *t0*/*t1* to paint a
+        time window only. Can be repeated to paint the positive and negative
+        portions separately, e.g., ``["p+gblack", "n+gred"]``. Set to ``True``
+        to paint the positive portion with the default fill.
+    amplitude_scale
+        Set the vertical scaling of the traces.
+
+        - If only *size* is given (optionally with a unit), all traces are
+          scaled to the given height on the map.
+        - If *size*/*alpha* is given with a negative *alpha*, all traces use
+          the same scaling factor determined by the first trace, which is
+          scaled to *size*.
+        - If *alpha* is 0, all traces are multiplied by *size* [no unit is
+          allowed].
+        - If *alpha* is positive, all traces are multiplied by *size* times the
+          epicentral distance (in km) raised to the power *alpha*.
+    vertical
+        Plot traces vertically, i.e., the y axis is time and the x axis is
+        amplitude.
+    time_scale
+        Set the time scale in seconds per unit while plotting on geographic
+        plots. Append a unit (``c``, ``i``, or ``p``); if omitted, the unit is
+        controlled by :gmt-term:`PROJ_LENGTH_UNIT`. Use the ``i`` prefix to give
+        the reciprocal scale, i.e., unit per second.
+    reduction_velocity
+        Set the reduction velocity in km/s. The trace times are reduced by
+        *distance*/*reduction_velocity*, where *distance* is the epicentral
+        distance in the SAC header.
+    time_shift
+        Shift all traces by the given number of seconds.
+    time_reference
+        Align all traces along the given time mark. Choose from ``"b"``,
+        ``"e"``, ``"o"``, ``"a"``, or ``"t0"`` to ``"t9"``.
+    $projection
+    $region
+    $frame
+    $verbose
+    $panel
+    $perspective
+    $transparency
+    """
+    aliasdict = AliasSystem(
+        C=Alias(time_window, name="time_window", sep="/", size=2),
+        D=Alias(offset, name="offset", sep="/", size=(1, 2)),
+        E=Alias(profile, name="profile"),
+        F=Alias(preprocess, name="preprocess"),
+        G=Alias(fill, name="fill"),
+        M=Alias(amplitude_scale, name="amplitude_scale", sep="/", size=(1, 2)),
+        Q=Alias(vertical, name="vertical"),
+        S=Alias(time_scale, name="time_scale"),
+        T=[
+            Alias(reduction_velocity, name="reduction_velocity", prefix="+r"),
+            Alias(time_shift, name="time_shift", prefix="+s"),
+            Alias(
+                time_reference,
+                name="time_reference",
+                prefix="+t",
+                mapping={
+                    "b": -5,
+                    "e": -4,
+                    "o": -3,
+                    "a": -2,
+                    "t0": 0,
+                    "t1": 1,
+                    "t2": 2,
+                    "t3": 3,
+                    "t4": 4,
+                    "t5": 5,
+                    "t6": 6,
+                    "t7": 7,
+                    "t8": 8,
+                    "t9": 9,
+                },
+            ),
+        ],
+        W=Alias(pen, name="pen"),
+    ).add_common(
+        B=frame,
+        J=projection,
+        R=region,
+        V=verbose,
+        c=panel,
+        p=perspective,
+        t=transparency,
+    )
+    aliasdict.merge(kwargs)
+
+    self._activate_figure()
+    with Session() as lib:
+        lib.call_module(module="sac", args=build_arg_list(aliasdict, infile=data))
diff --git a/pygmt/tests/baseline/test_sac.png b/pygmt/tests/baseline/test_sac.png
new file mode 100644
index 00000000000..bca144a6e15
--- /dev/null
+++ b/pygmt/tests/baseline/test_sac.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:14e0cd2cb6c55008a3cf76d68ee86ef9ad1e361752ac2673fab5cc17f21e3fd4
+size 24619
diff --git a/pygmt/tests/baseline/test_sac_amplitude_scale.png b/pygmt/tests/baseline/test_sac_amplitude_scale.png
new file mode 100644
index 00000000000..52c797fddfa
--- /dev/null
+++ b/pygmt/tests/baseline/test_sac_amplitude_scale.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:ae374193202586cf79c3ec1bbb1bab0c18810aa8149d6a0a394275ef701e8e0b
+size 17418
diff --git a/pygmt/tests/baseline/test_sac_fill.png b/pygmt/tests/baseline/test_sac_fill.png
new file mode 100644
index 00000000000..bc6a2119143
--- /dev/null
+++ b/pygmt/tests/baseline/test_sac_fill.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:5d32a7ec7e5be1f5eefe4232b764f706297f04914a26301790beb55758290954
+size 21182
diff --git a/pygmt/tests/baseline/test_sac_fill_positive_negative.png b/pygmt/tests/baseline/test_sac_fill_positive_negative.png
new file mode 100644
index 00000000000..3927e01885c
--- /dev/null
+++ b/pygmt/tests/baseline/test_sac_fill_positive_negative.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:9969d33d01bb7f0f54483170ec4cceac3f4f9ede10f1b07355b9c9fa26ccd281
+size 17104
diff --git a/pygmt/tests/baseline/test_sac_multiple_traces.png b/pygmt/tests/baseline/test_sac_multiple_traces.png
new file mode 100644
index 00000000000..8425503c4aa
--- /dev/null
+++ b/pygmt/tests/baseline/test_sac_multiple_traces.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:9784469c87512a1713676e493689988753b63f8fee5287baea673aa87f08e7c0
+size 19847
diff --git a/pygmt/tests/baseline/test_sac_offset.png b/pygmt/tests/baseline/test_sac_offset.png
new file mode 100644
index 00000000000..62f311af33a
--- /dev/null
+++ b/pygmt/tests/baseline/test_sac_offset.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:71322522d7394b7989fa193e0aa853f0568aec95a0303d2c48a4ef9c62d58c5c
+size 19019
diff --git a/pygmt/tests/baseline/test_sac_preprocess.png b/pygmt/tests/baseline/test_sac_preprocess.png
new file mode 100644
index 00000000000..d3f7929e0b0
--- /dev/null
+++ b/pygmt/tests/baseline/test_sac_preprocess.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:47cf61a63fbc8a8be105096a3f7ba1ce4538cd8c67873fec71c4fa4c014e4e6a
+size 20430
diff --git a/pygmt/tests/baseline/test_sac_profile.png b/pygmt/tests/baseline/test_sac_profile.png
new file mode 100644
index 00000000000..39a40da1889
--- /dev/null
+++ b/pygmt/tests/baseline/test_sac_profile.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:b6574526ec089585cf4d8757faacfbed0aefadbd81bfcfb052d3304424f3d6b9
+size 20796
diff --git a/pygmt/tests/baseline/test_sac_time_options.png b/pygmt/tests/baseline/test_sac_time_options.png
new file mode 100644
index 00000000000..d895a7dc14b
--- /dev/null
+++ b/pygmt/tests/baseline/test_sac_time_options.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:7a09482de41d2c8e4c7ed9f90578dd96214ff2e901d5e0b11a7b0f57bcf35eaa
+size 19292
diff --git a/pygmt/tests/baseline/test_sac_time_scale.png b/pygmt/tests/baseline/test_sac_time_scale.png
new file mode 100644
index 00000000000..b265287e0e1
--- /dev/null
+++ b/pygmt/tests/baseline/test_sac_time_scale.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:d0a66bd8a8ccaf4c51303e2ea3e799275bd7a454acdb1ea9874f8a53ace7a949
+size 11664
diff --git a/pygmt/tests/baseline/test_sac_time_window.png b/pygmt/tests/baseline/test_sac_time_window.png
new file mode 100644
index 00000000000..87d517d8a24
--- /dev/null
+++ b/pygmt/tests/baseline/test_sac_time_window.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:27d8b96799814a9cd7897443059344b5144071107ab05bc178e4cda17d1b2bf6
+size 19587
diff --git a/pygmt/tests/baseline/test_sac_time_window_bare.png b/pygmt/tests/baseline/test_sac_time_window_bare.png
new file mode 100644
index 00000000000..3f030c58b3c
--- /dev/null
+++ b/pygmt/tests/baseline/test_sac_time_window_bare.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:a5940df97c2f6e1d30c0ed9fb6ce35ca4d5e8261627a2c2d51df0f6ce968a5e8
+size 20820
diff --git a/pygmt/tests/baseline/test_sac_vertical.png b/pygmt/tests/baseline/test_sac_vertical.png
new file mode 100644
index 00000000000..760dd9d2327
--- /dev/null
+++ b/pygmt/tests/baseline/test_sac_vertical.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:17a4ecd30d97cbc23a0fee4ea84ba34b4b7ca0d713c9653cf2b308d673e397ef
+size 24242
diff --git a/pygmt/tests/data/seis.sac b/pygmt/tests/data/seis.sac
new file mode 100644
index 00000000000..2ed16d049f8
Binary files /dev/null and b/pygmt/tests/data/seis.sac differ
diff --git a/pygmt/tests/test_sac.py b/pygmt/tests/test_sac.py
new file mode 100644
index 00000000000..7980851a140
--- /dev/null
+++ b/pygmt/tests/test_sac.py
@@ -0,0 +1,260 @@
+"""
+Test Figure.sac.
+"""
+
+from pathlib import Path
+
+import pytest
+from packaging.version import Version
+from pygmt import Figure
+from pygmt.clib import __gmt_version__
+
+SAC_DATA = Path(__file__).parent / "data" / "seis.sac"
+
+# TODO(GMT>6.7.0): Remove the xfail marker when the minimum GMT version is 6.7.0.
+XFAIL_GMT_LE_6_6 = pytest.mark.xfail(
+    condition=Version(__gmt_version__) <= Version("6.6.0"),
+    reason="The sac module was added in GMT 6.7.0.",
+)
+
+
+@XFAIL_GMT_LE_6_6
+@pytest.mark.mpl_image_compare
+def test_sac():
+    """
+    Plot a single SAC waveform on a linear time plot.
+    """
+    fig = Figure()
+    fig.sac(
+        data=SAC_DATA,
+        region=[9, 20, -2, 2],
+        projection="X15c/5c",
+        frame=["xaf+lTime (s)", "yaf+lAmplitude", "WSen"],
+        pen="0.5p,red",
+    )
+    return fig
+
+
+@XFAIL_GMT_LE_6_6
+@pytest.mark.mpl_image_compare
+def test_sac_multiple_traces():
+    """
+    Plot multiple SAC waveforms on a trace number profile.
+    """
+    fig = Figure()
+    fig.sac(
+        data=[SAC_DATA, SAC_DATA],
+        # Plot two traces on a trace number profile, i.e., the y positions of the
+        # two traces are 0 and 1, respectively. The data amplitude is ~±1.6, so
+        # the region is set with some margins.
+        region=[9, 20, -0.5, 3],
+        projection="X15c/5c",
+        frame=True,
+    )
+    return fig
+
+
+@XFAIL_GMT_LE_6_6
+@pytest.mark.mpl_image_compare
+def test_sac_time_window():
+    """
+    Plot a SAC waveform in a given time window.
+    """
+    fig = Figure()
+    fig.sac(
+        data=SAC_DATA,
+        time_window=[10, 18],
+        region=[9, 20, -2, 2],
+        projection="X15c/5c",
+        frame=True,
+        pen="0.5p,red",
+    )
+    return fig
+
+
+@XFAIL_GMT_LE_6_6
+@pytest.mark.mpl_image_compare
+def test_sac_time_window_bare():
+    """
+    Plot a SAC waveform in the time window determined by the region.
+    """
+    fig = Figure()
+    fig.sac(
+        data=SAC_DATA,
+        time_window=True,
+        region=[9, 20, -2, 2],
+        projection="X15c/5c",
+        frame=True,
+        pen="0.5p,red",
+    )
+    return fig
+
+
+@XFAIL_GMT_LE_6_6
+@pytest.mark.mpl_image_compare
+def test_sac_fill_positive_negative():
+    """
+    Plot a SAC waveform with the positive and negative portions filled.
+    """
+    fig = Figure()
+    fig.sac(
+        data=SAC_DATA,
+        fill=["p+gblack", "n+gred"],
+        preprocess="r",
+        region=[9, 20, -2, 2],
+        projection="X15c/5c",
+        frame=True,
+    )
+    return fig
+
+
+@XFAIL_GMT_LE_6_6
+@pytest.mark.mpl_image_compare
+def test_sac_offset():
+    """
+    Plot SAC waveforms with an offset.
+    """
+    fig = Figure()
+    fig.sac(
+        data=[SAC_DATA, SAC_DATA],
+        # Offset the two traces (at y positions 0 and 1) by dy=1, so the y range
+        # of the traces becomes 1 and 2, respectively. The data amplitude is
+        # ~±1.6, so the region is set with some margins.
+        offset=[0, 1],
+        region=[9, 20, -1, 4],
+        projection="X15c/5c",
+        frame=True,
+        pen="0.5p,red",
+    )
+    return fig
+
+
+@XFAIL_GMT_LE_6_6
+@pytest.mark.mpl_image_compare
+def test_sac_fill():
+    """
+    Plot a SAC waveform with the positive portion filled.
+    """
+    fig = Figure()
+    fig.sac(
+        data=SAC_DATA,
+        fill="+ggray",
+        region=[9, 20, -2, 2],
+        projection="X15c/5c",
+        frame=True,
+        pen="0.5p,red",
+    )
+    return fig
+
+
+@XFAIL_GMT_LE_6_6
+@pytest.mark.mpl_image_compare
+def test_sac_amplitude_scale():
+    """
+    Plot a SAC waveform with vertical scaling.
+    """
+    fig = Figure()
+    fig.sac(
+        data=SAC_DATA,
+        amplitude_scale="1.5c",
+        region=[9, 20, -4, 4],
+        projection="X15c/5c",
+        frame=True,
+        pen="0.5p,red",
+    )
+    return fig
+
+
+@XFAIL_GMT_LE_6_6
+@pytest.mark.mpl_image_compare
+def test_sac_profile():
+    """
+    Plot a SAC waveform on a trace number profile.
+    """
+    fig = Figure()
+    fig.sac(
+        data=SAC_DATA,
+        profile="n1",
+        region=[9, 20, -1, 3],
+        projection="X15c/5c",
+        frame=True,
+        pen="0.5p,red",
+    )
+    return fig
+
+
+@XFAIL_GMT_LE_6_6
+@pytest.mark.mpl_image_compare
+def test_sac_preprocess():
+    """
+    Plot a SAC waveform with the mean removed.
+    """
+    fig = Figure()
+    fig.sac(
+        data=SAC_DATA,
+        preprocess="r",
+        region=[9, 20, -2, 2],
+        projection="X15c/5c",
+        frame=True,
+        pen="0.5p,red",
+    )
+    return fig
+
+
+@XFAIL_GMT_LE_6_6
+@pytest.mark.mpl_image_compare
+def test_sac_vertical():
+    """
+    Plot a SAC waveform vertically.
+    """
+    fig = Figure()
+    fig.sac(
+        data=SAC_DATA,
+        vertical=True,
+        region=[-2, 2, 9, 20],
+        projection="X5c/15c",
+        frame=True,
+        pen="0.5p,red",
+    )
+    return fig
+
+
+@XFAIL_GMT_LE_6_6
+@pytest.mark.mpl_image_compare
+def test_sac_time_scale():
+    """
+    Plot a SAC waveform on a geographic map with a time scale.
+    """
+    fig = Figure()
+    # The station is at (stlo, stla) = (-120, 48), so the region is extended
+    # westward to keep the waveform from sticking to the map boundary.
+    fig.basemap(region=[-135, -30, 35, 65], projection="M10c", frame=True)
+    fig.sac(
+        data=SAC_DATA,
+        amplitude_scale="1i",
+        # Use the reciprocal time scale, i.e., 0.5 cm per second, so that the
+        # 10-s waveform occupies 5 cm on the map.
+        time_scale="i0.5c",
+        pen="0.5p,red",
+    )
+    return fig
+
+
+@XFAIL_GMT_LE_6_6
+@pytest.mark.mpl_image_compare
+def test_sac_time_options():
+    """
+    Plot a SAC waveform with time alignment, shift, and reduction velocity.
+    """
+    fig = Figure()
+    fig.sac(
+        data=SAC_DATA,
+        reduction_velocity=8,
+        time_shift=2,
+        time_reference="o",
+        region=[5, 18, -2, 2],
+        projection="X15c/5c",
+        frame=True,
+        pen="0.5p,red",
+    )
+    return fig