Reputation: 67
I have some data that is on a consistent grid like the one bellow, where each point has some scalar value:
I'd like to interpolate the scalars onto a regular grid like the one bellow:
Currently I am using scipy.interpolate.griddata, but it is a bit slow. I have 162000 data files that need to be interpolated in the exact same way (and plenty more in the future), and so a much faster algorithm would be desirable. All data files have scalars that need to be interpolated with the exact same grid structure.
A simple toy example is provided bellow:
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import griddata
inputGrid = np.array([
[-9.00000000e+01, 4.92919128e-07],
[-9.00000000e+01, 6.34349493e+01],
[ 9.00000000e+01, -6.34349493e+01],
[ 9.00000000e+01, -4.92919128e-07],
[ 2.77323013e+01, -1.60450574e+01],
[-2.77323013e+01, 1.60450574e+01],
[ 1.52267699e+02, -1.60450574e+01],
[-1.52267699e+02, 1.60450574e+01],
[ 1.39613822e+02, 4.63530726e+01],
[ 4.03861780e+01, 4.63530726e+01],
[-1.39613822e+02, -4.63530726e+01],
[-4.03861780e+01, -4.63530726e+01]])
randomScalar = np.random.random(inputGrid.shape[0])
outputGrid = np.array([[i, j] for i in np.arange(-180, 180, 30) for j in np.arange(-90, 90, 20)])
interpolatedScalar = griddata(inputGrid, randomScalar, outputGrid, method='nearest')
plt.scatter(inputGrid[:, 0], inputGrid[:, 1], c=randomScalar, s=100)
plt.scatter(outputGrid[:, 0], outputGrid[:, 1], c=interpolatedScalar, s=10)
I feel like we could take advantage of the consistent shape that the data needs to be interpolated from to significantly speed things up, but I don't know if there are any python functions that can make use of this advantage.
Specifically, I am converting scalar data from an Icosphere onto a UV sphere. Perhaps there are functions out there that already do this? I've tried 3D interpolation algorithms, but they are much slower than the griddata approach.
Any suggestions?
Upvotes: 2
Views: 237
Reputation: 2315
You are using a nearest neighbor interpolation, and so my answer below applies to a nearest neighbor interpolation approach. One could probably apply the technique used to find the right methods for other interpolation approaches.
If you look at the source for scipy.interpolate.griddata
, you can see at line 256 that for the nearest neighbor interpolation, scipy
uses a class called NearestNDInterpolator
(in the same file). In turn, this class gets the nearest neighbors using scipy.spatial.KDTree
. We are going to use that class to get the indices of the points for each grid point. Following the example in the documentation and applying it to your toy model, we can get the points as
tree = KDTree(inputGrid)
dd, ii = tree.query(outputGrid, k=1)
where dd
is the distance (which you don't need for your application) and ii
is the index in inputGrid
.
We can now alter your second scatter
command using these indices:
plt.scatter(outputGrid[:, 0], outputGrid[:, 1], c=randomScalar[ii], s=10)
So, if you are able to predictably break up your input and output grid into sections, you can keep reapplying the indices in ii
for each section.
Upvotes: 1