Maximum value not always included in pcolormesh output when levels are requested

What happened?

When plotting data with "pcolormesh" and requesting levels, if the calculated "vmax" if equal to the upper bound of the last bin, it is not included in the map, as if it was "over" or "missing".

What did you expect to happen?

I expected values == vmax to be painted with the color of the last bin, which is what contourf does.

Minimal Complete Verifiable Example

# /// script
# requires-python = ">=3.11"
# dependencies = [
#   "xarray[complete]@git+https://github.com/pydata/xarray.git@main",
# ]
# ///
#
# This script automatically imports the development branch of xarray to check for issues.
# Please delete this header if you have _not_ tested this script with `uv run`!

import xarray as xr
xr.show_versions()
# your reproducer code ...

import xarray as xr
ds = xr.tutorial.open_dataset('air_temperature')
air = ds.air.isel(time=0).round(-1)
air.plot(levels=5)

Steps to reproduce

The example above gives the image below:

You can see that with the heavy discretization I imposed, the values at the bottom are equal to the calculated "vmax" and are not shown.

MVCE confirmation

• Minimal example — the example is as focused as reasonably possible to demonstrate the underlying issue in xarray.
• Complete example — the example is self-contained, including all data and the text of any traceback.
• Verifiable example — the example copy & pastes into an IPython prompt or Binder notebook, returning the result.
• New issue — a search of GitHub Issues suggests this is not a duplicate.
• Recent environment — the issue occurs with the latest version of xarray and its dependencies.

Relevant log output

Anything else we need to know?

Matplotlib's from_levels_and_colors, it is said that :

The quantization levels used to construct the BoundaryNorm. Value v is quantized to level i if lev[i] <= v < lev[i+1].

So this explains why xarray's implementation gives no color. However, the result seems to be unexpected and incoherent with the output of contourf. I think it is logic to expect the maximum value to be included in the plot when simply passing levels as a scalar.

Environment

Start by running the provided xarray pcolormesh example with scalar levels and inspect the plotting path that constructs the Matplotlib normalization. Compare the vmax boundary behavior with contourf and Matplotlib's from_levels_and_colors documentation. Done means values equal to the calculated vmax receive the final bin color, with a regression test covering the example.