How to plot nice 2d density plots of samples in python

contour.py

import numpy as np
import scipy.special as sp
from matplotlib.colors import LinearSegmentedColormap
import matplotlib.cm as cm
import matplotlib.pyplot as pl

def reshist2d(x, y, *args, **kwargs):
    ax = kwargs.pop('ax', pl.gca())

    extent = kwargs.pop('extent', [[x.min(), x.max()], [y.min(), y.max()]])
    bins   = kwargs.pop('bins', 50)
    color  = kwargs.pop('color', 'k')

    cmap = cm.get_cmap('gray')
    cmap._init()
    cmap._lut[:-3,:-1] = 0.
    cmap._lut[:-3, -1] = np.linspace(1,0,cmap.N)

    X = np.linspace(extent[0][0], extent[0][1], bins+1)
    Y = np.linspace(extent[1][0], extent[1][1], bins+1)
    H, X, Y = np.histogram2d(x.flatten(), y.flatten(), bins=(X,Y))

    V = sp.erf(np.arange(0.5, 2.1, 0.5)/np.sqrt(2))
    Hflat = H.flatten()
    inds = np.argsort(Hflat)[::-1]
    Hflat = Hflat[inds]
    sm = np.cumsum(Hflat)
    sm /= sm[-1]

    for i, v0 in enumerate(V):
        try:
            V[i] = Hflat[sm <= v0][-1]
        except:
            V[i] = Hflat[0]

    X, Y = 0.5*(X[1:]+X[:-1]), 0.5*(Y[1:]+Y[:-1])

    ax.plot(x, y, 'o', color=color, ms=1.5, zorder=-1, alpha=0.1,
            rasterized=True)
    ax.contourf(X, Y, H.T, [V[-1], 0.],
        cmap=LinearSegmentedColormap.from_list('cmap',([1]*3,[1]*3),N=2))
    ax.pcolor(X, Y, H.max()-H.T, cmap=cmap)
    ax.contour(X, Y, H.T, V, colors=color)

Dude, are you going to merge this with the plotting in pappy, or are we going to make a separate plotting package built on matplotlib for plotting samples and PDFs inferred from samples?