Kilian Liss
Kilian Liss

Reputation: 67

Python fast interpolation with consistent input shapes

I have some data that is on a consistent grid like the one bellow, where each point has some scalar value:

Input shape

I'd like to interpolate the scalars onto a regular grid like the one bellow:

Output shape

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)

enter image description here

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

Answers (1)

ramzeek
ramzeek

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.

replace

Upvotes: 1

Related Questions