我试图找到一条直线(红色虚线)与以红色突出显示的轮廓线的交点(见图)。我在第二个图中使用 .get_paths 将所述等高线与其他(第二个图)隔离。
我查看了轮廓相交问题 How to find all the intersection points between two contour-set in an efficient way ，并尝试将其用作基础，但未能重现任何有用的内容。
http://postimg.org/image/hz01fouvn/
http://postimg.org/image/m6utofwb7/
有没有人有任何想法？
重新创建绘图的相关功能，
#for contour
def p_0(num,t) :
esc_p = np.sum((((-1)**n)*(np.exp(t)**n)*((math.factorial(n)*((n+1)**0.5))**-1)) for n in range(1,num,1))
return esc_p+1
tau = np.arange(-2,3,0.1)
r=[]
p1 = p_0(51,tau)
p2 = p_0(51,tau)
for i in p1:
temp_r=i/p2
r.append(temp_r)
x,y= np.meshgrid(tau,tau)
cs = plt.contour(x, y, np.log(r),50,colors='k')
whichContour =20
pa = CS.collections[whichContour].get_paths()[0]
v = pa.vertices
xx = v[:, 0]
yy = v[:, 1]
plt.plot(xx, yy, 'r-', label='Crossing contour')
#straight line
p=0.75
logp = (np.log(p*np.exp(tau)))
plt.plot(tau,logp)
目前的尝试，
import matplotlib
import numpy as np
import matplotlib.cm as cm
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import math
def intercepting_line() :
matplotlib.rcParams['xtick.direction'] = 'out'
matplotlib.rcParams['ytick.direction'] = 'out'
#fake data
delta = 0.025
x = np.arange(-3.0, 3.0, delta)
y = np.arange(-2.0, 2.0, delta)
X, Y = np.meshgrid(x, y)
Z1 = mlab.bivariate_normal(X, Y, 1.0, 1.0, 0.0, 0.0)
Z2 = mlab.bivariate_normal(X, Y, 1.5, 0.5, 1, 1)
Z = 10.0 * (Z2 - Z1)
#plot
cs = plt.contour(X,Y,Z)
whichContour = 2 # change this to find the right contour lines
#get the vertices to calculate an intercept with a line
p = cs.collections[whichContour].get_paths()[0]
#see: http://matplotlib.org/api/path_api.html#module-matplotlib.path
v = p.vertices
xx = v[:, 0]
yy = v[:, 1]
#this shows the innermost ring now
plt.plot(xx, yy, 'r--', label='inner ring')
#fake line
x = np.arange(-2, 3.0, 0.1)
y=lambda x,m:(m*x)
y=y(x,0.9)
lineMesh = np.meshgrid(x,y)
plt.plot(x,y,'r' ,label='line')
#get the intercepts, two in this case
x, y = find_intersections(v, lineMesh[1])
print x
print y
#plot the intercepting points
plt.plot(x[0], y[0], 'bo', label='first intercept')
#plt.plot(x[1], y[1], 'rs', label='second intercept')
plt.legend(shadow=True, fancybox=True, numpoints=1, loc='best')
plt.show()
#now we need to calculate the intercept of the vertices and whatever line
#this is pseudo code but works in case of two intercepting contour vertices
def find_intersections(A, B):
# min, max and all for arrays
amin = lambda x1, x2: np.where(x1<x2, x1, x2)
amax = lambda x1, x2: np.where(x1>x2, x1, x2)
aall = lambda abools: np.dstack(abools).all(axis=2)
slope = lambda line: (lambda d: d[:,1]/d[:,0])(np.diff(line, axis=0))
x11, x21 = np.meshgrid(A[:-1, 0], B[:-1, 0])
x12, x22 = np.meshgrid(A[1:, 0], B[1:, 0])
y11, y21 = np.meshgrid(A[:-1, 1], B[:-1, 1])
y12, y22 = np.meshgrid(A[1:, 1], B[1:, 1])
m1, m2 = np.meshgrid(slope(A), slope(B))
m1inv, m2inv = 1/m1, 1/m2
yi = (m1*(x21-x11-m2inv*y21) + y11)/(1 - m1*m2inv)
xi = (yi - y21)*m2inv + x21
xconds = (amin(x11, x12) < xi, xi <= amax(x11, x12),
amin(x21, x22) < xi, xi <= amax(x21, x22) )
yconds = (amin(y11, y12) < yi, yi <= amax(y11, y12),
amin(y21, y22) < yi, yi <= amax(y21, y22) )
return xi[aall(xconds)], yi[aall(yconds)]
目前它找到了相交点，但只在线条一致的地方，我在这里找不到解决方案的主要原因是我不明白原作者在这里的思路，
yi = (m1*(x21-x11-m2inv*y21) + y11)/(1 - m1*m2inv)
xi = (yi - y21)*m2inv + x21
最佳答案
使用 shapely 可以找到交点，而不是使用该点作为 fsolve() 的初始猜测值来找到真正的解决方案:
#for contour
def p_0(num,t) :
esc_p = np.sum((((-1)**n)*(np.exp(t)**n)*((math.factorial(n)*((n+1)**0.5))**-1)) for n in range(1,num,1))
return esc_p+1
tau = np.arange(-2,3,0.1)
x,y= np.meshgrid(tau,tau)
cs = plt.contour(x, y, np.log(p_0(51, y)/p_0(51, x)),[0.2],colors='k')
p=0.75
logp = (np.log(p*np.exp(tau)))
plt.plot(tau,logp)
from shapely.geometry import LineString
v1 = cs.collections[0].get_paths()[0].vertices
ls1 = LineString(v1)
ls2 = LineString(np.c_[tau, logp])
points = ls1.intersection(ls2)
x, y = points.x, points.y
from scipy import optimize
def f(p):
x, y = p
e1 = np.log(0.75*np.exp(x)) - y
e2 = np.log(p_0(51, y)/p_0(51, x)) - 0.2
return e1, e2
x2, y2 = optimize.fsolve(f, (x, y))
plt.plot(x, y, "ro")
plt.plot(x2, y2, "gx")
print x, y
print x2, y2
这是输出:
0.273616328952 -0.0140657435002
0.275317387697 -0.0123646847549
和情节:
关于python - 寻找直线和轮廓的交点，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/22232311/