以下代码用python2运行
#!/usr/bin/python
#coding=utf-8
from __future__ import print_function
import sys
reload(sys)
sys.setdefaultencoding("utf-8")
import mathdef haversine_distance(lon1, lat1, lon2, lat2):# 将纬度和经度从度转换为弧度lat1_rad = math.radians(lat1)lon1_rad = math.radians(lon1)lat2_rad = math.radians(lat2)lon2_rad = math.radians(lon2)# 计算差值dlon = lon2_rad - lon1_raddlat = lat2_rad - lat1_rad# 应用 Haversine 公式a = math.sin(dlat / 2)**2 + math.cos(lat1_rad) * math.cos(lat2_rad) * math.sin(dlon / 2)**2c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))distance = 6371 * c * 1000 # 地球平均半径约为 6371 千米return distancedef vincenty(lon1, lat1, lon2, lat2):# 定义椭球体参数a = 6378137.0 # 长半轴（单位：米）f = 1 / 298.257223563 # 扁率b = (1 - f) * a # 短半轴（单位：米）# 将经纬度从度数转换为弧度lon1, lat1, lon2, lat2 = map(math.radians, [lon1, lat1, lon2, lat2])# 差值L = lon2 - lon1U1 = math.atan((1 - f) * math.tan(lat1))U2 = math.atan((1 - f) * math.tan(lat2))sinU1, cosU1 = math.sin(U1), math.cos(U1)sinU2, cosU2 = math.sin(U2), math.cos(U2)# 初始设定lamb = Liter_limit = 100 # 最大迭代次数for _ in range(iter_limit):sin_lamb = math.sin(lamb)cos_lamb = math.cos(lamb)sin_sigma = math.sqrt((cosU2 * sin_lamb) ** 2 + (cosU1 * sinU2 - sinU1 * cosU2 * cos_lamb) ** 2)if sin_sigma == 0:return 0 # 重合点cos_sigma = sinU1 * sinU2 + cosU1 * cosU2 * cos_lambsigma = math.atan2(sin_sigma, cos_sigma)sin_alpha = cosU1 * cosU2 * sin_lamb / sin_sigmacos2_alpha = 1 - sin_alpha ** 2cos2_sigma_m = cos_sigma - 2 * sinU1 * sinU2 / cos2_alphaif math.isnan(cos2_sigma_m):cos2_sigma_m = 0 # 赤道情况C = f / 16 * cos2_alpha * (4 + f * (4 - 3 * cos2_alpha))lamb_prev = lamblamb = L + (1 - C) * f * sin_alpha * (sigma + C * sin_sigma * (cos2_sigma_m + C * cos_sigma * (-1 + 2 * cos2_sigma_m ** 2)))if abs(lamb - lamb_prev) < 1e-12:breakelse:raise ValueError("Vincenty formula failed to converge")u2 = cos2_alpha * (a ** 2 - b ** 2) / (b ** 2)A = 1 + u2 / 16384 * (4096 + u2 * (-768 + u2 * (320 - 175 * u2)))B = u2 / 1024 * (256 + u2 * (-128 + u2 * (74 - 47 * u2)))delta_sigma = B * sin_sigma * (cos2_sigma_m + B / 4 * (cos_sigma * (-1 + 2 * cos2_sigma_m ** 2) -B / 6 * cos2_sigma_m * (-3 + 4 * sin_sigma ** 2) *(-3 + 4 * cos2_sigma_m ** 2)))s = b * A * (sigma - delta_sigma)return s # 距离，单位：米# 示例：计算北京和上海之间的距离
beijing_lon, beijing_lat = 116.407396, 39.904200
shanghai_lon, shanghai_lat = 121.473701, 31.230416distance = vincenty(beijing_lon, beijing_lat, shanghai_lon, shanghai_lat)
print("北京和上海之间的距离是 {:.2f} 米".format(distance))
haversine算法与vincenty的差距很小，对于500多米的近距离，误差在1米多。
将haversine用sql语句表达结果为：
select 6371*1000 * acos(sin(radians(a.latitude)) * sin(radians(b.latitude)) + cos(radians(a.latitude)) * cos(radians( b.latitude)) * cos(radians(b.longitude) - radians(a.longitude))) AS distance