diff --git a/examples/images_contours_and_fields/bivariate_demo.py b/examples/images_contours_and_fields/bivariate_demo.py
@@ -0,0 +1,44 @@
+"""
+==============
+Bivariate Demo
+==============
+
+Plotting bivariate data.
+
+imshow, pcolor, pcolormesh, pcolorfast allows you to plot bivariate data
+using a bivaraite colormap.
+
+In this example we use imshow to plot air temperature with surface pressure
+alongwith a color square.
+"""
+import matplotlib.colors as colors
+from matplotlib.cbook import get_sample_data
+import matplotlib.pyplot as plt
+import numpy as np
+
+
+###############################################################################
+# Bivariate plotting demo
+# -----------------------
+
+air_temp = np.load(get_sample_data('air_temperature.npy'))
+surf_pres = np.load(get_sample_data('surface_pressure.npy'))
+
+fig, ax = plt.subplots()
+
+bivariate = [air_temp, surf_pres]
+
+###############################################################################
+# To distinguish bivariate data either BivariateNorm or BivariateColormap must
+# be passed in as argument
+
+cax = ax.imshow(bivariate, norm=colors.BivariateNorm(),
+                cmap=colors.BivariateColormap())
+
+###############################################################################
+# If input data is bivariate then colorbar automatically draws colorsquare
+# instead of colorbar
+
+cbar = fig.colorbar(cax, xlabel='air_temp', ylabel='surf_pres')
+
+plt.show()
diff --git a/lib/matplotlib/axes/_axes.py b/lib/matplotlib/axes/_axes.py
@@ -4044,7 +4044,8 @@ def scatter(self, x, y, s=None, c=None, marker=None, cmap=None, norm=None,
 
         if colors is None:
             if norm is not None and not isinstance(norm, mcolors.Normalize):
-                msg = "'norm' must be an instance of 'mcolors.Normalize'"
+                msg = ("'norm' must be an instance of 'mcolors.Normalize' or "
+                       "'mcolors.BivariateNorm'")
                 raise ValueError(msg)
             collection.set_array(np.asarray(c))
             collection.set_cmap(cmap)
@@ -4403,7 +4404,8 @@ def hexbin(self, x, y, C=None, gridsize=100, bins=None,
             accum = bins.searchsorted(accum)
 
         if norm is not None and not isinstance(norm, mcolors.Normalize):
-            msg = "'norm' must be an instance of 'mcolors.Normalize'"
+            msg = ("'norm' must be an instance of 'mcolors.Normalize' or "
+                   "'mcolors.BivariateNorm'")
             raise ValueError(msg)
         collection.set_array(accum)
         collection.set_cmap(cmap)
@@ -5037,23 +5039,25 @@ def imshow(self, X, cmap=None, norm=None, aspect=None,
 
         Parameters
         ----------
-        X : array_like, shape (n, m) or (n, m, 3) or (n, m, 4)
+        X : array_like, shape (n, m) or (n, m, 3) or (n, m, 4) or (2, n, m)
             Display the image in `X` to current axes.  `X` may be an
             array or a PIL image. If `X` is an array, it
             can have the following shapes and types:
 
-            - MxN -- values to be mapped (float or int)
+            - MxN -- univariate values to be mapped (float or int)
             - MxNx3 -- RGB (float or uint8)
             - MxNx4 -- RGBA (float or uint8)
+            - 2xMxN -- bivariate values to be mapped (float or int)
 
             The value for each component of MxNx3 and MxNx4 float arrays
-            should be in the range 0.0 to 1.0. MxN arrays are mapped
+            should be in the range 0.0 to 1.0. MxN and 2xMxN arrays are mapped
             to colors based on the `norm` (mapping scalar to scalar)
             and the `cmap` (mapping the normed scalar to a color).
 
-        cmap : `~matplotlib.colors.Colormap`, optional, default: None
+        cmap : `~matplotlib.colors.Colormap`, \
+            `~matplotlib.colors.BivariateColormap`, optional, default: None
             If None, default to rc `image.cmap` value. `cmap` is ignored
-            if `X` is 3-D, directly specifying RGB(A) values.
+            if `X` is 3-D but not bivariate, directly specifying RGB(A) values.
 
         aspect : ['auto' | 'equal' | scalar], optional, default: None
             If 'auto', changes the image aspect ratio to match that of the
@@ -5077,7 +5081,8 @@ def imshow(self, X, cmap=None, norm=None, aspect=None,
             on the Agg, ps and pdf backends. Other backends will fall back to
             'nearest'.
 
-        norm : `~matplotlib.colors.Normalize`, optional, default: None
+        norm : `~matplotlib.colors.Normalize`, \
+            `matplotlib.colors.BivariateNorm`, optional, default: None
             A `~matplotlib.colors.Normalize` instance is used to scale
             a 2-D float `X` input to the (0, 1) range for input to the
             `cmap`. If `norm` is None, use the default func:`normalize`.
@@ -5137,16 +5142,29 @@ def imshow(self, X, cmap=None, norm=None, aspect=None,
         of pixel (0, 0).
 
         """
-
         if not self._hold:
             self.cla()
 
-        if norm is not None and not isinstance(norm, mcolors.Normalize):
-            msg = "'norm' must be an instance of 'mcolors.Normalize'"
+        if norm is not None and not isinstance(norm, mcolors.Norms):
+            msg = ("'norm' must be an instance of 'mcolors.Normalize' or "
+                   "'mcolors.BivariateNorm'")
             raise ValueError(msg)
+
+        temp = np.asarray(X)
+        is_bivari = (isinstance(norm, mcolors.BivariateNorm) or
+                     isinstance(cmap, mcolors.BivariateColormap))
+        if is_bivari:
+            if temp.ndim != 3 and temp.shape[0] != 2:
+                raise TypeError("Expected shape like (2, n, m)")
+            if cmap is None:
+                cmap = mcolors.BivariateColormap()
+            if norm is None:
+                norm = mcolors.BivariateNorm()
+
         if aspect is None:
             aspect = rcParams['image.aspect']
         self.set_aspect(aspect)
+
         im = mimage.AxesImage(self, cmap, norm, interpolation, origin, extent,
                               filternorm=filternorm, filterrad=filterrad,
                               resample=resample, **kwargs)
@@ -5173,7 +5191,6 @@ def imshow(self, X, cmap=None, norm=None, aspect=None,
 
     @staticmethod
     def _pcolorargs(funcname, *args, **kw):
-        # This takes one kwarg, allmatch.
         # If allmatch is True, then the incoming X, Y, C must
         # have matching dimensions, taking into account that
         # X and Y can be 1-D rather than 2-D.  This perfect
@@ -5186,10 +5203,17 @@ def _pcolorargs(funcname, *args, **kw):
         # is False.
 
         allmatch = kw.pop("allmatch", False)
+        norm = kw.pop("norm", None)
+        cmap = kw.pop("cmap", None)
 
         if len(args) == 1:
             C = np.asanyarray(args[0])
-            numRows, numCols = C.shape
+            is_bivari = (isinstance(norm, mcolors.BivariateNorm) or
+                         isinstance(cmap, mcolors.BivariateColormap))
+            if is_bivari:
+                numRows, numCols = C.shape[1:]
+            else:
+                numRows, numCols = C.shape
             if allmatch:
                 X, Y = np.meshgrid(np.arange(numCols), np.arange(numRows))
             else:
@@ -5200,7 +5224,12 @@ def _pcolorargs(funcname, *args, **kw):
 
         if len(args) == 3:
             X, Y, C = [np.asanyarray(a) for a in args]
-            numRows, numCols = C.shape
+            is_bivari = (isinstance(norm, mcolors.BivariateNorm) or
+                         isinstance(cmap, mcolors.BivariateColormap))
+            if is_bivari:
+                numRows, numCols = C.shape[1:]
+            else:
+                numRows, numCols = C.shape
         else:
             raise TypeError(
                 'Illegal arguments to %s; see help(%s)' % (funcname, funcname))
@@ -5235,7 +5264,7 @@ def _pcolorargs(funcname, *args, **kw):
     @docstring.dedent_interpd
     def pcolor(self, *args, **kwargs):
         """
-        Create a pseudocolor plot of a 2-D array.
+        Create a pseudocolor plot of a 2-D univariate or 3-D bivariate array.
 
         Call signatures::
 
@@ -5273,10 +5302,12 @@ def pcolor(self, *args, **kwargs):
             vectors, they will be expanded as needed into the appropriate 2-D
             arrays, making a rectangular grid.
 
-        cmap : `~matplotlib.colors.Colormap`, optional, default: None
+        cmap : `~matplotlib.colors.Colormap` or \
+            `matplotlib.colors.BivariateColormap`, optional, default: None
             If `None`, default to rc settings.
 
-        norm : `matplotlib.colors.Normalize`, optional, default: None
+        norm : `matplotlib.colors.Normalize` or \
+            `matplotlib.colors.BivariateNorm`, optional, default: None
             An instance is used to scale luminance data to (0, 1).
             If `None`, defaults to :func:`normalize`.
 
@@ -5382,9 +5413,20 @@ def pcolor(self, *args, **kwargs):
         vmin = kwargs.pop('vmin', None)
         vmax = kwargs.pop('vmax', None)
 
-        X, Y, C = self._pcolorargs('pcolor', *args, allmatch=False)
+        kw = {'norm': norm, 'cmap': cmap, 'allmatch': False}
+        X, Y, C = self._pcolorargs('pcolor', *args, **kw)
         Ny, Nx = X.shape
 
+        is_bivari = (isinstance(norm, mcolors.BivariateNorm) or
+                     isinstance(cmap, mcolors.BivariateColormap))
+        if is_bivari:
+            if C.ndim != 3 and C.shape[0] != 2:
+                raise TypeError("Expected shape like (2, n, m)")
+            if cmap is None:
+                cmap = mcolors.BivariateColormap()
+            if norm is None:
+                norm = mcolors.BivariateNorm()
+
         # unit conversion allows e.g. datetime objects as axis values
         self._process_unit_info(xdata=X, ydata=Y, kwargs=kwargs)
         X = self.convert_xunits(X)
@@ -5399,7 +5441,10 @@ def pcolor(self, *args, **kwargs):
         xymask = (mask[0:-1, 0:-1] + mask[1:, 1:] +
                   mask[0:-1, 1:] + mask[1:, 0:-1])
         # don't plot if C or any of the surrounding vertices are masked.
-        mask = ma.getmaskarray(C) + xymask
+        if isinstance(norm, mcolors.BivariateNorm):
+            mask = ma.getmaskarray(C[0]) + ma.getmaskarray(C[1]) + xymask
+        else:
+            mask = ma.getmaskarray(C) + xymask
 
         newaxis = np.newaxis
         compress = np.compress
@@ -5423,7 +5468,15 @@ def pcolor(self, *args, **kwargs):
                             axis=1)
         verts = xy.reshape((npoly, 5, 2))
 
-        C = compress(ravelmask, ma.filled(C[0:Ny - 1, 0:Nx - 1]).ravel())
+        if isinstance(norm, mcolors.BivariateNorm):
+            C = np.array([
+                            compress(
+                                ravelmask,
+                                ma.filled(c[0:Ny - 1, 0:Nx - 1]).ravel()
+                            ) for c in C
+                        ])
+        else:
+            C = compress(ravelmask, ma.filled(C[0:Ny - 1, 0:Nx - 1]).ravel())
 
         linewidths = (0.25,)
         if 'linewidth' in kwargs:
@@ -5450,9 +5503,12 @@ def pcolor(self, *args, **kwargs):
 
         collection.set_alpha(alpha)
         collection.set_array(C)
-        if norm is not None and not isinstance(norm, mcolors.Normalize):
-            msg = "'norm' must be an instance of 'mcolors.Normalize'"
+
+        if norm is not None and not isinstance(norm, mcolors.Norms):
+            msg = ("'norm' must be an instance of 'mcolors.Normalize' or "
+                   "'mcolors.BivariateNorm'")
             raise ValueError(msg)
+
         collection.set_cmap(cmap)
         collection.set_norm(norm)
         collection.set_clim(vmin, vmax)
@@ -5518,11 +5574,13 @@ def pcolormesh(self, *args, **kwargs):
         Keyword arguments:
 
           *cmap*: [ *None* | Colormap ]
-            A :class:`matplotlib.colors.Colormap` instance. If *None*, use
-            rc settings.
+            A :class:`matplotlib.colors.Colormap` or
+            :class:`matplotlib.colors.BivariateColormap` instance. If *None*,
+            use rc settings.
 
           *norm*: [ *None* | Normalize ]
-            A :class:`matplotlib.colors.Normalize` instance is used to
+            A :class:`matplotlib.colors.Normalize` or
+            :class:`matplotlib.colors.BivariateNorm` instance is used to
             scale luminance data to 0,1. If *None*, defaults to
             :func:`normalize`.
 
@@ -5582,25 +5640,41 @@ def pcolormesh(self, *args, **kwargs):
 
         allmatch = (shading == 'gouraud')
 
-        X, Y, C = self._pcolorargs('pcolormesh', *args, allmatch=allmatch)
+        kw = {'norm': norm, 'cmap': cmap, 'allmatch': allmatch}
+        X, Y, C = self._pcolorargs('pcolormesh', *args, **kw)
         Ny, Nx = X.shape
 
+        is_bivari = (isinstance(norm, mcolors.BivariateNorm) or
+                     isinstance(cmap, mcolors.BivariateColormap))
+        if is_bivari:
+            if C.ndim != 3 and C.shape[0] != 2:
+                raise TypeError("Expected shape like (2, n, m)")
+            if cmap is None:
+                cmap = mcolors.BivariateColormap()
+            if norm is None:
+                norm = mcolors.BivariateNorm()
+
         # unit conversion allows e.g. datetime objects as axis values
         self._process_unit_info(xdata=X, ydata=Y, kwargs=kwargs)
         X = self.convert_xunits(X)
         Y = self.convert_yunits(Y)
 
-        # convert to one dimensional arrays
-        C = C.ravel()
+        # convert to one dimensional arrays if univariate
+        if isinstance(norm, mcolors.BivariateNorm):
+            C = np.asarray([c.ravel() for c in C])
+        else:
+            C = C.ravel()
+
         coords = np.column_stack((X.flat, Y.flat)).astype(float, copy=False)
 
         collection = mcoll.QuadMesh(Nx - 1, Ny - 1, coords,
                                     antialiased=antialiased, shading=shading,
                                     **kwargs)
         collection.set_alpha(alpha)
         collection.set_array(C)
-        if norm is not None and not isinstance(norm, mcolors.Normalize):
-            msg = "'norm' must be an instance of 'mcolors.Normalize'"
+        if norm is not None and not isinstance(norm, mcolors.Norms):
+            msg = ("'norm' must be an instance of 'mcolors.Normalize' or "
+                   "'mcolors.BivariateNorm'")
             raise ValueError(msg)
         collection.set_cmap(cmap)
         collection.set_norm(norm)
@@ -5634,7 +5708,7 @@ def pcolormesh(self, *args, **kwargs):
     @docstring.dedent_interpd
     def pcolorfast(self, *args, **kwargs):
         """
-        pseudocolor plot of a 2-D array
+        pseudocolor plot of a 2-D univariate or 3-D bivariate array
 
         Experimental; this is a pcolor-type method that
         provides the fastest possible rendering with the Agg
@@ -5693,11 +5767,13 @@ def pcolorfast(self, *args, **kwargs):
         Optional keyword arguments:
 
           *cmap*: [ *None* | Colormap ]
-            A :class:`matplotlib.colors.Colormap` instance from cm. If *None*,
-            use rc settings.
+            A :class:`matplotlib.colors.Colormap` or
+            :class:`matplotlib.colors.BivariateColormap` instance from cm.
+            If *None*, use rc settings.
 
           *norm*: [ *None* | Normalize ]
-            A :class:`matplotlib.colors.Normalize` instance is used to scale
+            A :class:`matplotlib.colors.Normalize` or
+            :class:`matplotlib.colors.BivariateNorm` instance is used to scale
             luminance data to 0,1. If *None*, defaults to normalize()
 
           *vmin*/*vmax*: [ *None* | scalar ]
@@ -5723,12 +5799,26 @@ def pcolorfast(self, *args, **kwargs):
         cmap = kwargs.pop('cmap', None)
         vmin = kwargs.pop('vmin', None)
         vmax = kwargs.pop('vmax', None)
-        if norm is not None and not isinstance(norm, mcolors.Normalize):
-            msg = "'norm' must be an instance of 'mcolors.Normalize'"
+
+        if norm is not None and not isinstance(norm, mcolors.Norms):
+            msg = ("'norm' must be an instance of 'mcolors.Normalize' or "
+                   "'mcolors.BivariateNorm'")
             raise ValueError(msg)
 
-        C = args[-1]
-        nr, nc = C.shape
+        C = np.asarray(args[-1])
+
+        is_bivari = (isinstance(norm, mcolors.BivariateNorm) or
+                     isinstance(cmap, mcolors.BivariateColormap))
+        if is_bivari:
+            if C.ndim != 3 and C.shape[0] != 2:
+                raise TypeError("Expected shape like (2, n, m)")
+            if cmap is None:
+                cmap = mcolors.BivariateColormap()
+            if norm is None:
+                norm = mcolors.BivariateNorm()
+            nr, nc = C.shape[1:]
+        else:
+            nr, nc = C.shape
         if len(args) == 1:
             style = "image"
             x = [0, nc]

diff --git a/lib/matplotlib/axes/_subplots.py b/lib/matplotlib/axes/_subplots.py
@@ -177,7 +177,6 @@ def subplot_class_factory(axes_class=None):
                          (SubplotBase, axes_class),
                          {'_axes_class': axes_class})
         _subplot_classes[axes_class] = new_class
-
     return new_class
 
 # This is provided for backward compatibility