I'm trying to calculate the intersection points (lat/long, in degrees) between two great circle (so, I input two lat/long pairs, one for each great circle, representing its start and end points, and the code is supposed to return the lat/long intersections). I followed the steps in other posts, such as this and this. I'm posting my code below: it returns 2 intersection points for each great circle pair, since by definition any two distinct great circles will intersect in two places on the earth. Unfortunately, it's returning wrong results (just tested several times by plotting all points on google earth).
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
This returns (34.314378256910636, -164.51906362183226, -34.314378256910636, 15.480936378167755), but it should be returning values around (34.280552, -105.549375) for one of the intersection points.
Any ideas on what I'm doing wrong? Thank you very much for your help!
EDIT: For future reference (in case somebody runs into the same problem), I fixed it by following Ruei-Min Lin advice, so replacing lines 44 and 45 with:
N1 = np.cross([y1, x1, z1], [y2, x2, z2])
N2 = np.cross([cy1, cx1, cz1], [cy2, cx2, cz2])
In fact, you can start by setting:
(p1_lat1, p1_long1) = (0, 0)
(p1_lat2, p1_long2) = (90, 0)
(p2_lat1, p2_long1) = (0, 0)
(p2_lat2, p2_long2) = (0, 90)
to see what happens. By tracing the result, you should swap x1 and y1, and so on.
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