首页 > 编程语言 >经纬度距离计算Vincenty's formulae算法

经纬度距离计算Vincenty's formulae算法

时间:2022-12-24 17:44:39浏览次数:64  
标签:lam Vincenty 经纬度 cosSigma sinSigma uSq formulae cos2SigmaM math

维基百科公式
原始论文

import math
from geopy.distance import geodesic


# 计算两点之间的椭球面距离
# 使用Vincent算法,使用WGS84参考椭球体。
# 输入两个点的经纬度坐标,输出距离,单位为米
def distance(lat1, lon1, lat2, lon2):
    # WGS84 参考椭球体参数 Ellipsoid Parameters
    a = 6378137.000  # 长半轴 semi-major axis
    b = 6356752.3142  # 短半轴 semi-minor axis
    f = 1 / 298.257223563  # 扁率 flattening
    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, cosU1, sinU2, cosU2 = math.sin(U1), math.cos(U1), math.sin(U2), math.cos(U2)
    lam = L

    for i in range(100):
        sinLam, cosLam = math.sin(lam), math.cos(lam)
        sinSigma = math.sqrt((cosU2 * sinLam) ** 2 +
                             (cosU1 * sinU2 - sinU1 * cosU2 * cosLam) ** 2)
        if sinSigma == 0:
            return 0
        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  # equatorial line
        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)))
    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)
    return s


e = [63, 12]
f = [63, -12]
dis = geodesic(e, f).km
print(dis)
dispaper = distance(63, 12, 63, -12) / 1000
print(dispaper)

标签:lam,Vincenty,经纬度,cosSigma,sinSigma,uSq,formulae,cos2SigmaM,math
From: https://www.cnblogs.com/liuyechang/p/17003091.html

相关文章

  • python之调用高德、百度api解析经纬度地址
    调用高德#高德地图根据经纬度反查地址,每天只能调用5000次defgaode_excute_single_query(coordStrings,currentkey='你自己的api-key'):#1-将coordList中的经纬......
  • 高德地图地址转化为经纬度(php)
    高德地图地址转化为经纬度(php)高德地图开放平台:https://developer.amap.com/1、注册一个高德开放平台链接地址:https://lbs.amap.com/dev/id/choose2、创建一个应用,......
  • Qt之自定义输入框(度分秒、经纬度、格式化显示)
    相关资料:http://www.manongjc.com/detail/15-grpefyhtwdpbehh.html  Qt自定义文本输入框实现支持输入度分秒和度两种格式(简易无限制输入)PS:重要的文件我用粗体标注......
  • SQL Server根据地图坐标经纬度计算距离
    实战1-亲测,返回结果是米(m)如果嫌麻烦直接跳转至:“参考方案一”注意:该方法仅支持SQLServer2008和该版本以上的数据库 1、准备-工具百度坐标拾取器:​​http://api.map.ba......
  • 百度开发者根据地名查询经纬度api
    说明工作中,项目经理让我全国的省市区的经纬度进行统计到数据库中,在网上也没有找到好的方式,有些网站也只有一个个查,后来,我想到通过爬虫的形式调用api进行获取,然后找到了百......
  • selenium高级用法:获取经纬度
    导言获取经纬度的方法有很多,通过调用某地图API,模块geopy。但这两种方式都有一定的缺点,调用某地图API访问次数有限,使用模块geopy虽然次数不受限制,但是这个模块只能精确到镇,如......
  • 经纬度坐标转换为工程坐标
    1.绪论在施工和工程测量时,经常需要将GPS坐标转换为工程中所使用的坐标。在Bing上检索到的类似问题,基本描述为两个坐标系的转换,但实际上并非如此。本文将详细解释转换过......
  • 国内经纬度随机数据
    国内经纬度随机数据随机数据:[121.48,31.22,121.24,31.4,121.48,31.41,121.7,31.19,121.76,31.05,121.46,30.92,121.24,31,121.16,30.89,121.1,31.15,121.4,31.73,102.7......
  • 移动端H5获取用户当前经纬度,地理位置信息
    移动端浏览器有自带的获取用户经纬度的方法,但是位置会有偏差。如果只是要个大概,就当我没说。这里我用的是腾讯的api,直接https引用。在页面上加上 <script charset="ut......
  • 由经纬度通过接口获取到省市区
    exportconstgetAddress=function(wx){returnnewPromise((resolve,reject)=>{let_this=thiswx.getSetting({success:(res)=>{......