python庫geographiclib計算經緯度方位角距離

Geographiclib庫的使用

geographiclib庫的dochttps://geographiclib.sourceforge.io/html/python/examples.html#initializing有時進不去，帶來很大困擾，於是把它搬運到博客里。

安裝

pip install geographiclib

初始化

from geographiclib.geodesic import Geodesic
geod = Geodesic.WGS84  # define the WGS84 ellipsoid
geod = Geodesic(6378388, 1/297.0) # altanatively custom the ellipsoid

已知兩點經緯度，求方位角距離圓弧長

g = geod.Inverse( lat1, lon1, lat2, lon2 )
print("Distance=",g['s12'], "m and Great circle Arc=", g['a12'],"km") ###官網此處可能有誤，已更正
'''
注意g['s12']的單位應該是m不是km，此處不同於官網！詳見評論區
'''
print("From 1 to 2, azimuth=", g['azi1'],"deg and from 2 to 1, azimuth=", g['azi2'], "deg")

已知一點經緯度方位角距離，求另一點經緯度

#unit: azimuth (N=0deg), distance (m)
g = geod.Direct( lat, lon, azimuth, distance )
print("The new point is", g['lat2'], "N,", g['lon2'], "E.")

已知兩點，計算大圓弧路徑

l = geod.InverseLine(40.1, 116.6, 37.6, -122.4)
ds = 1000e3; n = int(math.ceil(l.s13 / ds))
for i in range(n + 1):
  if i == 0:
    print("distance latitude longitude azimuth")
  s = min(ds * i, l.s13)
  g = l.Position(s, Geodesic.STANDARD | Geodesic.LONG_UNROLL)
  print("{:.0f} {:.5f} {:.5f} {:.5f}".format(g['s12'], g['lat2'], g['lon2'], g['azi2']))

輸出結果如下

distance latitude longitude azimuth
0 40.10000 116.60000 42.91642
1000000 46.37321 125.44903 48.99365
2000000 51.78786 136.40751 57.29433
3000000 55.92437 149.93825 68.24573
4000000 58.27452 165.90776 81.68242
5000000 58.43499 183.03167 96.29014
6000000 56.37430 199.26948 109.99924
7000000 52.45769 213.17327 121.33210
8000000 47.19436 224.47209 129.98619
9000000 41.02145 233.58294 136.34359
9513998 37.60000 237.60000 138.89027

已知不規則圖形的經緯度，求面積

p = geod.Polygon()
antarctica = [...
  [-63.1, -58], [-72.9, -74], [-71.9,-102], [-74.9,-102], [-74.3,-131],...
  [-77.5,-163], [-77.4, 163], [-71.7, 172], [-65.9, 140], [-65.7, 113],...
  [-66.6,  88], [-66.9,  59], [-69.8,  25], [-70.0,  -4], [-71.0, -14],...
  [-77.3, -33], [-77.9, -46], [-74.7, -61]...
]
for pnt in antarctica:
  p.AddPoint(pnt[0], pnt[1])
num, perim, area = p.Compute()
print("Perimeter/area of Antarctica are {:.3f} m / {:.1f} m^2".format(perim, area))