查找2个轮廓的交点坐标

Question

使用"Finding intersection of two contour plots in Python" 作为指导，我收到以下错误消息（代码如下）：

<ipython-input-98-993ccd512742> in <module>
      9 
     10 
---> 11 intersection_example = findIntersection(c1,c2)
     12 
<ipython-input-97-995cbb9fd0d0> in findIntersection(contour1, contour2)
      2 
      3 def findIntersection(contour1,contour2):
----> 4   p1 = contour1.collections[0].get_paths()[0]
      5   v1 = p1.vertices
      6 
IndexError: list index out of range

下面的第一个代码示例给了我一个没有错误的 3D 等高线图：

import matplotlib.pyplot as plt
import numpy as np
    
fig = plt.figure()
fig.set_size_inches(18.5, 10.5)   #, forward=True)
ax = plt.axes(projection='3d')
x = np.linspace(0, 21, 20)
y = np.linspace(0, 21, 20)
X, Y = np.meshgrid(x, y)

ax.contour3D(X, Y, ((X - 5) * (Y - 5) - 25), 29, cmap='winter')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z');
#fig.add_subplot(ax.contour3D(X, Y, Z2, 70, cmap='winter')) #binary'))
ax.contour3D(X, Y, (X**2 + Y**2 - 400), 29, cmap='autumn')
ax.set_title('Contour3D: Using meshgrid X Y ')

以上产生：

下一个示例是使用contour（而不是contour3D）并导致错误的代码段：

IndexError: 列表索引超出范围

这也是findIntersection使用未定义参数调用时产生的错误。

from shapely import geometry

def findIntersection(contour1,contour2):
    p1 = contour1.collections[0].get_paths()[0]
    v1 = p1.vertices
    
    p2 = contour2.collections[0].get_paths()[0]
    v2 = p2.vertices
    
    poly1 = geometry.LineString(v1)
    poly2 = geometry.LineString(v2)
    
    intersection = poly1.intersection(poly2)
    
    return intersection

figtst2 = plt.figure()
figtst2.set_size_inches(18.5, 10.5)   #, forward=True)
ax2 = plt.axes(projection='3d')
c1 = ax2.contour(X,Y,((X - 5) * (Y - 5) - 25),1,colors='green', linewidths=3)
c2 = ax2.contour(X,Y,(X**2 + Y**2 - 400),1,colors='orange', linewidths=3)

ax2.set_title('Contour Using meshgrid X Y & looking for intersections')


# Error is generated on the next line
intersection_example = findIntersection(c1,c2)

# where coordinates can be accessed by

intersection_example.x ##get x points
intersection_example.y ##get y points
list(intersection_example.coords)  ## get in [x,y] formatting

绘制ax2 产生：

注意：如果我使用线性空间 x 和 y 而不是网格 X 和 Y 我得到：

TypeError: 输入 z 必须是二维数组。

Answer 1

我开发了一个有点软糖的答案，所以如果有办法在 3D 中完成这一切，请告诉我。 软糖是做一个 2D 绘图来获取我的交叉点坐标，然后我做一个 3D 绘图并使用我在 2D 绘图中找到的内容（请参阅 Jupyter notebook 代码）

from shapely import geometry
fig2 = plt.figure()
fig2.set_size_inches(18.5, 10.5)   #, forward=True)
#plt.axes(projection='3d')  # 3D gives no vertices "list index out of range"

tau = np.arange(0,23,1)

x,y= np.meshgrid(tau,tau)
cs = plt.contour(x, y, np.array((x - 5) * (y - 5) - 25), [1],colors='k')
cs2 = plt.contour(x, y, np.array((x**2 + y**2 - 400)), [1],colors='k')  #x**2 + y**2 - 400



from shapely.geometry import LineString
v1 = cs.collections[0].get_paths()[0].vertices
v2 = cs2.collections[0].get_paths()[0].vertices

ls1 = LineString(v1)
ls2 = LineString(v2)
points = ls1.intersection(ls2)
print('points', np.array(points))

这给出了以下图像

3D 码是；

def fhyp(x, y):
    return ((x - 5) * (y - 5) - 25)

x = np.linspace(0, 21, 20)
y = np.linspace(0, 21, 20)

X, Y = np.meshgrid(x, y)
Z = fhyp(X, Y)

def fcirc(x, y):
    return (x**2 + y**2 - 400)

x = np.linspace(0, 21, 20)
y = np.linspace(0, 21, 20)

X, Y = np.meshgrid(x, y)
Z2 = fcirc(X, Y)
fig = plt.figure()
fig.set_size_inches(18.5, 10.5)   #, forward=True)
ax = plt.axes(projection='3d')
x = np.linspace(0, 21, 20)
y = np.linspace(0, 21, 20)
X, Y = np.meshgrid(x, y)

ax.contour3D(X, Y, Z, 1, cmap='winter') #4th parm was 29 to give 29 contour lines
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z');
#
ax.contour3D(X, Y, Z2, 1, cmap='autumn') 
ax.set_title('Contour3D: Using meshgrid X Y ')


from scipy import optimize

def f(p):
    x, y = p
    e1 = ((x - 5) * (y - 5) - 25)
    e2 = (x**2 + y**2 - 400)
return e1, e2

x2, y2 = optimize.fsolve(f, (5, 5))

zval = ((x2 - 5) * (y2 - 5) - 25)

print(x2,y2, zval)
vals = [x2,y2,zval]
result = np.array(vals)
print(result)

plt.plot([result[0]],[result[1]],[result[2]], "rx")

for i in range(5):
    x = [18.80452816482531,18.80452816482531]
    y = [6.810999963323182, 6.810999963323182]
    z = [-400,400]
plt.plot(x,y,z,'k--',alpha=0.8, linewidth=0.5)