获取两个geopandas数据框几何点之间的距离

Question

我是第一次使用空间数据。我必须比较两个具有纬度和经度细节的数据框。我已将两者都转换为 GeoPandas 数据帧。

import pandas as pd
from pandas import DataFrame
import geopandas as gpd
from neighbors import nearest_neighbor


df = pd.DataFrame([[1973,22.525158,88.330775],[1976,72.85136,19.10840],[898,91.78523,26.15012]],columns=['id', 'lat', 'long'])
gdf1 = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.long,df.lat))

df2 = pd.DataFrame([['06c979eaa59f',29.873870,76.965620],['19aedbb2e743',20.087574,76.180045],['5060a3931a43',31.289770,75.572340]],columns=['id','lat','lon']) 
gdf2 = gpd.GeoDataFrame(df2, geometry=gpd.points_from_xy(df2.lon,df2.lat))

我的 DF1 有 100 万行，df2 有大约 7000 行。我正在尝试为 DF1 中的每条记录从 DF2 获取最近的邻居。

两种方法我都试过了。两者都运行得非常快并且结果可行。但是，它们并不准确。

方法一：

Please check this link

在此页面中，我使用了 sklearn.neighbors 中的最近邻方法。这 returns 结果以米为单位。但是，当我手动检查两个数据帧的经纬度之间的距离时，我总是找到最近的邻居 returns 1/4 的距离。

比如上述方法返回的距离是125米，那么google和https://www.geodatasource.com/distance-calculatorreturns都是500米左右的距离。距离的差异一直在返回结果的4倍左右波动。

方法二：

在第二种方法中，我遵循了 gis.stackexchange.com 中给出的代码。

https://gis.stackexchange.com/questions/222315/geopandas-find-nearest-point-in-other-dataframe

import itertools
from operator import itemgetter

import geopandas as gpd
import numpy as np
import pandas as pd

from scipy.spatial import cKDTree
from shapely.geometry import Point, LineString

df = pd.DataFrame([[1973,22.525158,88.330775],[1976,72.85136,19.10840],[898,91.78523,26.15012]],columns=['id', 'lat', 'long'])
gdf1 = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.long,df.lat))

df2 = pd.DataFrame([['06c979eaa59f',29.873870,76.965620],['19aedbb2e743',20.087574,76.180045],['5060a3931a43',31.289770,75.572340]],columns=['id','lat','lon']) 
gdf2 = gpd.GeoDataFrame(df2, geometry=gpd.points_from_xy(df2.lon,df2.lat))

这里，我用自己的数据框替换了gpd1和gpd2。

def ckdnearest(gdfA, gdfB, gdfB_cols=['id']):
    # resetting the index of gdfA and gdfB here.
    gdfA = gdfA.reset_index(drop=True)
    gdfB = gdfB.reset_index(drop=True)
    A = np.concatenate(
        [np.array(geom.coords) for geom in gdfA.geometry.to_list()])
    B = [np.array(geom.coords) for geom in gdfB.geometry.to_list()]
    B_ix = tuple(itertools.chain.from_iterable(
        [itertools.repeat(i, x) for i, x in enumerate(list(map(len, B)))]))
    B = np.concatenate(B)
    ckd_tree = cKDTree(B)
    dist, idx = ckd_tree.query(A, k=1)
    idx = itemgetter(*idx)(B_ix)
    gdf = pd.concat(
        [gdfA, gdfB.loc[idx, gdfB_cols].reset_index(drop=True),
         pd.Series(dist, name='dist')], axis=1)
    return gdf

c = ckdnearest(gdf1, gdf2)

以上运行速度非常快，returns结果。但是返回的距离值至少比我得到的低100倍。

乘数：107.655914

在上面的excel图中，第一列是python返回的结果，而第二列是上面给出的同一个网站返回的结果。虽然这些近似结果让我开始了，但我想要准确的结果。如何比较上面给出的两个数据框并获得 DF1 中每一行的最准确的最近距离。

Answer 1

在处理空间数据时，您应该知道点坐标是从球体投影到平面上的。在墨卡托投影中，经纬度点之间的距离以度为单位，而不是米。并且转换取决于点的纬度，因为赤道上的 1 度将小于高纬度上的 1 度。

您可以查看此讨论以了解此问题的可能解决方案： https://gis.stackexchange.com/questions/293310/how-to-use-geoseries-distance-to-get-the-right-answer

举个例子，一种可能是您将地理数据框转换为覆盖您所在区域的 UTM 投影。例如，比利时与 UTM 区域 31N EPSG:32631 相交。墨卡托投影有一个 epsg 代码 EPSG:4326。要转换 GeoDataFrame/GeoSeries，您需要在创建它时提供 CRS：

s = gpd.GeoSeries(points, crs=4326)

其中 points 是 shapely.geometry.Point

的列表

然后转换为给定的 UTM：

s_utm = s.to_crs(epsg=32631)

现在，您计算 s_utm 中各点之间的距离将以米为单位。

但是您需要确保您的点确实落在给定的 UTM 区域中，否则结果将不准确。我链接的答案表明其他方法也可能有效并且可以应用于整个点集合。

您也可以尝试转换为应该保留距离的 EPSG 32663（WGS 84/世界等距圆柱）。

另一个选项可以使用 geopy，它允许使用 geopy.geodesic.distance

计算测地线距离

获取两个geopandas数据框几何点之间的距离

Getting the distance between two geopandas data frames geometry points

python

geospatial

pandas