Python實踐 | 億級經緯度距離計算工具V2
- 2020 年 2 月 26 日
- 筆記
閱讀文本大概需要 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:pythongeosTable.csv' 2df1=pd.read_csv(file_name) 3file_name2 = r'D:pythongeotTable.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獲取完整代碼
