首页 > 脚本专栏 > python > python库geopy计算多组经纬度距离

python库geopy计算多组经纬度距离的实现方式

2023-08-16 10:06:32 作者：凡人叶枫

这篇文章主要介绍了python库geopy计算多组经纬度距离的实现方式,具有很好的参考价值,希望对大家有所帮助,如有错误或未考虑完全的地方,望不吝赐教

python库geopy计算多组经纬度距离

日常工作中有时会用到需要计算gnss定位模组的定位精确度，需要将被测设备和真值设备进行经纬度坐标之间的对比，由于经纬度坐标想要计算差值，需要涉及到坐标系的转换，计算方法比较复杂，geopy库很好的解决了这个问题，集成了大量的方法，可以做很多地理坐标相关的事情，其中就有计算两个坐标点之间距离的方法。

下面是我写的计算一系列坐标点之间的距离的python脚本，可以给大家提供参考。整理出两个设备输出的gps的utc时间以及经纬度，按照下面格式写到txt文本中，脚本读取txt文本，进行整秒的经纬度进行比较，输出某一时间点的定位误差，单位为m。

locationA.txt

1641541189000  38.018449    119.205511
1641541190000  38.018492    119.205554
1641541191000  38.018551    119.205591
1641541192000  38.018589    119.205625
1641541193000  38.018629    119.205657
1641541194000  38.018660    119.205693

locationB.txt

1641541189000  33.018449    119.205511
1641541190000  33.018492    114.205554
1641541191000  38.018551    116.205591
1641541192000  38.018589    116.205625
1641541193000  38.018629    119.205657
1641541194000  38.018660    119.205693

diff.txt

1641541189000	1.0
1641541190000	2.0
1641541191000	3.0
1641541192000	3.0
1641541193000	4.0
1641541194000	1.0

postion.py

# -*- coding:utf-8 -*-
#!/usr/bin/python3
from geopy.distance import distance
import time
dictpostionA = {}
dictpostionB = {}
dictdiff = {}
#转换时间戳 dt为字符串
def datetime_timestamp(dt):
    #中间过程，一般都需要将字符串转化为时间数组
    time.strptime(dt, '%Y/%m/%d %H:%M:%S')
    s = time.mktime(time.strptime(dt, '%Y/%m/%d %H:%M:%S'))
    return int(s)
def saveATempData(filename):
    with open(filename, "w") as f:
        for time, pos in dictpostionA.items():
            data = str(time) + "\t" + pos[0] + "\t" + pos[1] + "\n"
            #print(data)
            f.write(data)
#读取excel数据转成txt
def readAData(filename):
    with open(filename, "r") as f:
        for line in f.readlines():
            line = line.strip('\n')
            #print(line)
            poslist = line.split()
            #print(poslist[0])
            #print(poslist[1])
            #print(poslist[2])
            #print(poslist[3])
            #print(poslist[4])
            #print(poslist[5])
            dt = poslist[0]+' '+poslist[1]+':'+poslist[2]+':'+poslist[3]
            #print(dt)
            index = datetime_timestamp(dt) * 1000
            #print(index)
            sublist = poslist[4:6]
            #print(sublist)
            dictpostionA[index] = tuple(sublist)
def readBData(filename):
    with open(filename, "r") as f:
        for line in f.readlines():
            line = line.strip('\n')
            #print(line)
            poslist = line.split()
            #print(poslist[0])
            #print(poslist[1])
            #print(poslist[2])
            sublist = poslist[1:3]
            #print(sublist)
            index = int(poslist[0])
            dictpostionB[index] = tuple(sublist)
def processData():
    for timeA, posA in dictpostionA.items():
        for timeB, posB in dictpostionB.items():
            #print(timeA, timeB)
            if timeA == timeB:
                d = distance(posA, posB).m  # m是单位
                print(timeA, d)
                dictdiff[timeA] = d
                continue
def saveDiffData(filename):
    with open(filename, "w") as f:
        for time, diff in dictdiff.items():
            data = str(time) + "\t" + str(diff) + "\n"
            #print(data)
            f.write(data)
#读取数据
readAData("locationA.txt")
saveATempData("locationA_temp.txt")
readBData("locationB.txt")
#处理数据
processData()
#保存数据
saveDiffData("diff.txt")

python库Geopy用法，经纬度坐标转换、经纬度距离计算

Geopy库介绍

这里介绍一个Python 包 Geopy ，借助它也可以实现经纬度地理位置转换，

这款包之经纬度转换原理其实还是借助了第三方 API 平台，因为市面上提供经纬度转换第三方平台很多，为了方便， Geopy 把这些接口都分别封装在一个类中，借助 Geopy 模块来调用，支持的第三放平台如下

Geopy作为一个专注于地理处理包之外，除了能实现上面地理编码、逆地理编码功能之外，还有一个其它令我经验的功能， 提供两个经纬度坐标，计算他们在地球上的最短距离

下面将介绍一下 Geopy 的具体用法，

地理编码

使用 地理编码功能时，需要借助 Geopy 的 geocoders 模块，Geopy 把所有第三方API封装到 geocoders 中

这里选用 OpenStreetMap 平台上提供的 Nominatim 地理编码器，因为可以免费供我们使用，不需要申请 API ，但缺点是限流，限额，不能大规模频繁访问，否则会返回 403，429错误代码

from geopy.geocoders import Nominatim
geolocator=Nominatim()
location= geolocator.geocode("北京市海淀区西二旗北路")
print(location.address)
print(location.latitude,location.longitude)

结果如下

西二旗北路, 东北旺村, 海淀区, 北京市, 102208, 中国
40.056793 116.305811

逆地理编码

from geopy.geocoders import Nominatim
geolocator=Nominatim()
location= geolocator.reverse("40.056793 116.305811")
print(location.address)

结果如下

1#, 西二旗北路, 东北旺村, 海淀区, 北京市, 102208, 中国

结果看起来还不错，简单方便；但提醒一下，因为前面说过 Nominatim 模块是限额度的，不要频繁访问，否则会出现以下错误

根据经纬度计算距离

Geopy 最让我惊喜的是这个用法，提供两个经纬度坐标计算他们之间的距离，因为地球具体来说是椭圆，所以不能按照常规方法来计算，目前现有比较流行的几个模型有以下几个

              model             major (km)   minor (km)     flattening
ELLIPSOIDS = {'WGS-84':        (6378.137,    6356.7523142,  1 / 298.257223563),
              'GRS-80':        (6378.137,    6356.7523141,  1 / 298.257222101),
              'Airy (1830)':   (6377.563396, 6356.256909,   1 / 299.3249646),
              'Intl 1924':     (6378.388,    6356.911946,   1 / 297.0),
              'Clarke (1880)': (6378.249145, 6356.51486955, 1 / 293.465),
              'GRS-67':        (6378.1600,   6356.774719,   1 / 298.25),
              }

根据官方介绍，官网选择的是 WGS-84 模型，根据统计最终计算到的距离误差最高在0.5%左右；使用方法如下

from geopy import distance
newport_ri = (41.49008, -71.312796)
cleveland_oh = (41.499498, -81.695391)
print(distance.distance(newport_ri, cleveland_oh).miles)#最后以英里单位输出
#output
538.39044536
wellington = (-41.32, 174.81)
salamanca = (40.96, -5.50)
print(distance.distance(wellington, salamanca).km)# 以 km 作为单位输出
19959.6792674

总结

以上为个人经验，希望能给大家一个参考，也希望大家多多支持脚本之家。

python库geopy计算多组经纬度距离的实现方式