32

I have a simple problem in python and matplotlib. I have 3 lists : x, y and rho with rho[i] a density at the point x[i], y[i]. All values of x and y are between -1. and 1. but they are not in a specific order.

How to make a contour plot (like with imshow) of the density rho (interpolated at the points x, y).

Thank you very much.

EDIT : I work with large arrays : x, y and rho have between 10,000 and 1,000,000 elements

Vincent
  • 57,703
  • 61
  • 205
  • 388
  • did the code you accepted worked for you ? 'm having same-sort of list scenario but not being able to solve it. – diffracteD May 05 '15 at 11:32

2 Answers2

52

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 np
import matplotlib.pyplot as plt
import scipy.interpolate

# Generate data:
x, y, z = 10 * np.random.random((3,10))

# Set up a regular grid of interpolation points
xi, yi = np.linspace(x.min(), x.max(), 100), np.linspace(y.min(), y.max(), 100)
xi, yi = np.meshgrid(xi, yi)

# Interpolate
rbf = 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()

enter image description here

Joe Kington
  • 275,208
  • 71
  • 604
  • 463
  • Very interesting method but it doesn't work with large arrays : x, y and rho have around 10000 elements... – Vincent Jan 25 '12 at 21:15
  • You can still use an Rbf for large arrays, you just need to only include nearby points. I'll add an example in just a bit. Alternately, if you don't want to actually sample everything on a regular grid, you can use delaunay triangulation to draw the contours (which is just a very simple and not particularly smooth form of interpolation). With that many points, however, it's more practical to just use a local interpolation method. – Joe Kington Jan 25 '12 at 21:27
  • @JoeKington Hi, I'm having a problem with this above code, My data-set consists of lists x,y and z. x and y vary independently, z varies depending on (x,y). `x = (1.2 to 2.5)`, `y=(90 to 180)` and `z=(5 to -5)`. If I try the above code with my dataset i'm getting a collapsed-plot(nothing along x-axis). Please help. – diffracteD May 05 '15 at 11:36
  • 2
    @diffracteD - I'm guessing, but it may be the way you're plotting the output. By default, `imshow` will force the aspect ratio of the plot to be 1. In other words, one centimenter in the x-direction is the same number of units as one centimenter in the y-direction. This will force the axes to be very long and narrow with your data ranges. You're probably getting reasonable output, but plotting it so that it's difficult to see. Try passing `aspect="auto"` to `imshow`. – Joe Kington May 05 '15 at 12:34
  • 1
    @diffracteD - Also be aware that you may want to rescale your data ranges before interpolating. E.g. see the second and third figures in http://stackoverflow.com/a/3867302/325565 Because the x and y ranges of your data only vary by a factor of 10, though, you won't have particularly severe anisotropy problems. You can probably ignore this, but it's good to be aware of. – Joe Kington May 05 '15 at 12:38
  • @JoeKington The `aspect` worked fine for me, the graph went wide. But I need to limit the z-axis from `4 to -3`, any help on that would be very nice. thank you. – diffracteD May 05 '15 at 13:27
  • @diffracteD - Do you mean on the plotted image? That's the `vmin` and `vmax` kwargs to `imshow`. – Joe Kington May 05 '15 at 13:30
  • @JoeKington hi, I did that. But still the plot is not making any sense to me. The contour-color-gradient is not getting generated for me although the scatter plot is denser somewhere and lighter somewhere. I don't understand. Is there any way I can share results 'm getting ? – diffracteD May 05 '15 at 13:46
  • Ah, if you're running the code above more or less as-is, the colorbar will based on the output of `scatter`. Pass the same vmin and vmax kwargs to `scatter` as well as `imshow`. – Joe Kington May 05 '15 at 13:49
  • @JoeKington My scatter-dots are coming colored as per colorbar. But the background-contour is something like 'added citric acid with milk'. Although 'm using nearly same code you've written. any ideas why the-fill-color gradient is lost. ? – diffracteD May 05 '15 at 13:58
14

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

import numpy as np
import matplotlib.pyplot as plt
import scipy.interpolate

# Generate data: for N=1e6, the triangulation hogs 1 GB of memory
N = 1000000
x, y = 10 * np.random.random((2, N))
rho = np.sin(3*x) + np.cos(7*y)**3

# Set up a regular grid of interpolation points
xi, 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 here
zi = 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

Community
  • 1
  • 1
pv.
  • 33,875
  • 8
  • 55
  • 49