Finding n nearest data points to grid locations
===

I'm working on a problem where I have a large set (>4 million) of data points located in a three-dimensional space, each with a scalar function value.  This is represented by four arrays: XD, YD, ZD, and FD.  The tuple (XD[i], YD[i], ZD[i]) refers to the location of data point i, which has a value of FD[i].

I'd like to superimpose a rectilinear grid of, say, 100x100x100 points in the same space as my data.  This grid is set up as follows.

    [XGrid, YGrid, ZGrid] = np.mgrid[Xmin:Xmax:Xstep, Ymin:Ymax:Ystep, Zmin:Zmax:Zstep]
    XG = XGrid[:,0,0]
    YG = YGrid[0,:,0]
    ZG = ZGrid[0,0,:]

XGrid is a 3D array of the x-value at each point in the grid.  XG is a 1D array of the x-values going from Xmin to Xmax, separated by a distance of XStep.

I'd like to use an interpolation algorithm I have to find the value of the function at each grid point based on the data surrounding it.  In this algorithm I require 20 data points closest (or at least close) to my grid point of interest.  That is, for grid point (XG[i], YG[j], ZG[k]) I want to find the 20 closest data points.

The only way I can think of is to have one for loop that goes through each data point and a subsequent embedded for loop going through all (so many!) data points, calculating the Euclidean distance, and picking out the 20 closest ones.

    for i in range(0,XG.shape):
      for j in range(0,YG.shape):
        for k in range(0,ZG.shape):
        
          Distance = np.zeros([XD.shape])

          for a in range(0,XD.shape):
            Distance[a] = (XD[a] - XG[i])**2 + (YD[a] - YG[j])**2 + (ZD[a] - ZG[k])**2

          B = np.zeros([20], int)
          for a in range(0,20):
            indx = np.argmin(Distance)
            B[a] = indx
            Distance[indx] = float(inf)

This would give me an array, B, of the indices of the data points closest to the grid point.  I feel like this would take too long to go through each data point at each grid point.

I'm looking for any suggestions, such as how I might be able to organize the data points before calculating distances, which could cut down on computation time.