计算GPS WGS_84 两点的距离

转载注明来源: 本文链接来自osnosn的博客，写于 2021-03-09.

python3 的包 geographiclib

文档【GeographicLib API】
参考【python 计算地球上两点距离和方位角（bearing）的包geographiclib】【Geodesic.WGS84.Inverse通过两点经纬度计算两点间的方位角】
Inverse() 函数偶尔会返回 nan

from geographiclib.geodesic import Geodesic
Geodesic.WGS84.Inverse(-41.32, 174.81, 40.96, -5.50)  #(lat1,lon1,lat2,lon2)

geopy 包

pip install geopy
文档【GeoPy】找 Great-circle distance
https://github.com/geopy/geopy

vincenty 包

pip install vincenty
文档【vincenty】

使用公式自己编写

低精度: javascript【计算GPS WGS_84 两点的距离】
高精度: javascript【计算GPS WGS_84 两点的距离更加细腻的算法】
R语言【已知经纬度计算两点间地理距离】
python3 低精度【根据经纬度坐标计算距离-python】【Python实现Haversine公式计算两点（经纬度坐标）距离】
c 高精度【GPS坐标（大地坐标）转高斯平面坐标，并计算 GPS 坐标（大地坐标）两点间的距离】
py3【python 计算两个点之间的距离和方向角，经纬度】
java【最全JAVA地球上两点间的距离算法！包含球体、椭球体、百度算法】
py3【经纬度计算的普通球面距离与精度更高的椭球距离】
js【获取经纬度之间的距离】【用JavaScript计算两点经纬度距离】

其他参考

网上搜索后，整理出 Haversine, Vincenty 等，四个算法的python3实现

Haversine最快。Vincenty最慢(耗时约多18-22%)。LL2Dist()和FlatternDist()速度一样(耗时约多8-10%)。
使用 geographiclib.geodesic.Geodesic.WGS84.Inverse() 耗时是 Haversine 的3.5-3.8倍。

import math
#from math import radians, cos, sin, asin, sqrt
# 这个算法误差比较大。5千km差几十km，3千km差十几km，1千km差几km，几百km差几百米，1km差十几米，
def haversine(lat1,lng1,lat2,lng2):
    '圆球计算,也叫Great-Circle Distance/大圆圆弧。计算两个 WGS84 经纬度点的距离'
    RR = 6378.137 #为地球赤道半径，单位:千米
    lng1, lat1, lng2, lat2 = map(math.radians, [float(lng1), float(lat1), float(lng2), float(lat2)]) # 经纬度转换成弧度
    dlon=lng2-lng1
    dlat=lat2-lat1
    a=math.sin(dlat/2)**2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon/2)**2
    distance=2*math.asin(math.sqrt(a))*RR*1000
    distance=round(distance,1) ##单位:米
    return distance

import math
#from math import radians,cos,sin,asin,sqrt,pi,atan,tan,atan2
# 这个算法和 geographiclib 算的几乎一样，相差<0.2米
def distVincenty(lat1,lon1,lat2,lon2):
    '精度更高的椭球计算,计算两个 WGS84 经纬度点的距离'
    a=6378137.0       #vincentyConstantA(WGS84) ##单位:米
    b=6356752.3142451 #vincentyConstantB(WGS84) ##单位:米
    f=1/298.257223563 #vincentyConstantF(WGS84)
    L = math.radians(lon2 - lon1)
    U1 = math.atan((1 - f) *math.tan(math.radians(lat1)))
    U2 = math.atan((1 - f) *math.tan(math.radians(lat2)))
    sinU1 =math.sin(U1)
    cosU1 =math.cos(U1)
    sinU2 =math.sin(U2)
    cosU2 =math.cos(U2)
    lambda1 = L
    lambdaP = 2 * math.pi
    iterLimit = 20

    sinLambda = 0.0
    cosLambda = 0.0
    sinSigma = 0.0
    cosSigma = 0.0
    sigma = 0.0
    alpha = 0.0
    cosSqAlpha = 0.0
    cos2SigmaM = 0.0
    C = 0.0
    while (abs(lambda1 - lambdaP) > 1e-12 and --iterLimit > 0) :
        sinLambda =math.sin(lambda1)
        cosLambda =math.cos(lambda1)
        sinSigma = math.sqrt((cosU2 * sinLambda) * (cosU2 * sinLambda) + (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) * (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda))
        if (sinSigma == 0) :
            return 0
        cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda
        sigma = math.atan2(sinSigma, cosSigma)
        alpha = math.asin(cosU1 * cosU2 * sinLambda / sinSigma)
        cosSqAlpha = math.cos(alpha) * math.cos(alpha)
        cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha
        C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha))
        lambdaP = lambda1
        lambda1 = L + (1 - C) * f * math.sin(alpha)* (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)))

    if iterLimit == 0 :
        return 0.0

    uSq = cosSqAlpha * (a * a - b * b) / (b * b)
    A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)))
    B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)))
    deltaSigma = B * sinSigma * (cos2SigmaM + B / 4 * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) - B / 6 * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2SigmaM * cos2SigmaM)))
    s = b * A * (sigma - deltaSigma)
    d = s  ##单位:米
    return d

