python - 两个大圆的交点 (lat/long)
python - intersection point of two great circles (lat/long)
我正在尝试计算两个大圆之间的交点（lat/long，以度为单位）（因此，我输入了两个 lat/long 对，每个大圆一个，代表它的起点和终点，并且代码应该 return lat/long 交叉点）。我按照其他帖子中的步骤操作，例如this and this。我在下面发布我的代码：每个大圆对 returns 2 个交点，因为根据定义，任何两个不同的大圆都会在地球上的两个地方相交。不幸的是，它是 returning 错误的结果（通过绘制 google 地球上的所有点测试了几次）。
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
这个 returns (34.314378256910636, -164.51906362183226, -34.314378256910636, 15.480936378167755)，但它应该是 returning 值 (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])
其实可以从设置开始：
(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
我正在尝试计算两个大圆之间的交点（lat/long，以度为单位）（因此，我输入了两个 lat/long 对，每个大圆一个，代表它的起点和终点，并且代码应该 return lat/long 交叉点）。我按照其他帖子中的步骤操作，例如this and this。我在下面发布我的代码：每个大圆对 returns 2 个交点，因为根据定义，任何两个不同的大圆都会在地球上的两个地方相交。不幸的是，它是 returning 错误的结果（通过绘制 google 地球上的所有点测试了几次）。
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
这个 returns (34.314378256910636, -164.51906362183226, -34.314378256910636, 15.480936378167755)，但它应该是 returning 值 (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])
其实可以从设置开始：
(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