Python : 2d contour plot from 3 lists : x, y and rho?

You need to interpolate your rho values. There's no one way to do this, and the "best" method depends entirely on the a-priori information you should be incorporating into the interpolation.

Before I go into a rant on "black-box" interpolation methods, though, a radial basis function (e.g. a "thin-plate-spline" is a particular type of radial basis function) is often a good choice. If you have millions of points, this implementation will be inefficient, but as a starting point:

import numpy as npimport matplotlib.pyplot as pltimport scipy.interpolate# Generate data:x, y, z = 10 * np.random.random((3,10))# Set up a regular grid of interpolation pointsxi, yi = np.linspace(x.min(), x.max(), 100), np.linspace(y.min(), y.max(), 100)xi, yi = np.meshgrid(xi, yi)# Interpolaterbf = scipy.interpolate.Rbf(x, y, z, function='linear')zi = rbf(xi, yi)plt.imshow(zi, vmin=z.min(), vmax=z.max(), origin='lower',           extent=[x.min(), x.max(), y.min(), y.max()])plt.scatter(x, y, c=z)plt.colorbar()plt.show()

You can use scipy's griddata (requires Scipy >= 0.10), it's a triangulation-based method.

import numpy as npimport matplotlib.pyplot as pltimport scipy.interpolate# Generate data: for N=1e6, the triangulation hogs 1 GB of memoryN = 1000000x, y = 10 * np.random.random((2, N))rho = np.sin(3*x) + np.cos(7*y)**3# Set up a regular grid of interpolation pointsxi, yi = np.linspace(x.min(), x.max(), 300), np.linspace(y.min(), y.max(), 300)xi, yi = np.meshgrid(xi, yi)# Interpolate; there's also method='cubic' for 2-D data such as herezi = scipy.interpolate.griddata((x, y), rho, (xi, yi), method='linear')plt.imshow(zi, vmin=rho.min(), vmax=rho.max(), origin='lower',           extent=[x.min(), x.max(), y.min(), y.max()])plt.colorbar()plt.show()

There's also inverse distance weighed interpolation -- similar to RBF, but should work better for large # of points: Inverse Distance Weighted (IDW) Interpolation with Python