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

[ghiggi]

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:

• ENH: Don't infer pcolormesh interval breaks for unevenly spaced coordiantes #792
• plt.pcolormesh will infer interval breaks per default #4364
• Spurious lines of the pcolormesh example #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)

[headtr1ck]

The major problem I see is that the N+1 dimensions must have different names from the N dimensions so an automatic mapping is not possible unless one assumes some kind of naming convention which we try to avoid.

[ghiggi]

Hi @headtr1ck !

Thanks for raising this — I may have misunderstood your concern, so let me try to clarify what I had in mind.

I agree that the mesh is unlikely to be contained inside the Dataset object, since its dimensions are necessarily (N+1, M+1) and would therefore require dedicated dimension names distinct from the other data x and y dimensions.

I am not proposing that x_edges and y_edges be provided as strings referencing mesh dimensions within the Dataset.

Instead, what I’m proposing is that the user explicitly provides x_edges and y_edges in one of two forms:

• As NumPy arrays with shape (N+1, M+1), consistent with the declared x (N) and y(M) dimensions of the data.

• Preferably, as a DataArray, with the same dimension names as x and y, but with shape (N+1, M+1).

The second option is safer because the explicit dimension names allow us to align or transpose the provided edge arrays to match the declared x and y ordering if necessary, without relying just on the array shape.

In both cases, the resulting NumPy arrays (either provided directly or extracted from the DataArray) would be passed to matplotlib.pcolormesh, after validating that their shape is indeed (N+1, M+1).