Python实践 | 亿级经纬度距离计算工具V2
Python实践 | 亿级经纬度距离计算工具V2
阅读文本大概需要 8 分钟。
计算经纬度的代码网上一搜一大把，通常是单点距离的计算，无法实现批量计算，本文将利用pandas实现亿级经纬度距离代码的实现。 最短距离计算建议参考下文，mapinfo能够很好的实现。 MAPINFO 最小站间距统计
本文将实现两张表的任意点之间100、200、300、500、800、1000米范围内的距离计算。 首先导入需要使用的包
1import pandas as pd
2import numpy as np
3from math import radians, cos, sin, asin, sqrt, ceil
4import math
5import time经纬度计算自定义函数
1def geodistance(lng1,lat1,lng2,lat2):
2 lng1, lat1, lng2, lat2 = map(radians, [float(lng1), float(lat1), float(lng2), float(lat2)])
3 # 经纬度转换成弧度
4 dlon=lng2-lng1
5 dlat=lat2-lat1
6 a=sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
7 distance=2*asin(sqrt(a))*6371*1000 # 地球平均半径，6371km
8 distance=round(distance,0)
9 return distance实现不同范围内的距离计算，例如100、200、300、500、800、1000，适合做成一张参数表。
由于地球是球形，不同纬度下，同一经度差值对应的距离不同，纬度相同且纬度越大时，同一经度对应的距离越小，中国经纬度跨度约为73°33′E至 135°05′E；纬度范围：3°51′N至53°33′N，此处为了计算最大经度差值，我们选取纬度为54.0；不同经度下，同一纬度差异对应的距离相同
不同经纬度差异对应最小距离表格如下：
pandas分别导入源表和目标表，两个表关联得到原点与目标点的所有配对
1file_name = r'D:\python\geo\sTable.csv'
2df1=pd.read_csv(file_name)
3file_name2 = r'D:\python\geo\tTable.csv'
4df2=pd.read_csv(file_name2)
5m = pd.concat([pd.concat([df1]*len(df2)).sort_index().reset_index(drop=True),
6 pd.concat([df2]*len(df1)).reset_index(drop=True) ], 1)然后根据经度和纬度差值进行过滤（经纬度差值大于某个值，距离大于某个值，参见参数表
x = m[abs(m.lon-m.lon2) < diff_lon]
n = x[abs(x.lat-x.lat2) < diff_lat]得到下图表格：
然后针对每一行的4个参数应用geodistance自定义函数，此处使用pandas内置模块apply（比使用for循环要高效很多）。
nn = n.copy()
nn['distance'] = nn.apply(lambda ser: geodistance(ser['lon'], ser['lat'], ser['lon2'], ser['lat2']), axis=1)根据经纬度差值判断距离是一个大致的范围，我们选取纬度值54.0获取了最大的经度差值，随着纬度减小，此时计算的距离会大于该阈值，所以要对初次计算结果进行过滤，得出满足阈值的条目：
1distance = distance.append(nn[nn.distance <= minx_mile])经过调试，发现该方法计算量上限基本为1000万，当计算量大于1000万怎么办？
偶然间想起了之前自己将csv文件分割的文章，当计算量大于1000万，我们对原表进行分割，分割个数就是计算量/10000000，不能整除时，需要先上取整，多分割一个文件
1pieces = ceil(count_a * count_b / 10000000) # 计算量上限为1000万分割数目有了，文件分片大小也就有了
1linesPerFile = ceil(count_a / pieces)+1文件分割代码：
1filecount = 1
2# 以0为起点，文件行数为终点，分片大小为间隔，循环遍历文件，每次遍历行数即为分片大小，而不是每行遍历一次，处理效率极高，但是比较吃内存
3for i in range(0, len(csv_file), linesPerFile):
4 # 打开目标文件准备写入，不存在则创建
5 with open(file_name[:-4] + '_' + str(filecount) + '.csv', 'w+') as f:
6 # 判断是否为第一个文件，不是的话需要先写入标题行
7 if filecount > 1:
8 f.write(csv_file[0])
9 # 批量写入i至i+分片大小的多行数据，效率极高
10 f.writelines(csv_file[i:i+linesPerFile])
11 # 完成一个文件写入之后，文件编号增加1
12 filecount += 1详情可以参考如下文章。
将文件分割之后，我们便可以循环处理分片文件与目标文件，将得到的结果合并到一个空的Dataframe里st_time)))
distance = pd.DataFrame(columns=('name','lon','lat','name2', 'lon2', 'lat2', 'distance'))
for i in range(1, filecount):
df_temp = pd.read_csv(file_name[:-4] + '_' + str(i) + '.csv')
m = pd.concat([pd.concat([df_temp]*len(df2)).sort_index().reset_index(drop=True),
pd.concat([df2]*len(df_temp)).reset_index(drop=True)], 1)
# 避免链式赋值
x = m[abs(m.lon-m.lon2) < diff_lon]
n = x[abs(x.lat-x.lat2) < diff_lat]
nn = n.copy()
nn['distance'] = nn.apply(lambda ser: geodistance(ser['lon'], ser['lat'], ser['lon2'], ser['lat2']), axis=1)
distance = distance.append(nn[nn.distance <= minx_mile])
distance.to_csv('D:/python/geo/distance_result.csv')使用测试数据测算，经纬度距离亿次计算量耗时约88.73秒，秒杀VBA。
公众号回复distance获取完整代码