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(n1,n2):
    n = np.cross( n1, n2 )
    c = car2sph( *n )
    return c

p1  = np.array([ 156.8, 39.0 ])
p2  = np.array([ -85.2, 34.3 ])
n1  = sph2car( *p1 ) #p1 vector
n2  = sph2car( *p2 ) #p2 vector
p3  = np.array([ -56.8, -9.0 ])
p4  = np.array([ 115.2, 14.3 ])
n3  = sph2car( *p3 ) #p3 vector
n4  = sph2car( *p4 ) #p4 vector
# normal vector of plane 
np1  = np.cross( n1, n2 )
np2  = np.cross( n3, n4 )
pn01 = cross_greatcircles( npse, np01 )