A commonly asked question on the matplotlib mailing lists is "how do I make a contour plot of my irregularly spaced data?". The answer is, first you interpolate it to a regular grid. As of version 0.98.3, matplotlib provides a griddata function that behaves similarly to the matlab version. It performs "natural neighbor interpolation" of irregularly spaced data a regular grid, which you can then plot with contour, imshow or pcolor.
Example 1
This requires Scipy 0.9:
1 import numpy as np
2 from scipy.interpolate import griddata
3 import matplotlib.pyplot as plt
4 import numpy.ma as ma
5 from numpy.random import uniform, seed
6 # make up some randomly distributed data
7 seed(1234)
8 npts = 200
9 x = uniform(-2,2,npts)
10 y = uniform(-2,2,npts)
11 z = x*np.exp(-x**2-y**2)
12 # define grid.
13 xi = np.linspace(-2.1,2.1,100)
14 yi = np.linspace(-2.1,2.1,100)
15 # grid the data.
16 zi = griddata((x, y), z, (xi[None,:], yi[:,None]), method='cubic')
17 # contour the gridded data, plotting dots at the randomly spaced data points.
18 CS = plt.contour(xi,yi,zi,15,linewidths=0.5,colors='k')
19 CS = plt.contourf(xi,yi,zi,15,cmap=plt.cm.jet)
20 plt.colorbar() # draw colorbar
21 # plot data points.
22 plt.scatter(x,y,marker='o',c='b',s=5)
23 plt.xlim(-2,2)
24 plt.ylim(-2,2)
25 plt.title('griddata test (%d points)' % npts)
26 plt.show()
Example 2
1 import numpy as np
2 from matplotlib.mlab import griddata
3 import matplotlib.pyplot as plt
4 import numpy.ma as ma
5 from numpy.random import uniform
6 # make up some randomly distributed data
7 npts = 200
8 x = uniform(-2,2,npts)
9 y = uniform(-2,2,npts)
10 z = x*np.exp(-x**2-y**2)
11 # define grid.
12 xi = np.linspace(-2.1,2.1,100)
13 yi = np.linspace(-2.1,2.1,100)
14 # grid the data.
15 zi = griddata(x,y,z,xi,yi)
16 # contour the gridded data, plotting dots at the randomly spaced data points.
17 CS = plt.contour(xi,yi,zi,15,linewidths=0.5,colors='k')
18 CS = plt.contourf(xi,yi,zi,15,cmap=plt.cm.jet)
19 plt.colorbar() # draw colorbar
20 # plot data points.
21 plt.scatter(x,y,marker='o',c='b',s=5)
22 plt.xlim(-2,2)
23 plt.ylim(-2,2)
24 plt.title('griddata test (%d points)' % npts)
25 plt.show()
By default, griddata uses the scikits delaunay package (included in matplotlib) to do the natural neighbor interpolation. Unfortunately, the delaunay package is known to fail for some nearly pathological cases. If you run into one of those cases, you can install the matplotlib natgrid toolkit. Once that is installed, the griddata function will use it instead of delaunay to do the interpolation. The natgrid algorithm is a bit more robust, but cannot be included in matplotlib proper because of licensing issues.
The radial basis function module in the scipy sandbox can also be used to interpolate/smooth scattered data in n dimensions. See Cookbook/RadialBasisFunctions for details.
Example 3
A less robust but perhaps more intuitive method is presented in the code below. This function takes three 1D arrays, namely two independent data arrays and one dependent data array and bins them into a 2D grid. On top of that, the code also returns two other grids, one where each binned value represents the number of points in that bin and another in which each bin contains the indices of the original dependent array which are contained in that bin. These can be further used for interpolation between bins if necessary.
The is essentially an Occam's Razor approach to the matplotlib.mlab griddata function, as both produce similar results.
1 # griddata.py - 2010-07-11 ccampo
2 import numpy as np
3
4 def griddata(x, y, z, binsize=0.01, retbin=True, retloc=True):
5 """
6 Place unevenly spaced 2D data on a grid by 2D binning (nearest
7 neighbor interpolation).
8
9 Parameters
10 ----------
11 x : ndarray (1D)
12 The idependent data x-axis of the grid.
13 y : ndarray (1D)
14 The idependent data y-axis of the grid.
15 z : ndarray (1D)
16 The dependent data in the form z = f(x,y).
17 binsize : scalar, optional
18 The full width and height of each bin on the grid. If each
19 bin is a cube, then this is the x and y dimension. This is
20 the step in both directions, x and y. Defaults to 0.01.
21 retbin : boolean, optional
22 Function returns `bins` variable (see below for description)
23 if set to True. Defaults to True.
24 retloc : boolean, optional
25 Function returns `wherebins` variable (see below for description)
26 if set to True. Defaults to True.
27
28 Returns
29 -------
30 grid : ndarray (2D)
31 The evenly gridded data. The value of each cell is the median
32 value of the contents of the bin.
33 bins : ndarray (2D)
34 A grid the same shape as `grid`, except the value of each cell
35 is the number of points in that bin. Returns only if
36 `retbin` is set to True.
37 wherebin : list (2D)
38 A 2D list the same shape as `grid` and `bins` where each cell
39 contains the indicies of `z` which contain the values stored
40 in the particular bin.
41
42 Revisions
43 ---------
44 2010-07-11 ccampo Initial version
45 """
46 # get extrema values.
47 xmin, xmax = x.min(), x.max()
48 ymin, ymax = y.min(), y.max()
49
50 # make coordinate arrays.
51 xi = np.arange(xmin, xmax+binsize, binsize)
52 yi = np.arange(ymin, ymax+binsize, binsize)
53 xi, yi = np.meshgrid(xi,yi)
54
55 # make the grid.
56 grid = np.zeros(xi.shape, dtype=x.dtype)
57 nrow, ncol = grid.shape
58 if retbin: bins = np.copy(grid)
59
60 # create list in same shape as grid to store indices
61 if retloc:
62 wherebin = np.copy(grid)
63 wherebin = wherebin.tolist()
64
65 # fill in the grid.
66 for row in range(nrow):
67 for col in range(ncol):
68 xc = xi[row, col] # x coordinate.
69 yc = yi[row, col] # y coordinate.
70
71 # find the position that xc and yc correspond to.
72 posx = np.abs(x - xc)
73 posy = np.abs(y - yc)
74 ibin = np.logical_and(posx < binsize/2., posy < binsize/2.)
75 ind = np.where(ibin == True)[0]
76
77 # fill the bin.
78 bin = z[ibin]
79 if retloc: wherebin[row][col] = ind
80 if retbin: bins[row, col] = bin.size
81 if bin.size != 0:
82 binval = np.median(bin)
83 grid[row, col] = binval
84 else:
85 grid[row, col] = np.nan # fill empty bins with nans.
86
87 # return the grid
88 if retbin:
89 if retloc:
90 return grid, bins, wherebin
91 else:
92 return grid, bins
93 else:
94 if retloc:
95 return grid, wherebin
96 else:
97 return grid
The following example demonstrates a usage of this method.
1 import numpy as np
2 import matplotlib.pyplot as plt
3 import griddata
4
5 npr = np.random
6 npts = 3000. # the total number of data points.
7 x = npr.normal(size=npts) # create some normally distributed dependent data in x.
8 y = npr.normal(size=npts) # ... do the same for y.
9 zorig = x**2 + y**2 # z is a function of the form z = f(x, y).
10 noise = npr.normal(scale=1.0, size=npts) # add a good amount of noise
11 z = zorig + noise # z = f(x, y) = x**2 + y**2
12
13 # plot some profiles / cross-sections for some visualization. our
14 # function is a symmetric, upward opening paraboloid z = x**2 + y**2.
15 # We expect it to be symmetric about and and y, attain a minimum on
16 # the origin and display minor Gaussian noise.
17
18 plt.ion() # pyplot interactive mode on
19
20 # x vs z cross-section. notice the noise.
21 plt.plot(x, z, '.')
22 plt.title('X vs Z=F(X,Y=constant)')
23 plt.xlabel('X')
24 plt.ylabel('Z')
25
26 # y vs z cross-section. notice the noise.
27 plt.plot(y, z, '.')
28 plt.title('Y vs Z=F(Y,X=constant)')
29 plt.xlabel('Y')
30 plt.ylabel('Z')
31
32 # now show the dependent data (x vs y). we could represent the z data
33 # as a third axis by either a 3d plot or contour plot, but we need to
34 # grid it first....
35 plt.plot(x, y, '.')
36 plt.title('X vs Y')
37 plt.xlabel('X')
38 plt.ylabel('Y')
39
40 # enter the gridding. imagine drawing a symmetrical grid over the
41 # plot above. the binsize is the width and height of one of the grid
42 # cells, or bins in units of x and y.
43 binsize = 0.3
44 grid, bins, binloc = griddata.griddata(x, y, z, binsize=binsize) # see this routine's docstring
45
46
47 # minimum values for colorbar. filter our nans which are in the grid
48 zmin = grid[np.where(np.isnan(grid) == False)].min()
49 zmax = grid[np.where(np.isnan(grid) == False)].max()
50
51 # colorbar stuff
52 palette = plt.matplotlib.colors.LinearSegmentedColormap('jet3',plt.cm.datad['jet'],2048)
53 palette.set_under(alpha=0.0)
54
55 # plot the results. first plot is x, y vs z, where z is a filled level plot.
56 extent = (x.min(), x.max(), y.min(), y.max()) # extent of the plot
57 plt.subplot(1, 2, 1)
58 plt.imshow(grid, extent=extent, cmap=palette, origin='lower', vmin=zmin, vmax=zmax, aspect='auto', interpolation='bilinear')
59 plt.xlabel('X values')
60 plt.ylabel('Y values')
61 plt.title('Z = F(X, Y)')
62 plt.colorbar()
63
64 # now show the number of points in each bin. since the independent data are
65 # Gaussian distributed, we expect a 2D Gaussian.
66 plt.subplot(1, 2, 2)
67 plt.imshow(bins, extent=extent, cmap=palette, origin='lower', vmin=0, vmax=bins.max(), aspect='auto', interpolation='bilinear')
68 plt.xlabel('X values')
69 plt.ylabel('Y values')
70 plt.title('X, Y vs The No. of Pts Per Bin')
71 plt.colorbar()
The binned data: 
Raw data superimposed on top of binned data: 