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:


CategoryCookbookMatplotlib

Cookbook/Matplotlib/Gridding irregularly spaced data (last edited 2011-09-07 13:23:10 by BertrandLacoste)