timhoffm · Jul 27, 2026 · Sep 24, 2025 · Nov 30, 2025 · Dec 1, 2025 · Dec 1, 2025
diff --git a/doc/api/colors_api.rst b/doc/api/colors_api.rst
@@ -55,6 +55,7 @@ Multivariate Colormaps
    BivarColormap
    SegmentedBivarColormap
    BivarColormapFromImage
+   MultivarColormap
 
 Other classes
 -------------

diff --git a/lib/matplotlib/axes/_axes.py b/lib/matplotlib/axes/_axes.py
@@ -6238,6 +6238,7 @@ def imshow(self, X, cmap=None, norm=None, *, aspect=None,
             - (M, N): an image with scalar data. The values are mapped to
               colors using normalization and a colormap. See parameters *norm*,
               *cmap*, *vmin*, *vmax*.
+            - (K, M, N): if coupled with a cmap that supports K scalars
             - (M, N, 3): an image with RGB values (0-1 float or 0-255 int).
             - (M, N, 4): an image with RGBA values (0-1 float or 0-255 int),
               i.e. including transparency.
@@ -6247,15 +6248,16 @@ def imshow(self, X, cmap=None, norm=None, *, aspect=None,
 
             Out-of-range RGB(A) values are clipped.
 
-        %(cmap_doc)s
+
+        %(multi_cmap_doc)s
 
             This parameter is ignored if *X* is RGB(A).
 
-        %(norm_doc)s
+        %(multi_norm_doc)s
 
             This parameter is ignored if *X* is RGB(A).
 
-        %(vmin_vmax_doc)s
+        %(multi_vmin_vmax_doc)s
 
             This parameter is ignored if *X* is RGB(A).
 
@@ -6334,6 +6336,10 @@ def imshow(self, X, cmap=None, norm=None, *, aspect=None,
             See :doc:`/gallery/images_contours_and_fields/image_antialiasing` for
             a discussion of image antialiasing.
 
+            When using a `~matplotlib.colors.BivarColormap` or
+            `~matplotlib.colors.MultivarColormap`, 'data' is the only valid
+            interpolation_stage.
+
         alpha : float or array-like, optional
             The alpha blending value, between 0 (transparent) and 1 (opaque).
             If *alpha* is an array, the alpha blending values are applied pixel
@@ -6439,6 +6445,7 @@ def imshow(self, X, cmap=None, norm=None, *, aspect=None,
         if aspect is not None:
             self.set_aspect(aspect)
 
+        X = mcolorizer._ensure_multivariate_data(X, im.norm.n_components)
         im.set_data(X)
         im.set_alpha(alpha)
         if im.get_clip_path() is None:
@@ -6594,9 +6601,10 @@ def pcolor(self, *args, shading=None, alpha=None, norm=None, cmap=None,
 
         Parameters
         ----------
-        C : 2D array-like
+        C : 2D (M, N) or 3D (K, M, N) array-like
             The color-mapped values.  Color-mapping is controlled by *cmap*,
-            *norm*, *vmin*, and *vmax*.
+            *norm*, *vmin*, and *vmax*. 3D arrays are supported only if the
+            cmap supports K channels.
 
         X, Y : array-like, optional
             The coordinates of the corners of quadrilaterals of a pcolormesh::
@@ -6639,11 +6647,11 @@ def pcolor(self, *args, shading=None, alpha=None, norm=None, cmap=None,
             See :doc:`/gallery/images_contours_and_fields/pcolormesh_grids`
             for more description.
 
-        %(cmap_doc)s
+        %(multi_cmap_doc)s
 
-        %(norm_doc)s
+        %(multi_norm_doc)s
 
-        %(vmin_vmax_doc)s
+        %(multi_vmin_vmax_doc)s
 
         %(colorizer_doc)s
 
@@ -6718,8 +6726,19 @@ def pcolor(self, *args, shading=None, alpha=None, norm=None, cmap=None,
         if shading is None:
             shading = mpl.rcParams['pcolor.shading']
         shading = shading.lower()
-        X, Y, C, shading = self._pcolorargs('pcolor', *args, shading=shading,
-                                            kwargs=kwargs)
+
+        mcolorizer.ColorizingArtist._check_exclusionary_keywords(colorizer,
+                                                                 vmin=vmin, vmax=vmax,
+                                                                 norm=norm, cmap=cmap)
+        if colorizer is None:
+            colorizer = mcolorizer.Colorizer(cmap=cmap, norm=norm)
+
+        C = mcolorizer._ensure_multivariate_data(args[-1],
+                                                 colorizer.cmap.n_variates)
+
+        X, Y, C, shading = self._pcolorargs('pcolor', *args[:-1], C,
+                                            shading=shading, kwargs=kwargs)
+
         linewidths = (0.25,)
         if 'linewidth' in kwargs:
             kwargs['linewidths'] = kwargs.pop('linewidth')
@@ -6754,9 +6773,7 @@ def pcolor(self, *args, shading=None, alpha=None, norm=None, cmap=None,
         coords = stack([X, Y], axis=-1)
 
         collection = mcoll.PolyQuadMesh(
-            coords, array=C, cmap=cmap, norm=norm, colorizer=colorizer,
-            alpha=alpha, **kwargs)
-        collection._check_exclusionary_keywords(colorizer, vmin=vmin, vmax=vmax)
+            coords, array=C, colorizer=colorizer, alpha=alpha, **kwargs)
         collection._scale_norm(norm, vmin, vmax)
 
         coords = coords.reshape(-1, 2)  # flatten the grid structure; keep x, y
@@ -6794,6 +6811,7 @@ def pcolormesh(self, *args, alpha=None, norm=None, cmap=None, vmin=None,
             - (M, N) or M*N: a mesh with scalar data. The values are mapped to
               colors using normalization and a colormap. See parameters *norm*,
               *cmap*, *vmin*, *vmax*.
+            - (K, M, N): if coupled with a cmap that supports K scalars
             - (M, N, 3): an image with RGB values (0-1 float or 0-255 int).
             - (M, N, 4): an image with RGBA values (0-1 float or 0-255 int),
               i.e. including transparency.
@@ -6828,11 +6846,11 @@ def pcolormesh(self, *args, alpha=None, norm=None, cmap=None, vmin=None,
             expanded as needed into the appropriate 2D arrays, making a
             rectangular grid.
 
-        %(cmap_doc)s
+        %(multi_cmap_doc)s
 
-        %(norm_doc)s
+        %(multi_norm_doc)s
 
-        %(vmin_vmax_doc)s
+        %(multi_vmin_vmax_doc)s
 
         %(colorizer_doc)s
 
@@ -6956,16 +6974,24 @@ def pcolormesh(self, *args, alpha=None, norm=None, cmap=None, vmin=None,
         shading = mpl._val_or_rc(shading, 'pcolor.shading').lower()
         kwargs.setdefault('edgecolors', 'none')
 
-        X, Y, C, shading = self._pcolorargs('pcolormesh', *args,
+        mcolorizer.ColorizingArtist._check_exclusionary_keywords(colorizer,
+                                                                 vmin=vmin, vmax=vmax,
+                                                                 norm=norm, cmap=cmap)
+        if colorizer is None:
+            colorizer = mcolorizer.Colorizer(cmap=cmap, norm=norm)
+
+        C = mcolorizer._ensure_multivariate_data(args[-1],
+                                                 colorizer.cmap.n_variates)
+
+        X, Y, C, shading = self._pcolorargs('pcolormesh', *args[:-1], C,
                                             shading=shading, kwargs=kwargs)
         coords = np.stack([X, Y], axis=-1)
 
         kwargs.setdefault('snap', mpl.rcParams['pcolormesh.snap'])
 
         collection = mcoll.QuadMesh(
             coords, antialiased=antialiased, shading=shading,
-            array=C, cmap=cmap, norm=norm, colorizer=colorizer, alpha=alpha, **kwargs)
-        collection._check_exclusionary_keywords(colorizer, vmin=vmin, vmax=vmax)
+            array=C, colorizer=colorizer, alpha=alpha, **kwargs)
         collection._scale_norm(norm, vmin, vmax)
 
         coords = coords.reshape(-1, 2)  # flatten the grid structure; keep x, y
@@ -8911,6 +8937,9 @@ def matshow(self, Z, **kwargs):
 
         """
         Z = np.asanyarray(Z)
+        if Z.ndim != 2:
+            if Z.ndim != 3 or Z.shape[2] not in (1, 3, 4):
+                raise TypeError(f"Invalid shape {Z.shape} for image data")
         kw = {'origin': 'upper',
               'interpolation': 'nearest',
               'aspect': 'equal',          # (already the imshow default)

diff --git a/lib/matplotlib/axes/_axes.pyi b/lib/matplotlib/axes/_axes.pyi
@@ -12,7 +12,13 @@ from matplotlib.collections import (
     QuadMesh,
 )
 from matplotlib.colorizer import Colorizer
-from matplotlib.colors import Colormap, Normalize
+from matplotlib.colors import (
+    Colormap,
+    BivarColormap,
+    MultivarColormap,
+    Norm,
+    Normalize,
+)
 from matplotlib.container import (
     BarContainer, PieContainer, ErrorbarContainer, StemContainer)
 from matplotlib.contour import ContourSet, QuadContourSet
@@ -503,14 +509,14 @@ class Axes(_AxesBase):
     def imshow(
         self,
         X: ArrayLike | PIL.Image.Image,
-        cmap: str | Colormap | None = ...,
-        norm: str | Normalize | None = ...,
+        cmap: str | Colormap | BivarColormap | MultivarColormap | None = ...,
+        norm: str | Norm | None = ...,
         *,
         aspect: Literal["equal", "auto"] | float | None = ...,
         interpolation: str | None = ...,
         alpha: float | ArrayLike | None = ...,
-        vmin: float | None = ...,
-        vmax: float | None = ...,
+        vmin: float | tuple[float, ...] | None = ...,
+        vmax: float | tuple[float, ...] | None = ...,
         colorizer: Colorizer | None = ...,
         origin: Literal["upper", "lower"] | None = ...,
         extent: tuple[float, float, float, float] | None = ...,
@@ -527,10 +533,10 @@ class Axes(_AxesBase):
         *args: ArrayLike,
         shading: Literal["flat", "nearest", "auto"] | None = ...,
         alpha: float | None = ...,
-        norm: str | Normalize | None = ...,
-        cmap: str | Colormap | None = ...,
-        vmin: float | None = ...,
-        vmax: float | None = ...,
+        norm: str | Norm | None = ...,
+        cmap: str | Colormap | BivarColormap | MultivarColormap | None = ...,
+        vmin: float | tuple[float, ...] | None = ...,
+        vmax: float | tuple[float, ...] | None = ...,
         colorizer: Colorizer | None = ...,
         data=...,
         **kwargs
@@ -539,10 +545,10 @@ class Axes(_AxesBase):
         self,
         *args: ArrayLike,
         alpha: float | None = ...,
-        norm: str | Normalize | None = ...,
-        cmap: str | Colormap | None = ...,
-        vmin: float | None = ...,
-        vmax: float | None = ...,
+        norm: str | Norm | None = ...,
+        cmap: str | Colormap | BivarColormap | MultivarColormap | None = ...,
+        vmin: float | tuple[float, ...] | None = ...,
+        vmax: float | tuple[float, ...] | None = ...,
         colorizer: Colorizer | None = ...,
         shading: Literal["flat", "nearest", "gouraud", "auto"] | None = ...,
         antialiased: bool = ...,

diff --git a/lib/matplotlib/collections.py b/lib/matplotlib/collections.py
@@ -2386,6 +2386,8 @@ def set_array(self, A):
             h, w = height, width
         ok_shapes = [(h, w, 3), (h, w, 4), (h, w), (h * w,)]
         if A is not None:
+            if hasattr(self, 'norm'):
+                A = mcolorizer._ensure_multivariate_data(A, self.norm.n_components)
             shape = np.shape(A)
             if shape not in ok_shapes:
                 raise ValueError(
@@ -2643,7 +2645,7 @@ def _get_unmasked_polys(self):
         mask = (mask[0:-1, 0:-1] | mask[1:, 1:] | mask[0:-1, 1:] | mask[1:, 0:-1])
         arr = self.get_array()
         if arr is not None:
-            arr = np.ma.getmaskarray(arr)
+            arr = self._getmaskarray(arr)
             if arr.ndim == 3:
                 # RGB(A) case
                 mask |= np.any(arr, axis=-1)

diff --git a/lib/matplotlib/collections.pyi b/lib/matplotlib/collections.pyi
@@ -7,7 +7,12 @@ from numpy.typing import ArrayLike, NDArray
 from . import colorizer, transforms
 from .backend_bases import MouseEvent
 from .artist import Artist
-from .colors import Normalize, Colormap
+from .colors import (
+    Colormap,
+    BivarColormap,
+    MultivarColormap,
+    Norm,
+)
 from .lines import Line2D
 from .path import Path
 from .patches import Patch
@@ -29,8 +34,8 @@ class Collection(colorizer.ColorizingArtist):
         antialiaseds: bool | Sequence[bool] | None = ...,
         offsets: tuple[float, float] | Sequence[tuple[float, float]] | None = ...,
         offset_transform: transforms.Transform | None = ...,
-        norm: Normalize | None = ...,
-        cmap: Colormap | None = ...,
+        norm: Norm | None = ...,
+        cmap: Colormap | BivarColormap | MultivarColormap | None = ...,
         colorizer: colorizer.Colorizer | None = ...,
         pickradius: float = ...,
         hatch: str | None = ...,

diff --git a/lib/matplotlib/colorizer.py b/lib/matplotlib/colorizer.py
@@ -74,7 +74,7 @@ def _scale_norm(self, norm, vmin, vmax, A):
         """
         if vmin is not None or vmax is not None:
             self.set_clim(vmin, vmax)
-            if isinstance(norm, colors.Normalize):
+            if isinstance(norm, colors.Norm):
                 raise ValueError(
                     "Passing a Normalize instance simultaneously with "
                     "vmin/vmax is not supported.  Please pass vmin/vmax "
@@ -277,8 +277,16 @@ def set_clim(self, vmin=None, vmax=None):
     def get_clim(self):
         """
         Return the values (min, max) that are mapped to the colormap limits.
+
+        This function always returns min and max as tuples to ensure type consistency
+        when working with both scalar and multivariate color mapping.
+        See also `._ColorizerInterface.get_clim()` which returns scalars but is
+        unavailable for multivariate color mapping.
         """
-        return self.norm.vmin, self.norm.vmax
+        if self.norm.n_components == 1:
+            return (self.norm.vmin, ), (self.norm.vmax, )
+        else:
+            return self.norm.vmin, self.norm.vmax
 
     def changed(self):
         """
@@ -306,7 +314,10 @@ def vmax(self, vmax):
 
     @property
     def clip(self):
-        return self.norm.clip
+        if self.norm.n_components == 1:
+            return (self.norm.clip, )
+        else:
+            return self.norm.clip
 
     @clip.setter
     def clip(self, clip):
@@ -360,8 +371,14 @@ def to_rgba(self, x, alpha=None, bytes=False, norm=True):
     def get_clim(self):
         """
         Return the values (min, max) that are mapped to the colormap limits.
+
+        This function is not available for multivariate data.
+        Use `.Colorizer.get_clim` via the .colorizer property instead.
         """
-        return self._colorizer.get_clim()
+        if self._colorizer.norm.n_components > 1:
+            raise AttributeError("get_clim() cannot be used with a multi-component "
+                                 "colormap. Use .colorizer.get_clim() instead")
+        return self.colorizer.norm.vmin, self.colorizer.norm.vmax
 
     def set_clim(self, vmin=None, vmax=None):
         """
@@ -376,9 +393,14 @@ def set_clim(self, vmin=None, vmax=None):
              tuple (*vmin*, *vmax*) as a single positional argument.
 
              .. ACCEPTS: (vmin: float, vmax: float)
+
+        This function is not available for multivariate data.
         """
         # If the norm's limits are updated self.changed() will be called
         # through the callbacks attached to the norm
+        if self._colorizer.norm.n_components > 1:
+            raise AttributeError("set_clim() cannot be used with a multi-component "
+                                 "colormap. Use .colorizer.set_clim() instead")
         self._colorizer.set_clim(vmin, vmax)
 
     def get_alpha(self):
@@ -600,6 +622,18 @@ def get_array(self):
         """
         return self._A
 
+    def _getmaskarray(self, A):
+        """
+        Similar to np.ma.getmaskarray but also handles the case where
+        the data has multiple fields.
+
+        The return array always has the same shape as the input, and dtype bool
+        """
+        mask = np.ma.getmaskarray(A)
+        if isinstance(self.norm, colors.MultiNorm):
+            mask = np.any(mask.view('bool').reshape((*A.shape, -1)), axis=-1)
+        return mask
+
     def changed(self):
         """
         Call this whenever the mappable is changed to notify all the