python - 获取两个 geopandas 数据框几何点之间的距离
问题描述
我第一次使用空间数据。我必须比较两个具有经纬度和经度详细信息的数据框。我已将两者都转换为 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 获取最近的邻居。
我试过两种方法。两者都运行得非常快,结果可行。但是,它们并不准确。
方法一:
在这个页面中,我使用了最近邻方法sklearn.neighbors
。这将返回以米为单位的结果。但是,当我从两个数据帧手动检查 lat long 之间的距离时,我总是发现最近的邻居返回 1/4 的距离。
比如上面方法返回的距离是125米,google map和https://www.geodatasource.com/distance-calculator都返回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)
以上运行速度非常快并返回结果。然而,返回的距离值至少比我得到的低 100 倍。
乘数:107.655914
在上面的excel图片中,第一列表示python返回的结果,而第二列表示上面给出的同一网站返回的结果。虽然结果中的这些近似值让我开始,但我想要准确的结果。如何比较上面给出的两个数据框,并为 DF1 中的每一行获取最准确的最近距离。
解决方案
处理空间数据时,您应该注意点坐标是从球体投影到平面中的。在墨卡托投影中,纬度点之间的距离以度为单位,而不是米。并且转换取决于点的纬度,因为赤道的 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
推荐阅读
- xml - XSLT,将一些带有前缀的标签从 lt/gt 更改为 <>,而有些标签仍然带有 lt/gt
- java - 带有图像数据的 Java Volley 请求多部分表单
- java - 使用powershell修改部分JAVA_HOME文件夹路径
- html - 该页面在移动模式下显示屏幕外的隐藏项目
- java - 如何搜索 FTS3 表,以便搜索到的单词显示在列表顶部
- php - 找不到 php 命令,但 php70 正在工作
- observable - takeFirst 方法的替代方法是什么
- javascript - 在javascript中自动插入两个值的总和
- youtube-api - 偶尔会收到 YouTube iframe 播放器错误消息“发生错误。请稍后再试”
- node.js - 如何在 mongodb 中合并两个查询?