经纬度求实际距离

计算距离时有3种地球模型可供选择.

• 平面
• 球体
• 椭球体

import math

#勾股定理求取 误差最大
def CalcPythagorean(Lat_A,Lng_A,Lat_B,Lng_B):
    # lat 维度
    # lng 经度
    EARTH_RADIUS=6378137
    x1=Lat_A
    y1=Lng_A
    x2=Lat_B
    y2=Lng_B
    x_dist=math.radians(x1-x2)  #角度转成弧度
    y_dist=math.radians(y1-y2)  #角度转成弧度
    dist_sq=x_dist**2+y_dist**2
    dist_rad=math.sqrt(dist_sq)
    distance=dist_rad*EARTH_RADIUS
    return distance

#半正矢公式 误差在1km以内
def calcDistance(Lat_A,Lng_A,Lat_B,Lng_B):
    # lat 维度
    # lng 经度
    EARTH_RADIUS=6378137
    x1=Lat_A
    y1=Lng_A
    x2=Lat_B
    y2=Lng_B
    x_dist=math.radians(x1-x2)
    y_dist=math.radians(y1-y2)
    y1_rad=math.radians(y1)
    y2_rad=math.radians(y2)
    a=math.sin(y_dist/2)**2+(math.sin(x_dist/2)**2)*math.cos(y1_rad)*math.cos(y2_rad)
    c=2*math.asin(math.sqrt(a))
    distance=c*EARTH_RADIUS
    return distance

#Vincenty公式 精度最高 远远小于1km
def VincentyDistance(Lat_A,Lng_A,Lat_B,Lng_B):
    # lat 维度
    # lng 经度
    x1=Lat_A
    y1=Lng_A
    x2=Lat_B
    y2=Lng_B
    f=1/298.257222101
    a = 6378137  #半长轴
    b=abs((f*a)-a) #半短轴
    L=math.radians(x2-x1)
    U1=math.atan((1-f)*math.tan(math.radians(y1)))
    U2=math.atan((1-f)*math.tan(math.radians(y2)))
    sinU1=math.sin(U1)
    cosU1=math.cos(U1)
    sinU2=math.sin(U2)
    cosU2=math.cos(U2)
    lam=L
    for i in range(100):
        sinLam=math.sin(lam)
        cosLam=math.cos(lam)
        sinSigma=math.sqrt((cosU2*sinLam)**2+(cosU1*sinU2-sinU1*cosU2*cosLam)**2)
        if(sinSigma==0):
            distance=0
            break
        cosSigma=sinU1*sinU2+cosU1*cosU2*cosLam
        sigma=math.atan2(sinSigma,cosSigma)
        sinAlpha=cosU1*cosU2*sinLam/sinSigma
        cosSqAlpha=1-sinAlpha**2
        cos2SigmaM=cosSigma-2*sinU1*sinU2/cosSqAlpha
        if math.isnan(cos2SigmaM):
            cos2SigmaM=0 #赤道线
        C=f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha))
        LP=lam
        lam=L+(1-C)*f*sinAlpha*(sigma+C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)))
        if not abs(lam-LP)>1e-12:
            break
    uSq=cosSqAlpha*(a**2-b**2)/b**2
    A=1+uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)))
    B=uSq/1024*(256+uSq*(-128+uSq*(74-47*uSq)))
    delteSigma=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-delteSigma)
    distance=s
    return distance