Fast 3D interpolation of atmospheric data in Numpy / Scipy
I am trying to interpolate 3D atmospheric data from one vertical coordinate to another using Numpy / Scipy. For example, I have temperature and relative humidity cubes, both of which are on constant, constant pressure surfaces. I want to interpolate relative humidity to constant temperature surface (s).
The exact problem I am trying to solve was asked earlier here , however the solution is very slow there. In my case, I have roughly 3M points in my cube (30x321x321) and this method takes about 4 minutes to work with one dataset.
This post is almost 5 years old. Newer versions of Numpy / Scipy perhaps have methods that handle this faster? Perhaps new perspectives on the problem are better suited? I am open to suggestions.
EDIT:
Slow = 4 minutes for one set of data cubes. I'm not sure how else I can quantify it.
Code used ...
def interpLevel(grid,value,data,interp='linear'):
"""
Interpolate 3d data to a common z coordinate.
Can be used to calculate the wind/pv/whatsoever values for a common
potential temperature / pressure level.
grid : numpy.ndarray
The grid. For example the potential temperature values for the whole 3d
grid.
value : float
The common value in the grid, to which the data shall be interpolated.
For example, 350.0
data : numpy.ndarray
The data which shall be interpolated. For example, the PV values for
the whole 3d grid.
kind : str
This indicates which kind of interpolation will be done. It is directly
passed on to scipy.interpolate.interp1d().
returns : numpy.ndarray
A 2d array containing the *data* values at *value*.
"""
ret = np.zeros_like(data[0,:,:])
for yIdx in xrange(grid.shape[1]):
for xIdx in xrange(grid.shape[2]):
# check if we need to flip the column
if grid[0,yIdx,xIdx] > grid[-1,yIdx,xIdx]:
ind = -1
else:
ind = 1
f = interpolate.interp1d(grid[::ind,yIdx,xIdx], \
data[::ind,yIdx,xIdx], \
kind=interp)
ret[yIdx,xIdx] = f(value)
return ret
EDIT 2: I could share npy sampled data dumps if anyone was interested enough to see what I am working with.
source to share
Since this is atmospheric data, I assume your mesh is not evenly spaced; however, if your grid is straight-line (such that every vertical column has the same set of z-coordinates), then you have some options.
For example, if you only want linear interpolation (say for simple visualization), you can simply do something like:
# Find nearest grid point
idx = grid[:,0,0].searchsorted(value)
upper = grid[idx,0,0]
lower = grid[idx - 1, 0, 0]
s = (value - lower) / (upper - lower)
result = (1-s) * data[idx - 1, :, :] + s * data[idx, :, :]
(Of course, you will need to add value
out-of-range checks ). For a grid of your size, this will be very fast (as in tiny fractions of a second)
You can easily modify the above to do cubic interpolation if needed; the challenge is to select the correct weights for the uneven vertical distance.
The problem with using scipy.ndimage.map_coordinates is that while it provides higher order interpolation and can handle arbitrary sample points, it assumes that the input will be evenly distributed. It will still get smooth results, but it will not be a reliable approximation.
If your grid is not linear, so the z value for a given index varies for different x and y indices, then the approach you are using now is probably best if you don't get a fair analysis of your particular problem.
UPDATE:
One neat trick (again, assuming each column has the same, not necessarily regular coordinates) is to use interp1d
to extract the weights doing something like the following:
NZ = grid.shape[0]
zs = grid[:,0,0]
ident = np.identity(NZ)
weight_func = interp1d(zs, ident, 'cubic')
You only need to do it above one grid at a time; you can even reuse weight_func as long as the vertical coordinates don't change.
When the time comes to interpolate, weight_func(value)
will give you weights that you can use to compute one interpolated value in (x_idx, y_idx)
with:
weights = weight_func(value)
interp_val = np.dot(data[:, x_idx, y_idx), weights)
If you want to compute the entire plane of the interpolated values, you can use np.inner
though, since your z-coordinate comes first, you would need to do:
result = np.inner(data.T, weights).T
Again, the calculation should be almost immediate.
source to share