保持 2D 插值在形状内

Question

我有一组测量点，我想在它们之间进行插值，为此我使用 SciPy 的 griddata()：

import matplotlib.pyplot as plt 
import numpy as np
import pandas as pd 
from scipy.interpolate import griddata 

data = pd.DataFrame({
    'time': [0, 1, 2, 3, 4, 0.1, 0.9, 2, 3.05, 4, 0, 1, 2.2, 3, 3.95], 
    'force': [1, 2, 4, 9, 16, 0, 0, 0, 0, 0, -1, -2, -4, -9, -16]
})

Times, Forces = np.meshgrid(
    np.linspace(0, 4, 100),
    np.linspace(-16, 16, 100)
)

data['work'] = data['time'] * data['force']

interpolation = griddata(
    (data['time'], data['force']), 
    data['work'], 
    (Times, Forces), 
    method= 'linear'
)

fig, ax = plt.subplots() 

contour = ax.contourf(
    Times, Forces, interpolation
)
ax.scatter(data['time'], data['force'])

fig.show()

我的问题是，我的测量点已经遵循物理可能性的边界，但是插值仍然会针对每个可跨区域进行插值，包括那些无法测量的区域。

如何将插值或至少插值图限制在外部点的形状“范围内”？不幸的是，“时间”测量值有小的偏差。

任何建议将不胜感激！

Answer 1

首先，您需要获取所有点的凹多边形。其次，使用多边形来裁剪轮廓填充。虽然这可能有点复杂，但一些有用的包可以帮助完成这些任务。

下面是代码。

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy.interpolate import griddata
import matplotlib.path as mpath

import alphashape
from descartes import PolygonPatch

Path = mpath.Path

data = pd.DataFrame({
    'time': [0, 1, 2, 3, 4, 0.1, 0.9, 2, 3.05, 4, 0, 1, 2.2, 3, 3.95],
    'force': [1, 2, 4, 9, 16, 0, 0, 0, 0, 0, -1, -2, -4, -9, -16]
})

Times, Forces = np.meshgrid(
    np.linspace(0, 4, 100),
    np.linspace(-16, 16, 100)
)

data['work'] = data['time'] * data['force']

interpolation = griddata(
    (data['time'], data['force']),
    data['work'],
    (Times, Forces),
    method= 'linear')

fig, ax = plt.subplots()
plt.xlim(-.3,4.3)
plt.ylim(-18,18)

#contour = ax.contourf( Times, Forces, interpolation)
contour = ax.tricontourf( data['time'], data['force'], data['work'])
ax.scatter(data['time'], data['force'])


x = data['time']
y =  data['force']

points = np.vstack([x, y]).T
alpha = 0.95 * alphashape.optimizealpha(points)
hull = alphashape.alphashape(points, alpha)
hull_pts = hull.exterior.coords.xy

ax.scatter(hull_pts[0], hull_pts[1], color='red')
ax.add_patch(PolygonPatch(hull, fill=False, color='red'))

plt.savefig("clip_before.png")
#make clip path
vertices = []
codes = []
xpts,ypts = hull_pts

## convert polygon to path for cliping contour fill
for ix,iy in zip(xpts,ypts):
    vertices.append((ix,iy))
codes += [Path.MOVETO]
codes += [Path.LINETO] * (len(xpts) -2)
codes += [Path.CLOSEPOLY]
clip = Path(vertices, codes)

for collection in contour.collections:
    collection.set_clip_path(clip,transform=ax.transData)
plt.savefig("clip_after.png")
plt.show()

这是输出数字。

剪裁前。

剪裁后。