从纬度经度坐标以矢量化方式（无循环）构建距离矩阵

Question

我想提出一种更快的方法来创建所有经纬度对之间的距离矩阵。这个QA 解决了使用标准线性代数的矢量化方式，但没有经纬度坐标。

就我而言，这些经纬度是农场。这是我的 Python 代码，对于完整的数据集（4000 (lat, lon)'s）至少需要五分钟。有什么想法吗？

> def slowdistancematrix(df, distance_calc=True, sparse=False, dlim=100):
    """
    inputs: df

    returns:
    1.) distance between all farms in miles
    2.) distance^2

    """

    from scipy.spatial import distance_matrix
    from geopy.distance import geodesic

    unique_farms = pd.unique(df.pixel)
    df_unique = df.set_index('pixel')
    df_unique = df_unique[~df_unique.index.duplicated(keep='first')] # only keep unique index values
    distance = np.zeros((unique_farms.size,unique_farms.size))

    for i in range(unique_farms.size):
        lat_lon_i = df_unique.Latitude.iloc[i],df_unique.Longitude.iloc[i]
        for j in range(i):
            lat_lon_j = df_unique.Latitude.iloc[j],df_unique.Longitude.iloc[j]
            if distance_calc == True:
                distance[i,j] = geodesic(lat_lon_i, lat_lon_j).miles
                distance[j,i] = distance[i,j] # make use of symmetry

    return distance, np.power(distance, 2)

Answer 1

我的解决方案是this implementation的矢量化版本：

import numpy as np

def dist(v):
    v = np.radians(v)

    dlat = v[:, 0, np.newaxis] - v[:, 0]
    dlon = v[:, 1, np.newaxis] - v[:, 1]

    a = np.sin(dlat / 2.0) ** 2 + np.cos(v[:, 0]) * np.cos(v[:, 0]) * np.sin(dlon / 2.0) ** 2

    c = 2 * np.arcsin(np.sqrt(a))
    result = 3956 * c

    return result

但是，您需要使用属性 values 将数据框转换为 numpy 数组。例如：

df = pd.read_csv('some_csv_file.csv')
distances = dist(df[['lat', 'lng']].values)

Answer 2

这不是一个纯 python 解决方案，而是依赖于安装 r 与 geodist 包和 rpy2 接口：

import rpy2.robjects as ro
from rpy2.robjects.packages import importr
from rpy2.robjects import pandas2ri

from rpy2.robjects.conversion import localconverter


def pygeodist(pd_df):
    """
    pd_df must have columns 'x' and 'y' such that 'x' is the lng coordinate
    and 'y' is the lat coordinate
    """
    geodist=importr('geodist')
    with localconverter(ro.default_converter + pandas2ri.converter):
      return geodist.geodist(pd_df, measure = "geodesic")