python求交点坐标_Python计算两个大圆的交点

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).

解决方案

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)