python已知經緯度，求大圓弧距離以及兩條大圓弧交點的經緯度(numpy實現)

python球面幾何的幾個實用的子函數(numpy實現)

如果不考慮地球的橢率，把地球近似為球體，下面的python程序可以快速計算出求球面上兩點之間的大圓弧距離，以及兩條大圓弧的交點。

原理：球坐標系轉為笛卡爾坐標系，然后進行矢量叉乘

import numpy as np

def sph2car(lon,lat):
    # input deg output 1
    i = np.deg2rad( lat )
    a = np.deg2rad( lon )
    z = np.sin( lat )
    x = np.cos( i ) * np.cos( a )
    y = np.cos( i ) * np.sin( a )
    return np.array([x,y,z])

def car2sph(x,y,z):
    # input 1 output deg
    l3  = np.sqrt( x*x+y*y+z*z )
    l2  = np.sqrt( x*x+y*y )
    lat = np.arccos( l2/l3 )
    lon = np.arctan( y / x )
    lat = np.rad2deg( lat  )
    lon = np.rad2deg( lon  )
    return np.array([lon,lat])

求兩點之間的大圓弧長度

def points2arc(p1,p2):
    # p1/2=(lon,lat), get their arc distance in deg
    n1 = sph2car( *p1 )
    n2 = sph2car( *p2 )
    co = np.dot( n1 , n2 ) # |n1|=1 and |n2|=1
    return np.rad2deg(np.arccos(co))

p1  = np.array([ 156.8, 39.0 ])
p2  = np.array([ -85.2, 34.3 ])
dis = points2arc(p1,p2)
print("The distance between ", p1, "and", p2, "is", dis, "deg") 

求兩條大圓弧之間的交點

def cross_greatcircles(p1,p2,p3,p4):
    # get the cross point of the line p1-p2 and line p3-p4
    # get the normal vector of the plane O-P1-P2
    n1  = sph2car( *p1 ) #p1 vector
    n2  = sph2car( *p2 ) #p2 vector
    n12 = np.cross( n1, n2 )
    # get the normal vector of the plane O-P3-P4
    n3  = sph2car( *p3 ) #p3 vector
    n4  = sph2car( *p4 ) #p4 vector
    n34  = np.cross( n3, n4 )
    # get their normal vector perpendicular to n12 and n34
    n = np.cross( n12, n34 )
    c = car2sph( *n )
    return c

p1  = np.array([ 156.8, 39.0 ])
p2  = np.array([ -85.2, 34.3 ])
p3  = np.array([ -56.8, -9.0 ])
p4  = np.array([ 115.2, 14.3 ])
pn  = cross_greatcircles(p1,p2,p3,p4)
print("Their cross point is", pn)