How to optimize interpolation of large set of scattered points?

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);
    }

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 distan开发者_JAVA技巧ce 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.

---

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)

---

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.