Interpolation over an irregular grid

Question

So, I have three numpy arrays which store latitude, longitude, and some property value on a grid -- that is, I have LAT(y,x), LON(y,x), and, say temperature T(y,x), for some limits of x and y. The grid isn't necessarily regular -- in fact, it's tripolar.

I then want to interpolate these property (temperature) values onto a bunch of different lat/lon points (stored as lat1(t), lon1(t), for about 10,000 t...) which do not fall on the actual grid points. I've tried matplotlib.mlab.griddata, but that takes far too long (it's not really designed for what I'm doing, after all). I've also tried scipy.interpolate.interp2d, but I get a MemoryError (my grids are about 400x400).

Is there any sort of slick, preferably fast way of doing this? I can't help but think the answer is something obvious... Thanks!!

Answer 1

I suggest you taking a look at GRASS (an open source GIS package) interpolation features (http://grass.ibiblio.org/gdp/html_grass62/v.surf.bspline.html). It's not in python but you can reimplement it or interface with C code.

Answer 2

Am I right in thinking your data grids look something like this (red is the old data, blue is the new interpolated data)?

This might be a slightly brute-force-ish approach, but what about rendering your existing data as a bitmap (opengl will do simple interpolation of colours for you with the right options configured and you could render the data as triangles which should be fairly fast). You could then sample pixels at the locations of the new points.

Alternatively, you could sort your first set of points spatially and then find the closest old points surrounding your new point and interpolate based on the distances to those points.

Answer 3

There's a bunch of options here, which one is best will depend on your data... However I don't know of an out-of-the-box solution for you

You say your input data is from tripolar data. There are three main cases for how this data could be structured.

1. Sampled from a 3d grid in tripolar space, projected back to 2d LAT, LON data.
2. Sampled from a 2d grid in tripolar space, projected into 2d LAT LON data.
3. Unstructured data in tripolar space projected into 2d LAT LON data

The easiest of these is 2. Instead of interpolating in LAT LON space, "just" transform your point back into the source space and interpolate there.

Another option that works for 1 and 2 is to search for the cells that maps from tripolar space to cover your sample point. (You can use a BSP or grid type structure to speed up this search) Pick one of the cells, and interpolate inside it.

Finally there's a heap of unstructured interpolation options .. but they tend to be slow. A personal favourite of mine is to use a linear interpolation of the nearest N points, finding those N points can again be done with gridding or a BSP. Another good option is to Delauney triangulate the unstructured points and interpolate on the resulting triangular mesh.

Personally if my mesh was case 1, I'd use an unstructured strategy as I'd be worried about having to handle searching through cells with overlapping projections. Choosing the "right" cell would be difficult.

Answer 4

Try the combination of inverse-distance weighting and scipy.spatial.KDTree described in SO inverse-distance-weighted-idw-interpolation-with-python. Kd-trees work nicely in 2d 3d ..., inverse-distance weighting is smooth and local, and the k= number of nearest neighbours can be varied to tradeoff speed / accuracy.