在Python中测量两个经纬度点之间的距离（以千为单位）

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次迭代。我担心第二个数据集中的记录数量会增加。

有更好的方法可以节省时间吗？

- Apricot

---

@0x6d64 for循环迭代。谢谢您的建议...我会检查一下。 - Apricot

for循环的哪一部分？或者你有机会找出距离计算中哪一部分是昂贵的吗？ - 0x6d64

你可不可以简单地去掉第二个循环，采用矢量化操作呢？ - Mercury

你能否分享数据（仅经纬度部分），以便有兴趣解决此问题的人可以用于测试？我的做法是分两步进行：（1）使用不太准确的公式进行粗略筛选（2）使用准确的公式进行实际距离计算。 - Niko Pasanen

顺便说一下，您正在使用的公式是Haversine距离，不太准确。对于400米的距离，它可能会产生1米左右的误差，但我想在您的情况下这不是问题，因为您只需要粗略过滤点。 - fdermishin

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)
        

- fdermishin

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 方法（空间连接）也很快。我怀疑这与算法包括边界框和/或排序坐标有关... 我使用这种方法取得了一些不错的结果。

---

警告：

我将数据集投影到EPSG 3857中，该投影具有全球性且具有笛卡尔坐标（以米为单位）。关于任何“真实”项目，您必须仔细选择投影（即在您的区域选择最好的“度量系统友好型”投影），以避免缓冲区的任何扭曲...

- tgrandje

---

使用专门的库会更好。可能有更聪明的人花费了很多时间来确保它的正确性和速度！ - 0x6d64

使用sjoin的方式很有趣！然而，6371 * 1000 * 400的半径肯定是错误的，因为它具有平方米的维度，但缓冲区接受距离作为参数。此外，使用epsg 4326似乎会扭曲距离，应该使用度量投影，正如您所提到的那样。因此，如果坐标在全球范围内分布，则这种解决方案似乎不太适合。 - fdermishin

关于重新投影的一些线索在这里 - 这完全取决于数据集的真实位置以及该地区是否存在合适的投影方式；-) - tgrandje

我在查找其他主题的资源时偶然发现了这个。该 Github 仓库似乎在进行任何缓冲之前将投影转换为 epsg 3857（确实是以米为单位）。我刚刚编辑了答案，以此投影作为默认投影。 - tgrandje

抱歉，我打错了...（已经更正：您应该使用gpd.sjoin(gdf1，gdf2)进行调用） - tgrandje

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的那个。

- HrkBrkkl

Answer 4

你可以使用BallTree和HaversineDistance度量来处理。首先，使用第一个表中的坐标构建树，然后从该树查询第二个表中的坐标。

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]])

- fdermishin