Proposal: Support explicit mesh coordinates in `xarray.plot.pcolormesh`

Hi @dcherian @jhamman, @kmuehlbauer

I would like to open a PR to address a current limitation in xarray.plot.pcolormesh regarding the correct display of:

• Irregular 1D grids (non-uniform bin spacing)
• Curvilinear 2D grids (e.g., distorted or geophysical meshes)

Other linked issues:

• https://github.com/pydata/xarray/issues/792
• https://github.com/pydata/xarray/issues/4364
• https://github.com/pydata/xarray/issues/5901

---

Problem

xarray.plot.pcolormesh assumes that coordinates represent cell centers and internally reconstructs edges (infer_intervals=True logic).

This leads to incorrect rendering when:

• Grid spacing is irregular
• True cell boundaries are known and should be respected
• The grid is curvilinear and defined by 2D vertex coordinates

Such grids are common in many domains:

• Numerical Weather Prediction (NWP) model grids
• Satellite swath data
• Radar products
• Ocean and climate model outputs
• ...

This limitation significantly impacts accurate scientific and cartographic visualization.

---

Current Workaround

A workaround exists: directly calling matplotlib.pyplot.pcolormesh with explicit edge arrays.

However, this comes with significant drawbacks:

• Loss of xarray’s high-level plotting API
• Loss of automatic labeling and metadata handling
• Loss of FacetGrid support (e.g., faceting over time or ensemble dimensions)
• More verbose code

In practice, this makes advanced visualization workflows much harder, especially when working with multi-dimensional datasets.

Describe the solution you'd like

I propose adding support for explicitly passing edge coordinates to pcolormesh, following the matplotlib.pyplot.pcolormesh(..., shading="flat") convention.

1D Case

• x_edges: shape (n+1,)
• y_edges: shape (m+1,)

2D Curvilinear Case

• x_edges: shape (m+1, n+1)
• y_edges: shape (m+1, n+1)

When provided, these edge coordinates would:

• Bypass edge inference
• Be passed directly to Matplotlib
• Ensure geometrically correct rendering

I suggest therefore to introduce x_edges and y_edges keyword arguments to pcolormesh

da.pcolormesh(..., x_edges=..., y_edges=...)

Under this approach:

• x and y must still be specified
• If x_edges / y_edges are provided, edge inference is skipped
• The edges /mesh arrays are passed directly to Matplotlib

Before starting work on the PR, I would appreciate your feedback on:

• Whether you are in favor of adding this capability
• Preferred API design direction
• Any architectural considerations to keep in mind

---

Minimal, reproducible example of the problem

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt

####-----------------------------------------------------------------
#### Rectilinear grid (1D case)
# Regular edges
x_edges = np.arange(0, 4)      # 4 edges
y_edges = np.arange(6, 9)      # 3 edges

# Irregular edges  
x_edges = np.array([0, 1, 4, 10]) # 4 edges
y_edges = np.array([30, 50, 100]) # 3 edges

# Plot limits 
xlim = x_edges.min(), x_edges.max()
ylim = y_edges.min(), y_edges.max()

# Compute centers
x_centers = 0.5 * (x_edges[:-1] + x_edges[1:])
y_centers = 0.5 * (y_edges[:-1] + y_edges[1:])

# Create 2D data
Xc, Yc = np.meshgrid(x_centers, y_centers)
data = np.sin(Xc/2) * np.cos(Yc/5)

# Build Dataset
ds = xr.Dataset(
    {
        "data": (("y", "x"), data)
    },
    coords={
        "x": x_centers,
        "y": y_centers,
    }
)

####-----------------------------------------------------------------
#### Curvilinear grid (2D case)
# Create 2D irregular edges from your 1D irregular edge
x_edges, y_edges = np.meshgrid(x_edges, y_edges)

# Add distortion (to make geometry clearly non-rectangular)
x_edges = x_edges + 0.8 * np.sin(y_edges / 40)
y_edges = y_edges + 3.0 * np.cos(x_edges / 5)

# Compute cell centers from corner average
Xc2 = 0.25 * (
    x_edges[:-1, :-1] +
    x_edges[1:, :-1] +
    x_edges[:-1, 1:] +
    x_edges[1:, 1:]
)

Yc2 = 0.25 * (
    y_edges[:-1, :-1] +
    y_edges[1:, :-1] +
    y_edges[:-1, 1:] +
    y_edges[1:, 1:]
)

data2 = np.sin(Xc2/2) * np.cos(Yc2/5)

ds = xr.Dataset(
    {
        "data": (("y", "x"), data2)
    },
    coords={
        "x": (("y", "x"), Xc2),
        "y": (("y", "x"), Yc2),
    }
)

# Plot limits from true edges
xlim = x_edges.min(), x_edges.max()
ylim = y_edges.min(), y_edges.max()

#--------------------------------------------------------------------
#### XARRAY PCOLORMESH PLOTTING IS CURRENTLY WRONG 
p = ds["data"].plot.pcolormesh(x="x", y="y", add_colorbar=False)
p.axes.set_title("")
p.axes.set_xlabel("")
p.axes.set_ylabel("")
p.axes.set_xlim(*xlim)
p.axes.set_ylim(*ylim)

#--------------------------------------------------------------------
#### PLOT CORRECTLY WITH MATPLOTLIB PCOLORMESH 
fig, ax = plt.subplots(1,1)
ax.pcolormesh(x_edges, y_edges, data)
p.axes.set_title("")
p.axes.set_xlabel("")
p.axes.set_ylabel("")
p.axes.set_xlim(*xlim)
p.axes.set_ylim(*ylim)

Start at the xarray.plot.pcolormesh entry point and review how its infer_intervals logic handles coordinate arrays. Define the API and tests for 1D and 2D x_edges/y_edges, ensuring explicit edges bypass inference and are passed to Matplotlib while existing labeling and FacetGrid behavior remains intact.