5

I came across this really cool demonstration of CloughTocher2DInterpolator as part of the response to this question: How to get a non-smoothing 2D spline interpolation with scipy, and am interested in getting the derivatives along the x and y axis. According to scipy.org, the extrapolant is continuously differentiable.

Couldn't find any help in the documentation. Getting dx and dy seems more straightforward with scipy's bivariatesplines but cloughTocher seems to fit my purpose best.

Any help appreciated!

Here is the example code from the link:

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

import numpy as np
from scipy.interpolate import CloughTocher2DInterpolator
from scipy.spatial import Delaunay

# Example unstructured mesh:
nodes = np.array([[-1.        , -1.        ],
       [ 1.        , -1.        ],
       [ 1.        ,  1.        ],
       [-1.        ,  1.        ],
       [ 0.        ,  0.        ],
       [-1.        ,  0.        ],
       [ 0.        , -1.        ],
       [-0.5       ,  0.        ],
       [ 0.        ,  1.        ],
       [-0.75      ,  0.4       ],
       [-0.5       ,  1.        ],
       [-1.        , -0.6       ],
       [-0.25      , -0.5       ],
       [-0.5       , -1.        ],
       [-0.20833333,  0.5       ],
       [ 1.        ,  0.        ],
       [ 0.5       ,  1.        ],
       [ 0.36174242,  0.44412879],
       [ 0.5       , -0.03786566],
       [ 0.2927264 , -0.5411368 ],
       [ 0.5       , -1.        ],
       [ 1.        ,  0.5       ],
       [ 1.        , -0.5       ]])

# Theoretical function:
def F(x, y):
    return x + y -  x*y - (x*y)**2 - 2*x*y**2 + x**2*y + 3*np.exp( -((x+1)**2 + (y+1)**2)*5 )

z = F(nodes[:, 0], nodes[:, 1])

# Finer regular grid:
N2 = 19
x2, y2 = np.linspace(-1, 1, N2), np.linspace(-1, 1, N2)
X2, Y2 = np.meshgrid(x2, y2)

# Interpolation:
tri = Delaunay(nodes)
CT_interpolator = CloughTocher2DInterpolator(tri, z)
z_interpolated = CT_interpolator(X2, Y2)

# Plot
fig = plt.figure(1, figsize=(8,14))

ax = fig.add_subplot(311, projection='3d')
ax.scatter3D(nodes[:, 0], nodes[:, 1], z, s=15, color='red', label='points')

ax.plot_wireframe(X2, Y2, z_interpolated, color='black', label='interpolated')
plt.legend();
annon
  • 63
  • 4

0 Answers0