import math
# 这个算法,几百km误差几十米,几千km误差几百米，不知道叫什么算法
# 两个经纬度之间的距离,椭球
def LL2Dist(Lat1,Lng1,Lat2,Lng2):
    ra = 6378137.0        # radius of equator: meter
    rb = 6356752.3142451  # radius of polar: meter
    flatten = (ra - rb) / ra  # Partial rate of the earth
    if Lat1==Lat2 and Lng1==Lng2:
        return 0
    # change angle to radians
    radLatA = math.radians(Lat1)
    radLonA = math.radians(Lng1)
    radLatB = math.radians(Lat2)
    radLonB = math.radians(Lng2)

    pA = math.atan(rb / ra * math.tan(radLatA))
    pB = math.atan(rb / ra * math.tan(radLatB))
    x = math.acos(math.sin(pA) * math.sin(pB) + math.cos(pA) * math.cos(pB) * math.cos(radLonA - radLonB))
    c1 = (math.sin(x) - x) * (math.sin(pA) + math.sin(pB))**2 / math.cos(x / 2)**2
    c2 = (math.sin(x) + x) * (math.sin(pA) - math.sin(pB))**2 / math.sin(x / 2)**2
    dr = flatten / 8 * (c1 - c2)
    distance = ra * (x + dr)  ##单位:米
    return distance

import math
# 这个算法,几百km误差几十米,几千km误差几百米，比LL2Dist()误差大十几米，也不知道叫什么算法
# 两个经纬度之间的距离,椭球
def FlatternDist(lat1,lng1,lat2,lng2):
    ra = 6378137.0        # radius of equator: meter
    #rb = 6356752.3142451  # radius of polar: meter
    flatten=1/298.257223563 #vincentyConstantF(WGS84)
    #flatten = (ra - rb) / ra  # Partial rate of the earth
    if lat1==lat2 and lng1==lng2:
        return 0
    # change angle to radians
    f = math.radians((lat1+lat2)/2)
    g = math.radians((lat1-lat2)/2)
    l = math.radians((lng1-lng2)/2)
    sf = math.sin(f)
    sg = math.sin(g)
    sl = math.sin(l)
    sg = sg * sg
    sl = sl * sl
    sf = sf * sf

    s = sg * (1 - sl)+(1 - sf) * sl
    c = (1 - sg) * (1 - sl) + sf * sl

    w = math.atan(math.sqrt(s / c))
    r = math.sqrt(s * c) / w
    d = 2 * w * ra
    h1 = (3 * r - 1) / 2 / c
    h2 = (3 * r + 1) / 2 / s

    ##单位:米
    return d * (1 + flatten * (h1 * sf * (1 - sg) - h2 * (1 - sf) * sg))

转载注明来源: 本文链接来自osnosn的博客.

相关阅读:
简单了解一下：var 、let、const
C# FlagAttriute 的小妙招
项目经验面试题
linux面试题详解
jvm面试题详解
数据库面试详解
微服务框架面试题
框架面试题（maven、ZooKeeper、Dubbo、Nginx、Redis、Lucene、Solr、ActiveMQ、JMS
设计模式面试题详解
WEB方面面试题详解

原文地址：https://www.cnblogs.com/osnosn/p/14505778.html