Search code examples
pythonnumpyinterpolation

Interpolate points over a surface


I have a function f(x,y) that I know for certain values of X and Y. Then, I have a new function x,y= g(t) that produces a new series of points X and Y. for those points, using the pre-calculated data I need interpolate.

import numpy as np
import timeit
import matplotlib.pyplot as plt

def f(X,Y):
    return np.sin(X)+np.cos(Y)

def g(t):
    return 3*np.cos(t)-0.5, np.sin(5*t)-1

x=      np.linspace(-4,4,100)
y=      np.linspace(-4,4,100)
X, Y=   np.meshgrid(x,y)
Z=      f(X,Y)
   
t=      np.linspace(0,1,50)
x_t, y_t=   g(t)
    
plt.figure()
plt.imshow(Z, extent=(-4, 4, -4, 4), origin='lower')
plt.scatter(x_t,y_t)
plt.show()

enter image description here In other words, I need to obtain the shown curve with the previously calculated Z values, doing interpolation since in real life I dont have access to the actual function

Many thanks!

EDIT:

I found a function that does exactly what I want, but it produces the wrongs samples.

import numpy as np
import timeit
import matplotlib.pyplot as plt
from scipy.interpolate import interpn
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

def f(X,Y):
    return np.sin(X)+np.cos(Y)

def g(t):
    return -np.cos(t), np.sin(2*t)

x=      np.linspace(-4,4,1001)
y=      np.linspace(-4,4,1001)
X, Y=   np.meshgrid(x,y)
Z=      f(X,Y)
grid= (x,y) 

t=      np.linspace(0,2*np.pi, 50)
x_t, y_t=   g(t)
xy_t=   np.stack([x_t, y_t], axis=-1)
Z_real= f(x_t, y_t)
Z_inter= interpn(grid, Z, xy_t)

fig= plt.figure()
gs=     fig.add_gridspec(1,2, hspace=0.1, wspace=0.1)
ax1=        fig.add_subplot(gs[0,0],projection="3d")
surf=   ax1.plot_surface(X,Y,Z, cmap=cm.jet, alpha=0.3)
ax1.scatter(x_t, y_t, Z_real)
fig.colorbar(surf, shrink=0.5,  aspect=5)

ax2=    fig.add_subplot(gs[0,1])
ax2.plot(t, Z_real)
ax2.plot(t, Z_inter, '+')

plt.show()

enter image description here

Can anyone tell me what I am doing wrong?


Solution

  • I recently shared a solution to a comparable question that is applicable to your case as well. In brief, scipy.interpolate provides a range of 2D interpolation options. In this instance, I utilized RegularGridInterpolator or CloughTocher2DInterpolator, but you also have the option to employ another methods within scipy.interpolate.

    RegularGridInterpolator outperforms CloughTocher2DInterpolator in terms of speed. However, the latter does not necessitate data to be arranged on a grid.

    import scipy.interpolate
    
    # Option 1: RegularGridInterpolator
    z_reg_interpolator = RegularGridInterpolator(
        (x, y), 
        f(*np.meshgrid(x, y, indexing='ij', sparse=True))
    )
    Z_reg_inter= z_reg_interpolator(np.array(g(t)).T)
    
    # Option 2: CloughTocher2DInterpolator with less datapoints to speed it up
    n = 101
    x_n= np.linspace(-4, 4, n)
    y_n= np.linspace(-4, 4, n)
    XY_n = np.meshgrid(x_n,y_n)
    
    z_interpolator = scipy.interpolate.CloughTocher2DInterpolator(
        np.dstack(XY_n).reshape(-1,2),
        f(*XY_n).reshape(-1,1))
    Z_clo_inter = z_interpolator(np.array(g(t)).T)
    

    enter image description here