How to optimize interpolation of large set of scattered points?

Question

I am currently working with a set of coordinate points (longitude, latitude, about 60000 of them) and the temperature at that location. I need to do a interpolation on them to compute the values at some points with unknown temperature as to map certain regions. As to respect the influence that the points have between them I have converted every (long, lat) point to a unit sphere point (x, y, z). I have started applying the generalized multidimension Shepard interpolation from "Numerical recipes 3rd Edition":

Doub interp(VecDoub_I &pt) 
    {
        Doub r, w, sum=0., sumw=0.;
        if (pt.size() != dim) 
            throw("RBF_interp bad pt size");
        for (Int i=0;i<n;i++) 
        {
            if ((r=rad(&pt[0],&pts[i][0])) == 0.) 
                return vals[i];
            sum += (w = pow(r,pneg));
            sumw += w*vals[i];
        }
        return sumw/sum;
    }

    Doub rad(const Doub *p1, const Doub *p2) 
    {
        Doub sum = 0.;
        for (Int i=0;i<dim;i++) 
            sum += SQR(p1[i]-p2[i]);
        return sqrt(sum);
    }

enter image description here As you can see, for the interpolation of one point, the algorithm computes the distance of that point to each of the other points and taking it as a weight in the final value. Even though this algorithm works it is much too slow compared to what I need since I will be computing a lot of points to map a grid of a certain region. One way of optimizing this is that I could leave out the points than are beyond a certain radius, but would pose a problem for areas with few or no points. Another thing would be to reduce the computing of the distance between each 2 points by only computing once a Look-up Table and storing the distances. The problem with this is that it is impossible to store such a large matrix (60000 x 60000). The grid of temperatures that is obtained, will be used to compute contours for different temperature values. If anyone knows a way to optimize this algorithm or maybe help with a better one, I will appreciate it.

You could get rid of the square root and replace `pneg` by `2 * pneg`. — Kerrek SB, Oct 10 '11 at 10:09

score 1 · Accepted Answer · edited May 23 '17 at 12:13

Radial basis functions with infinite support is probably not what you want to be using if you have a large number of data points and will be taking a large number of interpolation values.

There are variants that use N nearest neighbours and finite support to reduce the number of points that must be considered for each interpolation value. A variant of this can be found in the first solution mentioned here Inverse Distance Weighted (IDW) Interpolation with Python. (though I have a nagging suspicion that this implementation can be discontinuous under certain conditions - there are certainly variants that are fine)

score 0 · Answer 2 · answered Oct 10 '11 at 10:18

Your look-up table doesn't have to store every point in the 60k square, only those once which are used repeatedly. You can map any coordinate x to int(x*resolution) to improve the hit rate by lowering the resolution.

A similar lookup table for the power function might also help.

How to optimize interpolation of large set of scattered points?

2 Answers2