为何Cartopy与Basemap的逆函数计算地球表面两点距离结果不同？

你遇到的结果差异，核心原因是两个库对经纬度参数的顺序要求完全相反，结合你给出的坐标定义，我们来修正代码并解释问题：

先明确你的坐标定义 #

你提到c0和c1是[纬度, 经度]格式的数组：

import numpy as np
c0 = np.array([77.343750, 22.593726])  # [纬度, 经度]
c1 = np.array([86.945801, 23.684774])  # [纬度, 经度]

两个库的参数顺序差异 #

1. Basemap（基于pyproj）的Geod.inv()
   这个方法要求参数顺序是 (经度0, 纬度0, 经度1, 纬度1)，也就是必须先传经度，再传纬度。
   你原代码里的传参逻辑是对的，但实际得到的结果明显偏离正确值，大概率是写代码时误把参数顺序搞反了（比如写成了c0[0], c0[1], c1[0], c1[1]）。修正后的正确代码：

   import mpl_toolkits.basemap.pyproj as pyproj
   k = pyproj.Geod(ellps="WGS84")
   # 严格按照：经度0, 纬度0, 经度1, 纬度1 的顺序传参
   _, _, distance_basemap = k.inv(c0[1], c0[0], c1[1], c1[0])
   distance_basemap = distance_basemap / 1000
   print(distance_basemap)  # 现在会输出和Cartopy一致的 ~990.6公里
   

2. Cartopy的Geodesic.inverse()
   这个方法接受的是**[纬度, 经度]**格式的点数组，刚好和你的坐标定义匹配，所以你的原代码是完全正确的，得到的结果是准确的：

   import cartopy.geodesic as gd
   k = gd.Geodesic()  # 默认使用WGS84椭球，无需额外指定
   distance_cartopy = k.inverse(c0, c1).base[0,0]/1000
   print(distance_cartopy)  # 输出：990.6094719605074
   

验证正确结果 #

用WGS84椭球计算这两个北极附近点的球面距离，正确值确实在990.6公里左右，所以Cartopy的结果是准确的，Basemap的问题完全源于参数顺序的错误。

内容的提问来源于stack exchange，提问作者Anveshan Lal