python4oceanographers

Turning ripples into waves

To Pcolor or not to Pcolor

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.

Creative Commons License
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/.

Comments