diff --git a/examples/gallery/images/gridding_with_surface.py b/examples/gallery/images/gridding_with_surface.py
new file mode 100644
index 00000000000..93c5916e094
--- /dev/null
+++ b/examples/gallery/images/gridding_with_surface.py
@@ -0,0 +1,85 @@
+"""
+Gridding scattered data
+=======================
+
+The :func:`pygmt.surface` function interpolates irregularly spaced (x, y, z)
+observations onto a regular grid. Before interpolation, it's usually recommended
+to preprocess the data using one of :func:`pygmt.blockmean`,
+:func:`pygmt.blockmedian`, and :func:`pygmt.blockmode`, to avoid aliasing short
+wavelengths, by returning one value for every occupied block.
+
+"""
+
+# %%
+import pygmt
+from pygmt.params import Axis, Frame, Position
+
+# Load irregular ship-track observations of bathymetry off Baja California.
+data = pygmt.datasets.load_sample_data(name="bathymetry")
+
+region = [245, 255, 20, 30]
+spacing = "5m"  # 5 arc-minutes
+
+# Calculate one representative median value for each occupied block.
+reduced_data = pygmt.blockmedian(
+    data=data, region=region, spacing=spacing, center=True
+)
+
+# Interpolate the reduced observations onto a regular grid. The value "d"
+# constrains the grid to the minimum and maximum input data values.
+grid = pygmt.surface(
+    data=reduced_data, region=region, spacing=spacing, lower="d", upper="d"
+)
+
+fig = pygmt.Figure()
+pygmt.makecpt(cmap="gmt/geo", series=[-7500, 0, 500])
+pygmt.config(FONT_TITLE="12p", MAP_TITLE_OFFSET="4p")
+
+with fig.subplot(nrows=1, ncols=2, figsize=("20c", "9c")):
+    # Block-median observations.
+    with fig.set_panel(panel=0):
+        fig.basemap(
+            region=region,
+            projection="M?",
+            frame=Frame(
+                axes="WSne",
+                axis=Axis(annot=2, tick=1),
+                title="Block-median observations",
+            ),
+        )
+        fig.plot(
+            x=reduced_data.longitude,
+            y=reduced_data.latitude,
+            fill=reduced_data.bathymetry,
+            style="c0.05c",
+            cmap=True,
+        )
+        fig.coast(land="gray", shorelines="0.5p,black")
+
+    # Interpolated surface grid.
+    with fig.set_panel(panel=1):
+        fig.grdimage(
+            grid=grid,
+            region=region,
+            projection="M?",
+            frame=Frame(
+                axes="WSne",
+                axis=Axis(annot=2, tick=1),
+                title="Interpolated surface",
+            ),
+            cmap=True,
+        )
+        fig.coast(land="gray", shorelines="0.5p,black")
+
+# Position the color bar using explicit plot coordinates:
+# 10 cm is the horizontal center of the 20 cm wide subplot.
+fig.colorbar(
+    position=Position(("10c", "-1.2c"), cstype="plotcoords", anchor="TC"),
+    length=12,
+    width=0.4,
+    orientation="horizontal",
+    annot=1000,
+    label="Bathymetry (m)",
+)
+
+fig.show()