蟒蛇 - 兩個大圓的交點（緯度/長）

Question

我想要計算兩個大圓之間的交點（lat/long，以度爲單位）（所以，我輸入兩個緯/長對，一個用於每個大圓代表其起點和終點，代碼應該返回經緯度）。我遵循其他帖子中的步驟，例如this和this。我在下面發佈我的代碼：它爲每個大圓圈對返回2個交點，因爲根據定義，任何兩個不同的大圓圈將相交於地球上的兩個地方。不幸的是，它返回了錯誤的結果（只通過在谷歌地球上繪製所有點來測試多次）。

import numpy as np 
import math 

# Define points in great circle 1 
p1_lat1 = 32.498520 
p1_long1 = -106.816846 
p1_lat2 = 38.199999 
p1_long2 = -102.371389 

# Define points in great circle 2 
p2_lat1 = 34.086771 
p2_long1 = -107.313379 
p2_lat2 = 34.910553 
p2_long2 = -98.711786 

# Convert points in great circle 1, degrees to radians 
p1_lat1_rad = ((math.pi * p1_lat1)/180.0) 
p1_long1_rad = ((math.pi * p1_long1)/180.0) 
p1_lat2_rad = ((math.pi * p1_lat2)/180.0) 
p1_long2_rad = ((math.pi * p1_long2)/180.0) 

# Convert points in great circle 2, degrees to radians 
p2_lat1_rad = ((math.pi * p2_lat1)/180.0) 
p2_long1_rad = ((math.pi * p2_long1)/180.0) 
p2_lat2_rad = ((math.pi * p2_lat2)/180.0) 
p2_long2_rad = ((math.pi * p2_long2)/180.0) 

# Put in polar coordinates 
x1 = math.cos(p1_lat1_rad) * math.sin(p1_long1_rad) 
y1 = math.cos(p1_lat1_rad) * math.cos(p1_long1_rad) 
z1 = math.sin(p1_lat1_rad) 
x2 = math.cos(p1_lat2_rad) * math.sin(p1_long2_rad) 
y2 = math.cos(p1_lat2_rad) * math.cos(p1_long2_rad) 
z2 = math.sin(p1_lat2_rad) 
cx1 = math.cos(p2_lat1_rad) * math.sin(p2_long1_rad) 
cy1 = math.cos(p2_lat1_rad) * math.cos(p2_long1_rad) 
cz1 = math.sin(p2_lat1_rad) 
cx2 = math.cos(p2_lat2_rad) * math.sin(p2_long2_rad) 
cy2 = math.cos(p2_lat2_rad) * math.cos(p2_long2_rad) 
cz2 = math.sin(p2_lat2_rad) 

# Get normal to planes containing great circles 
# np.cross product of vector to each point from the origin 
N1 = np.cross([x1, y1, z1], [x2, y2, z2]) 
N2 = np.cross([cx1, cy1, cz1], [cx2, cy2, cz2]) 

# Find line of intersection between two planes 
L = np.cross(N1, N2) 

# Find two intersection points 
X1 = L/np.sqrt(L[0]**2 + L[1]**2 + L[2]**2) 
X2 = -X1 
i_lat1 = math.asin(X1[2]) * 180./np.pi 
i_long1 = math.atan2(X1[1], X1[0]) * 180./np.pi 
i_lat2 = math.asin(X2[2]) * 180./np.pi 
i_long2 = math.atan2(X2[1], X2[0]) * 180./np.pi 

# Print results 
print i_lat1, i_long1, i_lat2, i_long2 

這將返回(34.314378256910636, -164.51906362183226, -34.314378256910636, 15.480936378167755)，但它應該圍繞(34.280552, -105.549375)的交叉點中的一個返回值。

關於我在做什麼錯的任何想法？非常感謝您的幫助！

---

編輯：供將來參考（如果某人運行到同樣的問題），我固定它由以下Ruei-Min Lin建議，所以更換線44和45與：

N1 = np.cross([y1, x1, z1], [y2, x2, z2]) 
N2 = np.cross([cy1, cx1, cz1], [cy2, cx2, cz2]) 

Answer 1

事實上，你可以從設置開始：

(p1_lat1, p1_long1) = (0, 0) 
(p1_lat2, p1_long2) = (90, 0) 
(p2_lat1, p2_long1) = (0, 0) 
(p2_lat2, p2_long2) = (0, 90) 

看看會發生什麼。 通過追蹤結果，您應該交換x1和y1，依此類推。

import numpy as np 
import math 

# Define points in great circle 1 
p1_lat1 = 32.498520 
p1_long1 = -106.816846 
p1_lat2 = 38.199999 
p1_long2 = -102.371389 

# Define points in great circle 2 
p2_lat1 = 34.086771 
p2_long1 = -107.313379 
p2_lat2 = 34.910553 
p2_long2 = -98.711786 

# Convert points in great circle 1, degrees to radians 
p1_lat1_rad = ((math.pi * p1_lat1)/180.0) 
p1_long1_rad = ((math.pi * p1_long1)/180.0) 
p1_lat2_rad = ((math.pi * p1_lat2)/180.0) 
p1_long2_rad = ((math.pi * p1_long2)/180.0) 

# Convert points in great circle 2, degrees to radians 
p2_lat1_rad = ((math.pi * p2_lat1)/180.0) 
p2_long1_rad = ((math.pi * p2_long1)/180.0) 
p2_lat2_rad = ((math.pi * p2_lat2)/180.0) 
p2_long2_rad = ((math.pi * p2_long2)/180.0) 

# Put in polar coordinates 
x1 = math.cos(p1_lat1_rad) * math.cos(p1_long1_rad) 
y1 = math.cos(p1_lat1_rad) * math.sin(p1_long1_rad) 
z1 = math.sin(p1_lat1_rad) 
x2 = math.cos(p1_lat2_rad) * math.cos(p1_long2_rad) 
y2 = math.cos(p1_lat2_rad) * math.sin(p1_long2_rad) 
z2 = math.sin(p1_lat2_rad) 
cx1 = math.cos(p2_lat1_rad) * math.cos(p2_long1_rad) 
cy1 = math.cos(p2_lat1_rad) * math.sin(p2_long1_rad) 
cz1 = math.sin(p2_lat1_rad) 
cx2 = math.cos(p2_lat2_rad) * math.cos(p2_long2_rad) 
cy2 = math.cos(p2_lat2_rad) * math.sin(p2_long2_rad) 
cz2 = math.sin(p2_lat2_rad) 

# Get normal to planes containing great circles 
# np.cross product of vector to each point from the origin 
N1 = np.cross([x1, y1, z1], [x2, y2, z2]) 
N2 = np.cross([cx1, cy1, cz1], [cx2, cy2, cz2]) 

# Find line of intersection between two planes 
L = np.cross(N1, N2) 

# Find two intersection points 
X1 = L/np.sqrt(L[0]**2 + L[1]**2 + L[2]**2) 
X2 = -X1 
i_lat1 = math.asin(X1[2]) * 180./np.pi 
i_long1 = math.atan2(X1[1], X1[0]) * 180./np.pi 
i_lat2 = math.asin(X2[2]) * 180./np.pi 
i_long2 = math.atan2(X2[1], X2[0]) * 180./np.pi 

# Print results 
print i_lat1, i_long1, i_lat2, i_long2