diff --git a/CHANGELOG.md b/CHANGELOG.md
index 8138f07..ba90f0e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,11 @@
 # Changelog
 
+## [unreleased]
+
+### Added
+
+- `gridplot` with `PyPlot` has now a color bar with annotation `boundary regions`
+
 ## [1.9.0] - 2024-11-25
 
 ### Added
diff --git a/src/pyplot.jl b/src/pyplot.jl
index 4c39c4d..8f63837 100644
--- a/src/pyplot.jl
+++ b/src/pyplot.jl
@@ -236,30 +236,45 @@ function gridplot!(ctx, TP::Type{PyPlotType}, ::Type{Val{2}}, grid)
     brflag = ones(Bool, nbfaceregions)
     ax.set_aspect(ctx[:aspect])
     tridat = tridata(grid, ctx[:gridscale])
-    cmap = region_cmap(max(ncellregions, 5))
-    cdata = ax.tripcolor(
+    # PyPlot.ColorMap cannot handle n ≤ 1, TODO: use single color instead of a color map
+    cmap = region_cmap(max(2, ncellregions))
+    bcmap = bregion_cmap(max(2, nbfaceregions))
+
+    cell_colors = cellcolors(grid, ctx[:cellcoloring])
+
+    # dummy plot to get a correct color bar for the boundary data
+    bcdata = ax.tripcolor(
+        tridat[1][1:3], tridat[2][1:3]; # extract a single point from the original triangulation
+        facecolors = cell_colors[1:1],  # only one triangle!
+        cmap = PyPlot.ColorMap(bcmap, length(bcmap)),
+        vmin = 0.5,
+        vmax = length(bcmap) + 0.5,
+    )
+
+    ax.tripcolor(
         tridat...;
-        facecolors = cellcolors(grid, ctx[:cellcoloring]),
+        facecolors = cell_colors,
         cmap = PyPlot.ColorMap(cmap, length(cmap)),
-        vmin = 1.0,
-        vmax = length(cmap),
     )
+
     if ctx[:show_colorbar]
         if ctx[:colorbar] == :horizontal
             cbar = fig.colorbar(
-                cdata;
+                bcdata;
                 ax = ax,
-                ticks = collect(1:length(cmap)),
+                ticks = collect(1:length(bcmap)),
                 orientation = "horizontal",
+                label = "boundary regions",
             )
         end
 
         if ctx[:colorbar] == :vertical
             cbar = fig.colorbar(
-                cdata;
+                bcdata;
                 ax = ax,
-                ticks = collect(1:length(cmap)),
+                ticks = collect(1:length(bcmap)),
                 orientation = "vertical",
+                label = "boundary regions",
             )
         end
     end