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