Scipy 2D interpolation not accomodating every point

Question

Context: I have a sparse 2D array representing measures of ice thickness along flight traces. However, because of instrument errors, intersects between traces can have different values. I want to interpolate the 2D grid between the traces.

Problem: I am comparing 2 interpolations from Scipy: griddata and CloughTocher2DInterpolator. I would like to introduce some tolerance in the interpolation: if two points close together have vastly different values, then the interpolation should not try to fit them exactly.

Question: Is there a way to use any of the two functions to that extent ? I know CloughTocher2DInterpolator has the tol argument but it does not remove artifacts from the interpolation.

Real Case: I plotted below the two interpolations. As you can notice, the CloughTocher2DInterpolator has weirdness induced from the overlapping traces, while the linear interpolation does not. Eventually, I would like something halfway between the two methods: no apparent weirdness, but also a better interpolation than just linear fit (because topography is better represented by splines than straight lines).

Example Code:

import numpy as np
import matplotlib.pyplot as plt

# Create a 100x100 array of NaNs
array = np.full((100, 100), np.nan)

# Store indices of points along the lines
line_indices = []
vals = []

# Generate random lines with sinusoidal shape
num_lines = 40
for _ in range(num_lines):
    x_start = np.random.randint(0, 100)
    y_start = np.random.randint(0, 100)
    length = np.random.randint(20, 50)
    angle = np.random.uniform(0, 2*np.pi)
    values = np.linspace(0, 400, length) + np.random.normal(0, 50, length)
    line_x = (np.arange(length) * np.cos(angle)).astype(int) + x_start
    line_y = (np.arange(length) * np.sin(angle)).astype(int) + y_start
    line_x = np.clip(line_x, 0, 99)
    line_y = np.clip(line_y, 0, 99)
    line_indices.append((line_y, line_x))  # Store indices
    array[line_y, line_x] = values
    vals.append(values)

# Plot the array
plt.imshow(array, cmap='jet')
plt.colorbar(label='Value')
plt.title('Random Lines with Sinusoidal Shape')
plt.show()


# Prepare the interpolation
line_indices = np.hstack(line_indices)
indices_x = line_indices[1]
indices_y = line_indices[0]
Y = np.arange(0, array.shape[0])
X = np.arange(0, array.shape[1])
X, Y = np.meshgrid(X,Y)

# Interpolation
interp = CloughTocher2DInterpolator(list(zip(indices_x, indices_y)), np.hstack(vals))
Z = interp(X, Y) # Interpolated bed 

plt.imshow(Z)