I have several values that are defined on the same irregular grid (x, y, z)
that I want to interpolate onto a new grid (x1, y1, z1)
. i.e., I have f(x, y, z), g(x, y, z), h(x, y, z)
and I want to calculate f(x1, y1, z1), g(x1, y1, z1), h(x1, y1, z1)
.
At the moment I am doing this using scipy.interpolate.griddata
and it works well. However, because I have to perform each interpolation separately and there are many points, it is quite slow, with a great deal of duplication in the calculation (i.e finding which points are closest, setting up the grids etc...).
Is there a way to speedup the calculation and reduce the duplicated calculations? i.e something along the lines of defining the two grids, then changing the values for the interpolation?
There are several things going on every time you make a call to
scipy.interpolate.griddata
:sp.spatial.qhull.Delaunay
is made to triangulate the irregular grid coordinates.The first three steps are identical for all your interpolations, so if you could store, for each new grid point, the indices of the vertices of the enclosing simplex and the weights for the interpolation, you would minimize the amount of computations by a lot. This is unfortunately not easy to do directly with the functionality available, although it is indeed possible:
The function
interp_weights
does the calculations for the first three steps I listed above. Then the functioninterpolate
uses those calcualted values to do step 4 very fast:So first, it does the same as
griddata
, which is good. Second, setting up the interpolation, i.e. computingvtx
andwts
takes roughly the same as a call togriddata
. But third, you can now interpolate for different values on the same grid in virtually no time.The only thing that
griddata
does that is not contemplated here is assigningfill_value
to points that have to be extrapolated. You could do that by checking for points for which at least one of the weights is negative, e.g.:Great thanks to Jaime for his solution (even if I don't really understand how the barycentric computation is done ...)
Here you will find an example adapted from his case in 2D :
It is possible to applied image transformation such as image mapping with a udge speed-up
You can't use the same function definition as the new coordinates will change at every iteration but you can compute triangulation Once for all.
On my laptop the speed-up is between 20 and 40x !
Hope that can help someone
You can try to use Pandas, as it provides high-performance data structures.
It is true that the interpolation method is a wrapper of the scipy interpolation BUT maybe with the improved structures you obtain better speed.
interpolate()
fills the NaN values in the Panel dataset using different methods. Hope it is faster than Scipy.If it doesn't work, there is one way to improve the performance (instead of using a parallelized version of your code): use Cython and implement small routine in C to use inside your Python code. Here you have an example about this.