diff --git a/examples/gallery/histograms/blockm.py b/examples/gallery/histograms/blockm.py
new file mode 100644
index 00000000000..87c35347392
--- /dev/null
+++ b/examples/gallery/histograms/blockm.py
@@ -0,0 +1,63 @@
+"""
+Blockmean
+---------
+The :meth:`pygmt.blockmean` method calculates different quantities
+inside blocks/bins whose dimensions are defined via the ``spacing`` parameter.
+The following examples show how to calculate the averages of the given values
+inside each bin and how to report the number of points inside each bin.
+"""
+
+import pygmt
+
+# Load sample data
+data = pygmt.datasets.load_japan_quakes()
+# Select only needed columns
+data = data[["longitude", "latitude", "depth_km"]]
+
+# Set the region for the plot
+region = [130, 152.5, 32.5, 52.5]
+# Define spacing in x and y direction (150 by 150 minute blocks)
+spacing = "150m"
+
+fig = pygmt.Figure()
+
+# Calculate mean depth in km from all events within 150x150 minute
+# bins using blockmean
+df = pygmt.blockmean(data=data, region=region, spacing=spacing)
+# convert to grid
+grd = pygmt.xyz2grd(data=df, region=region, spacing=spacing)
+
+fig.grdimage(
+    grid=grd,
+    region=region,
+    frame=["af", '+t"Mean earthquake depth inside each block"'],
+    cmap="batlow",
+)
+# plot slightly transparent landmasses on top
+fig.coast(land="darkgray", transparency=40)
+# plot original data points
+fig.plot(
+    x=data.longitude, y=data.latitude, style="c0.3c", color="white", pen="1p,black"
+)
+fig.colorbar(frame=["x+lkm"])
+
+fig.shift_origin(xshift="w+5c")
+
+# Calculate number of total locations within 150x150 minute bins via
+# blockmean's summary parameter
+df = pygmt.blockmean(data=data, region=region, spacing=spacing, summary="n")
+grd = pygmt.xyz2grd(data=df, region=region, spacing=spacing)
+
+fig.grdimage(
+    grid=grd,
+    region=region,
+    frame=["af", '+t"Number of points inside each block"'],
+    cmap="batlow",
+)
+fig.coast(land="darkgray", transparency=40)
+fig.plot(
+    x=data.longitude, y=data.latitude, style="c0.3c", color="white", pen="1p,black"
+)
+fig.colorbar(frame=["x+lcount"])
+
+fig.show()