我正在尝试计算两个大圆之间的交点（纬度/经度，以度为单位）（因此，我输入了两个纬度/经度对，每个大圆一个，代表它的起点和终点，并且代码应该返回纬度/经度交叉点）。我按照其他帖子中的步骤进行操作，例如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])