0

I am trying to calculate the point of intersection (lat and lon in degrees) of two great circles that are each defined by two points on the circle. I have been trying to follow method outlined here. But the answer I get is incorrect, my code is below does anyone see where I went wrong?

import numpy as np
from numpy import cross
from math import cos, sin, atan2, asin, asinh

################################################
#### Intersection of two great circles.
# Points on great circle 1.
glat1 = 54.8639587
glon1 = -8.177818

glat2 = 52.65297082
glon2 = -10.78064876

# Points on great circle 2.
cglat1 = 51.5641564
cglon1 = -9.2754284

cglat2 = 53.35422063
cglon2 = -12.5767799

# 1. Put in polar coords.

x1 = cos(glat1) * sin(glon1)
y1 = cos(glat1) * cos(glon1)
z1 = sin(glat1)

x2 = cos(glat2) * sin(glon2)
y2 = cos(glat2) * cos(glon2)
z2 = sin(glat2)


cx1 = cos(cglat1) * sin(cglon1)
cy1 = cos(cglat1) * cos(cglon1)
cz1 = sin(cglat1)

cx2 = cos(cglat2) * sin(cglon2)
cy2 = cos(cglat2) * cos(cglon2)
cz2 = sin(cglat2)


# 2. Get normal to planes containing great circles.
#    It's the cross product of vector to each point from the origin.

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


# 3. Find line of intersection between two planes.
#    It is normal to the poles of each plane.

L = cross(N1, N2)


# 4. Find intersection points.

X1 = L / abs(L)
X2 = -X1


ilat = asin(X1[2]) * 180./np.pi
ilon = atan2(X1[1], X1[0]) * 180./np.pi

I should also mention this is on the Earth's surface (assuming a sphere).

2
  • 3
    I don't have time to check to see if this is the only problem, but it looks like your angles are in degrees while sin and cos expect radians.
    – DSM
    Feb 21, 2013 at 17:30
  • Ah, how did I not see that. Also X1 should be X1 = L / np.sqrt(L[0]**2 + L[1]**2 + L[2]**2)
    – Dave
    Feb 21, 2013 at 20:28

2 Answers 2

1

Solution from DSM in comments above, your angles are in degrees while sin and cos expect radians. Also the line

X1 = L / abs(L)

should be,

X1 = L / np.sqrt(L[0]**2 + L[1]**2 + L[2]**2) 
0

One more correction that needs to be done is to change cos/sin before "lon" in x and y dimensions:

x = cos(lat) * cos(lon)
y = cos(lat) * sin(lon)
z = sin(lat)

This is because original conversion from angle to spherical system is done using polar/azimuthal sphere angles and they are not the same as lat/lon angles (wiki it https://en.wikipedia.org/wiki/Spherical_coordinate_system).

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.