在 Python 中测量 1000 中两个经纬度点之间的距离

Question

我有两个数据框。

df1 有 580 条唯一记录 - 带有纬度和经度信息

df2 有 490000 条唯一记录 - 带有经纬度信息

我正在尝试获取 - 在这 580 个位置中，490000 个位置的 400 米半径范围内有多少个位置。

我正在使用以下代码并且它正在运行。

from __future__ import print_function
from config import conn
from pandas import DataFrame
import pandas as pd
import math

def distance(origin, destination):
    lat1, lon1 = origin
    lat2, lon2 = destination
    radius = 6371 *1000# km

    dlat = math.radians(lat2-lat1)
    dlon = math.radians(lon2-lon1)
    a = math.sin(dlat/2) * math.sin(dlat/2) + math.cos(math.radians(lat1)) \
        * math.cos(math.radians(lat2)) * math.sin(dlon/2) * math.sin(dlon/2)
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))
    d = radius * c
    return d

def convertTuple(tup): 
    str =  ''.join(tup) 
    return str


df1 = pd.read_csv("/home/ubuntu/maid80.csv")
df2 = pd.read_csv("/home/ubuntu/iodr.csv")
ll = []
for index,rows in df2.iterrows():
        lat1 = rows['latitude']
        lon1 = rows['longitude']
        for i,r in df1.iterrows():
                k = distance((lat1,lon1),(r['latitude'],r['longitude']))
                if (k <= 400):
                        ll.append(rows['id'])
#                       print(ll)
        print(index)
        myset = set(ll)
        print(myset)

我运行这不是我的笔记本电脑，完成所有 580 次迭代需要 2 个多小时。我担心第二个数据集中的记录数量会膨胀。

有没有更好的方法，可以节省时间。

Answer 1

按纬度对两个数据框进行排序。如果它们的纬度差异很大，这将允许不计算点对之间的距离。在最好的情况下，您可以获得 580 倍的加速。

想法是遍历 df1 的行，并为该数组的每一行找到第二个数组的左右索引，这些索引的纬度离该行不远

df1.sort_values(by='latitude')
df2.sort_values(by='latitude')
n1 = df1.shape[0]
n2 = df2.shape[0]
left = 0
right = 0
threshold = 400
lat_threshold = threshold / radius # latitude difference that corresponds to 400 m
for i in range(n1):
    row1 = df1.iloc[[i]]
    lat1 = row1['latitude']
    lon1 = row1['longitude']
    while left < n2 and df2.iloc[[left]]['latitude'] < lat1 - lat_threshold:
        left += 1
    while right < n2 and df2.iloc[[right]]['latitude'] < lat1 + lat_threshold:
        right += 1
    for j in range(left, right):
        row2 = df2.iloc[[j]]
        lat2 = row2['latitude']
        lon2 = row2['longitude']
        k = distance((lat1, lon1), (lat2, lon2))
        if (k <= threshold):
            ll.append(row2)

Answer 2

您可以使用 geopandas 进行尝试：

import geopandas as gpd
import pandas as pd
import pyproj

df1 = pd.read_csv("/home/ubuntu/maid80.csv")
df2 = pd.read_csv("/home/ubuntu/iodr.csv")

gdf1 = gpd.GeoDataFrame(df1, geometry=gpd.points_from_xy(df1['longitude'], df1['latitude']), crs=pyproj.CRS.from_epsg(4326))
gdf2 = gpd.GeoDataFrame(df2, geometry=gpd.points_from_xy(df2['longitude'], df2['latitude']), crs=pyproj.CRS.from_epsg(4326))

radius = 400
for gdf in [gdf1, gdf2]:
  gdf.to_crs(pyproj.CRS.from_epsg(3857), inplace=True)

gdf1['geometry'] = gdf1['geometry'].buffer(radius)
gdf2['IS_WITHIN_400M'] = 1

gdf = gpd.sjoin(gdf1, gdf2['geometry'], how='left')
print(gdf[gdf.IS_WITHIN_400M_right==1].head())

一些解释：

Geopandas 将允许您使用 GeoDataFrame，您可以在其上使用半径（非常快）“缓冲”您的点。 points_from_xy 函数也非常快，可以让您高效地构建这些对象。

sjoin 方法（代表空间连接）也很快。我怀疑这与包含边界框 and/or 排序坐标的算法有关...我使用此方法取得了一些不错的结果。

警告：

我将数据集投影到 EPSG 3857 中，它是全局的 AND 具有笛卡尔坐标（以米为单位）。对于任何“真实”项目，您必须仔细选择投影（即选择您所在地区最好的“公制友好”投影）以避免缓冲区失真...

Answer 3

您只能对距离函数使用 numpy 函数并将其矢量化。那应该快很多：

from __future__ import print_function

import pandas as pd
import math

import numpy as np


def distance(origin: pd.DataFrame, lat2, lon2):
'''Measure distance not for a pair but for the whole dataframa against one point'''
    lat1 = origin['latitude']
    lon1 = origin['longitude']
    radius = 6371 * 1000  # km
    dlat = np.radians(lat2 - lat1)
    dlon = np.radians(lon2 - lon1)
    a = np.sin(dlat / 2) * np.sin(dlat / 2) + np.cos(np.radians(lat1)) \
        * np.cos(np.radians(lat2)) * np.sin(dlon / 2) * np.sin(dlon / 2)
    c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))
    d = radius * c
    return d


def main():
    df1 = pd.read_csv("/home/ubuntu/maid80.csv")
    df2 = pd.read_csv("/home/ubuntu/iodr.csv")
    ll = []
    for index, row in df2.iterrows():
        #because you can test the whole dataframe gainst one point you can remove    one loop.
        mask= distance(df1,row['latitude'],row['longitude'])<400.0
        ll.extend(df1.index[mask].to_list()) #only add points to the list where the distance is <400

    
    myset = set(ll)
    print(myset)

也许您必须切换数据帧。不知道你要从哪一个收集id。

Answer 4

您可以使用 BallTree with HaversineDistance 指标。首先使用第一个 table 的坐标构建树，然后从该树

的第二个 table 查询坐标

from sklearn.neighbors import BallTree, DistanceMetric

radius = 6371 * 1000
max_distance = 400 / radius

# ensure that format of array is [latitude, longitude]
rows1 = np.deg2rad(df1[['latitude', 'longitude']].to_numpy())
rows2 = np.deg2rad(df2[['latitude', 'longitude']].to_numpy())

# haversine metric accepts latitude and longitude only in radians and returns distance
# on unit sphere
tree = BallTree(rows1, metric=DistanceMetric.get_metric('haversine'))

count = tree.query_radius(rows2, r=max_distance, count_only=True)
print(df2['id'].iloc[np.nonzero(count)[0]])

在 Python 中测量 1000 中两个经纬度点之间的距离

Measuring distance between two lat long points in 1000's in Python

python

latitude-longitude