根据经纬度计算范围_亿级经纬度距离计算88.73秒,秒杀VBA!
计算经纬度的代码网上一搜一大把,通常是单点距离的计算,无法实现批量计算,本文将利用pandas实现亿级经纬度距离代码的实现。最短距离计算建议参考下文,mapinfo能够很好的实现。MAPINFO 最小站间距统计本文将实现两张表的任意点之间100、200、300、500、800、1000米范围内的距离计算。首先导入需要使用的包1importpandasaspd2importnump...
·
计算经纬度的代码网上一搜一大把,通常是单点距离的计算,无法实现批量计算,本文将利用pandas实现亿级经纬度距离代码的实现。最短距离计算建议参考下文,mapinfo能够很好的实现。
MAPINFO 最小站间距统计本文将实现两张表的任意点之间100、200、300、500、800、1000米范围内的距离计算。首先导入需要使用的包
偶然间想起了之前自己将
MAPINFO 最小站间距统计本文将实现两张表的任意点之间100、200、300、500、800、1000米范围内的距离计算。首先导入需要使用的包
1import pandas as pd2import numpy as np3from math import radians, cos, sin, asin, sqrt, ceil4import math5import 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-lng15 dlat=lat2-lat16 a=sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**27 distance=2*asin(sqrt(a))*6371*1000 # 地球平均半径,6371km8 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)
然后根据经度和纬度差值进行过滤(经纬度差值大于某个值,距离大于某个值,参见参数表
1 x = m[abs(m.lon-m.lon2) 2 n = x[abs(x.lat-x.lat2)
得到下图表格:
geodistance
自定义函数,此处使用pandas
内置模块apply
(比使用for
循环要高效很多)。
1 nn = n.copy()
2 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 # 完成一个文件写入之后,文件编号增加112 filecount += 1
详情可以参考如下文章。Python工具开发实践-csv文件分割将文件分割之后,我们便可以循环处理分片文件与目标文件,将得到的结果合并到一个空的Dataframe
里st_time)))
1distance = pd.DataFrame(columns=('name','lon','lat','name2', 'lon2', 'lat2', 'distance'))
2for i in range(1, filecount):
3 df_temp = pd.read_csv(file_name[:-4] + '_' + str(i) + '.csv')
4 m = pd.concat([pd.concat([df_temp]*len(df2)).sort_index().reset_index(drop=True),
5 pd.concat([df2]*len(df_temp)).reset_index(drop=True)], 1)
6 # 避免链式赋值
7 x = m[abs(m.lon-m.lon2) 8 n = x[abs(x.lat-x.lat2) 9 nn = n.copy()
10 nn['distance'] = nn.apply(lambda ser: geodistance(ser['lon'], ser['lat'], ser['lon2'], ser['lat2']), axis=1)
11 distance = distance.append(nn[nn.distance <= minx_mile])
12distance.to_csv('D:/python/geo/distance_result.csv')
使用测试数据测算,经纬度距离亿次计算量耗时约88.73
秒,秒杀VBA。

DAMO开发者矩阵,由阿里巴巴达摩院和中国互联网协会联合发起,致力于探讨最前沿的技术趋势与应用成果,搭建高质量的交流与分享平台,推动技术创新与产业应用链接,围绕“人工智能与新型计算”构建开放共享的开发者生态。
更多推荐
所有评论(0)