How to set a maximum distance between points for interpolation when using scipy.interpolate.griddata?

Question

I have a spatial set of data with Z values I want to interpolate using some matplotlib or scipy module. My XY points have a concave shape and I don't want interpolated values in the empty zone. Is there a method that easily allow user to set a maximum distance between points to avoid interpolation in the empty zone?

Answer 1

I struggled with the same question and found a work around by re-using the kd-tree implementation that scipy itself uses for the nearest neighbour interpolation, masking the interpolated result array with the result of the kd-tree querying result.

Consider the example code below:

import numpy as np
import scipy.interpolate
import matplotlib.pyplot as plt

# Generate some random data
xy = np.random.random((2**15, 2))
z = np.sin(10*xy[:,0]) * np.cos(10*xy[:,1])

grid = np.meshgrid(
    np.linspace(0, 1, 512),
    np.linspace(0, 1, 512)
)

# Interpolate
result1 = scipy.interpolate.griddata(xy, z, tuple(grid), 'linear')

# Show
plt.figimage(result1)
plt.show()

# Remove rectangular window
mask = np.logical_and.reduce((xy[:,0] > 0.2, xy[:,0] < 0.8, xy[:,1] > 0.2, xy[:,1] < 0.8))
xy, z = xy[~mask], z[~mask]

# Interpolate
result2 = scipy.interpolate.griddata(xy, z, tuple(grid), 'linear')

# Show
plt.figimage(result2)
plt.show()

This generates the following two images. Notices the strong interpolation artefacts because of the missing rectangle window in the centre of the data.

Now if we run the code below on the same example data, the following image is obtained.

THRESHOLD = 0.01

from scipy.interpolate.interpnd import _ndim_coords_from_arrays
from scipy.spatial import cKDTree

# Construct kd-tree, functionality copied from scipy.interpolate
tree = cKDTree(xy)
xi = _ndim_coords_from_arrays(tuple(grid), ndim=xy.shape[1])
dists, indexes = tree.query(xi)

# Copy original result but mask missing values with NaNs
result3 = result2[:]
result3[dists > THRESHOLD] = np.nan

# Show
plt.figimage(result3)
plt.show()

I realize it may not be the visual effect you're after exactly. Especially if your dataset is not very dense you'll need to have a high distance threshold value in order for legitimately interpolated data not to be masked. If your data is dense enough, you might be able to get away with a relatively small radius, or maybe come up with a smarter cut-off function. Hope that helps.