Matplotlib's pcolor is not like MatlabTM pcolor, and that leads to some confusion. In fact, matplotlib's pcolor is an improved way to "pcolor" things.

Here are some plots to illustrate what I mean. First, some fake data

In [2]:

import numpy as np

data = np.array([[1, 2, 3], [4, 5, 6]])
x = np.array([1, 2, 3])
y = np.array([1, 2])

now let's plot the data passing the x and y coordinates to pcolor:

In [3]:

import matplotlib.pyplot as plt

fig, ax = plt.subplots()
cs = ax.pcolor(x, y, data)

OK, that is identical to Matlab's pcolor. But wait, what if we omit the coordinates?

In [4]:

fig, ax = plt.subplots()
cs = ax.pcolor(data)

Aha! That is more like it!

The secret is hidden in plain site in the docstring,

"""
X and Y, if given, specify the (x, y) coordinates of the colored quadrilaterals; the quadrilateral for C[i,j] has corners at:

(X[i,   j],   Y[i,   j]),
(X[i,   j+1], Y[i,   j+1]),
(X[i+1, j],   Y[i+1, j]),
(X[i+1, j+1], Y[i+1, j+1]).

Ideally the dimensions of X and Y should be one greater than those of C; if the dimensions are the same, then the last row and column of C will be ignored.
"""

That means you can say goodbye to hacked versions of MatlabTM pcolor like pcolorjw.m and variants.

If we pass x and y with dimensions greater than C, like the docstring says, we can get a proper plot as expected:

(Note that something like this in Matlab would result in an error ;-)

In [5]:

import matplotlib.pyplot as plt

fig, ax = plt.subplots()
cs = ax.pcolor(np.r_[x, 4], np.r_[y, 3], data)

If your data is an iris cube you will get proper pcolor rendering for free plus the right coordinate ticklabels, if the coordinate bounds are defined.

In [6]:

import numpy as np
from iris.cube import Cube
from iris.coords import DimCoord


def create_cube():
    data = np.array([[1, 2, 3], [4, 5, 6]])
    cube = Cube(data)
    
    lon = DimCoord(np.array([1, 2, 3]), standard_name='longitude',
                   units='degrees', circular=False)
    lat = DimCoord(np.array([1, 2]), standard_name='latitude',
                   units='degrees')
    cube.add_dim_coord(lon, 1)
    cube.add_dim_coord(lat, 0)
    return cube

cube = create_cube()

In [7]:

x = cube.coord(axis='X')
x.guess_bounds()
x

Out[7]:

DimCoord(array([1, 2, 3]), bounds=array([[ 0.5,  1.5],
       [ 1.5,  2.5],
       [ 2.5,  3.5]]), standard_name='longitude', units=Unit('degrees'))

In [8]:

y = cube.coord(axis='Y')
y.guess_bounds()
y

Out[8]:

DimCoord(array([1, 2]), bounds=array([[ 0.5,  1.5],
       [ 1.5,  2.5]]), standard_name='latitude', units=Unit('degrees'))

That works with plain iris,

In [9]:

import iris.quickplot as qplt

cs = qplt.pcolormesh(cube)

or cartopy with cartopy:

In [10]:

import cartopy.crs as ccrs
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter

fig, ax = plt.subplots(subplot_kw=dict(projection=ccrs.PlateCarree()))
cs = qplt.pcolormesh(cube)
ax.set_xticks(x.points, crs=ccrs.PlateCarree())
ax.set_yticks(y.points, crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)

Note that when we added the extra point in x and y previously we changed the center of the quadrilateral. The correct x and y coordinates in that example should be something like the guessed bounds from iris:

In [11]:

xx = np.unique(x.bounds)
yy = np.unique(y.bounds)

fig, ax = plt.subplots()
cs = ax.pcolor(xx, yy, cube.data)

In [12]:

HTML(html)

Out[12]:

This post was written as an IPython notebook. It is available for download or as a static html.

python4oceanographers by Filipe Fernandes is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
Based on a work at https://ocefpaf.github.io/.

python4oceanographers

Turning ripples into waves

To Pcolor or not to Pcolor

Comments