我有一个可能很简单的问题，这让我已经安静了一会儿。有没有一种简单的方法可以在 python matplotlib 中返回两个绘制（非分析）数据集的交集？
为了详细说明，我有这样的事情：
x=[1.4,2.1,3,5.9,8,9,23]
y=[2.3,3.1,1,3.9,8,9,11]
x1=[1,2,3,4,6,8,9]
y1=[4,12,7,1,6.3,8.5,12]
plot(x1,y1,'k-',x,y,'b-')
本例中的数据完全是任意的。我现在想知道是否有一个我一直缺少的简单内置函数，它可以返回两个图之间的精确交集。
希望我说清楚了，也希望我没有遗漏一些非常明显的东西......
【问题讨论】：
标签： python matplotlib intersection
我们可以使用scipy.interpolate.PiecewisePolynomial 创建由您的分段线性数据定义的函数。
p1=interpolate.PiecewisePolynomial(x1,y1[:,np.newaxis])
p2=interpolate.PiecewisePolynomial(x2,y2[:,np.newaxis])
然后我们可以取这两个函数的不同，
def pdiff(x):
return p1(x)-p2(x)
并使用optimize.fsolve 查找pdiff 的根：
import scipy.interpolate as interpolate
import scipy.optimize as optimize
import numpy as np
x1=np.array([1.4,2.1,3,5.9,8,9,23])
y1=np.array([2.3,3.1,1,3.9,8,9,11])
x2=np.array([1,2,3,4,6,8,9])
y2=np.array([4,12,7,1,6.3,8.5,12])
p1=interpolate.PiecewisePolynomial(x1,y1[:,np.newaxis])
p2=interpolate.PiecewisePolynomial(x2,y2[:,np.newaxis])
def pdiff(x):
return p1(x)-p2(x)
xs=np.r_[x1,x2]
xs.sort()
x_min=xs.min()
x_max=xs.max()
x_mid=xs[:-1]+np.diff(xs)/2
roots=set()
for val in x_mid:
root,infodict,ier,mesg = optimize.fsolve(pdiff,val,full_output=True)
# ier==1 indicates a root has been found
if ier==1 and x_min<root<x_max:
roots.add(root[0])
roots=list(roots)
print(np.column_stack((roots,p1(roots),p2(roots))))
产量
[[ 3.85714286 1.85714286 1.85714286]
[ 4.60606061 2.60606061 2.60606061]]
第一列是 x 值，第二列是在 x 处计算的第一个 PiecewisePolynomial 的 y 值，第三列是第二个 PiecewisePolynomial 的 y 值。
【讨论】：
参数解决方案
如果序列 {x1,y1} 和 {x2,y2} 定义任意 (x,y) 曲线，而不是 y(x) 曲线，我们需要一种参数化方法来寻找交点。由于如何做到这一点并不完全清楚，而且@unutbu 的解决方案在 SciPy 中使用了一个已失效的插值器，我认为重新审视这个问题可能会很有用。
import numpy as np
from numpy.linalg import norm
from scipy.optimize import fsolve
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
x1_array = np.array([1,2,3,4,6,8,9])
y1_array = np.array([4,12,7,1,6.3,8.5,12])
x2_array = np.array([1.4,2.1,3,5.9,8,9,23])
y2_array = np.array([2.3,3.1,1,3.9,8,9,11])
s1_array = np.linspace(0,1,num=len(x1_array))
s2_array = np.linspace(0,1,num=len(x2_array))
# Arguments given to interp1d:
# - extrapolate: to make sure we don't get a fatal value error when fsolve searches
# beyond the bounds of [0,1]
# - copy: use refs to the arrays
# - assume_sorted: because s_array ('x') increases monotonically across [0,1]
kwargs_ = dict(fill_value='extrapolate', copy=False, assume_sorted=True)
x1_interp = interp1d(s1_array,x1_array, **kwargs_)
y1_interp = interp1d(s1_array,y1_array, **kwargs_)
x2_interp = interp1d(s2_array,x2_array, **kwargs_)
y2_interp = interp1d(s2_array,y2_array, **kwargs_)
xydiff_lambda = lambda s12: (np.abs(x1_interp(s12[0])-x2_interp(s12[1])),
np.abs(y1_interp(s12[0])-y2_interp(s12[1])))
s12_intercept, _, ier, mesg \
= fsolve(xydiff_lambda, [0.5, 0.3], full_output=True)
xy1_intercept = x1_interp(s12_intercept[0]),y1_interp(s12_intercept[0])
xy2_intercept = x2_interp(s12_intercept[1]),y2_interp(s12_intercept[1])
plt.plot(x1_interp(s1_array),y1_interp(s1_array),'b.', ls='-', label='x1 data')
plt.plot(x2_interp(s2_array),y2_interp(s2_array),'r.', ls='-', label='x2 data')
if s12_intercept[0]>0 and s12_intercept[0]<1:
plt.plot(*xy1_intercept,'bo', ms=12, label='x1 intercept')
plt.plot(*xy2_intercept,'ro', ms=8, label='x2 intercept')
plt.legend()
print('intercept @ s1={}, s2={}\n'.format(s12_intercept[0],s12_intercept[1]),
'intercept @ xy1={}\n'.format(np.array(xy1_intercept)),
'intercept @ xy2={}\n'.format(np.array(xy2_intercept)),
'fsolve apparent success? {}: "{}"\n'.format(ier==1,mesg,),
'is intercept really good? {}\n'.format(s12_intercept[0]>=0 and s12_intercept[0]<=1
and s12_intercept[1]>=0 and s12_intercept[1]<=1
and np.isclose(0,norm(xydiff_lambda(s12_intercept)))) )
返回，对于这个初始猜测的特定选择 [0.5,0.3]：
intercept @ s1=0.4761904761904762, s2=0.3825944170771757
intercept @ xy1=[3.85714286 1.85714286]
intercept @ xy2=[3.85714286 1.85714286]
fsolve apparent success? True: "The solution converged."
is intercept really good? True
这个方法只找到一个交集：我们需要迭代几个初始猜测（就像@unutbu 的代码所做的那样），检查它们的真实性，并使用np.close 消除重复项。请注意，fsolve 可能在返回值ier 中错误地指示成功检测到交叉点，这就是为什么在这里进行额外检查的原因。
【讨论】：