在pandas中使用两个地理数据框获取最近距离

14

这是我的第一个地理数据框:

!pip install geopandas
import pandas as pd
import geopandas

city1 = [{'City':"Buenos Aires","Country":"Argentina","Latitude":-34.58,"Longitude":-58.66},
           {'City':"Brasilia","Country":"Brazil","Latitude":-15.78 ,"Longitude":-70.66},
         {'City':"Santiago","Country":"Chile ","Latitude":-33.45 ,"Longitude":-70.66 }]
city2 =  [{'City':"Bogota","Country":"Colombia ","Latitude":4.60 ,"Longitude":-74.08},
           {'City':"Caracas","Country":"Venezuela","Latitude":10.48  ,"Longitude":-66.86}]
city1df = pd.DataFrame(city1)
city2df = pd.DataFrame(city2)
gcity1df = geopandas.GeoDataFrame(
    city1df, geometry=geopandas.points_from_xy(city1df.Longitude, city1df.Latitude))
gcity2df = geopandas.GeoDataFrame(
    city2df, geometry=geopandas.points_from_xy(city2df.Longitude, city2df.Latitude))

城市1

           City    Country  Latitude  Longitude                     geometry
0  Buenos Aires  Argentina    -34.58     -58.66  POINT (-58.66000 -34.58000)
1      Brasilia     Brazil    -15.78     -47.91  POINT (-47.91000 -15.78000)
2      Santiago      Chile    -33.45     -70.66  POINT (-70.66000 -33.45000)

我的第二个地理数据框:

城市2 :

         City    Country  Latitude  Longitude                     geometry
1        Bogota   Colombia      4.60     -74.08    POINT (-74.08000 4.60000)
2       Caracas  Venezuela     10.48     -66.86   POINT (-66.86000 10.48000)

我希望获得第三个数据框,其中包含从城市1到城市2最近的城市和距离:

           City    Country  Latitude  Longitude                     geometry    Nearest    Distance
0  Buenos Aires  Argentina    -34.58     -58.66  POINT (-58.66000 -34.58000)    Bogota    111 Km

这是我使用geodjango和dict的实际解决方案(但太长了):

from django.contrib.gis.geos import GEOSGeometry
result = []
dict_result = {}
for city01 in city1 :
  dist = 99999999
  pnt = GEOSGeometry('SRID=4326;POINT( '+str(city01["Latitude"])+' '+str(city01['Longitude'])+')')
  for city02 in city2:
    pnt2 = GEOSGeometry('SRID=4326;POINT('+str(city02['Latitude'])+' '+str(city02['Longitude'])+')')
    distance_test = pnt.distance(pnt2) * 100
    if distance_test < dist :
      dist = distance_test
  result.append(dist)
  dict_result[city01['City']] = city02['City']

这是我的尝试:

from shapely.ops import nearest_points
# unary union of the gpd2 geomtries 
pts3 = gcity2df.geometry.unary_union
def Euclidean_Dist(df1, df2, cols=['x_coord','y_coord']):
    return np.linalg.norm(df1[cols].values - df2[cols].values,
                   axis=1)
def near(point, pts=pts3):
     # find the nearest point and return the corresponding Place value
     nearest = gcity2df.geometry == nearest_points(point, pts)[1]

     return gcity2df[nearest].City
gcity1df['Nearest'] = gcity1df.apply(lambda row: near(row.geometry), axis=1)
gcity1df
    City    Country     Latitude    Longitude   geometry    Nearest
0   Buenos Aires    Argentina   -34.58  -58.66  POINT (-58.66000 -34.58000)     Bogota
1   Brasilia    Brazil  -15.78  -70.66  POINT (-70.66000 -15.78000)     Bogota
2   Santiago    Chile   -33.45  -70.66  POINT (-70.66000 -33.45000)     Bogota

敬意


你好,欢迎来到StackOverflow!你似乎认为StackOverflow是一个你发布问题然后得到一些代码的网站。实际上并非如此。你的问题很可能会很快被关闭甚至删除。为了避免这种情况发生,请参观一下我们的导览查看一下帮助中心。特别是,让自己熟悉一下这里所认可的话题 - azro
此外,当您发布有关DF的内容时,请附上包含DF内容的Python代码,以便所有愿意帮助您的人不必自己编写。 - azro
@azro 我已经编辑并添加了我的解决方案以及初始数据。 - Bussiere
你们的城市只在南美洲吗?如果不是,它们之间可以有多远的距离?city1中可以有多少个城市,city2中可以有多少个城市?找到最快的解决方案很重要吗,还是一个在合理时间内运行的简单解决方案就可以了?如果是后者,那么合理的时间是多久? - Walter Tross
@WalterTross 我的城市遍布全球,我正在寻找最快的解决方案。谢谢。 - Bussiere
显示剩余2条评论
3个回答

13

首先,我通过交叉连接合并了两个数据框。然后,在Python中使用 map 来查找两点之间的距离。我使用 map,因为大多数情况下它比 applyitertuplesiterrows 等方法更快。(参考:https://dev59.com/aFQJ5IYBdhLWcg3w_7EN#52674448

最后,我按数据框进行分组,并获取距离的最小值。

这是相关的库:

import pandas as pd
import geopandas
import geopy.distance
from math import radians, cos, sin, asin, sqrt

这里展示了使用的函数:

def dist1(p1, p2):
    lon1, lat1, lon2, lat2 = map(radians, [p1.x, p1.y, p2.x, p2.y])

    dlon = lon2 - lon1 
    dlat = lat2 - lat1 
    a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
    c = 2 * asin(sqrt(a)) 

    return c * 6373

def dist2(p1, p2):
    lon1, lat1, lon2, lat2 = map(radians, [p1[0], p1[1], p2[0], p2[1]])

    dlon = lon2 - lon1 
    dlat = lat2 - lat1 
    a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
    c = 2 * asin(sqrt(a)) 

    return c * 6373

def dist3(p1, p2):
    x = p1.y, p1.x
    y = p2.y, p2.x

    return geopy.distance.geodesic(x, y).km

def dist4(p1, p2):
    x = p1[1], p1[0]
    y = p2[1], p2[0]

    return geopy.distance.geodesic(x, y).km

还有数据,

city1 = [
  {
    'City': 'Buenos Aires',
    'Country': 'Argentina',
    'Latitude': -34.58,
    'Longitude': -58.66
  },
  {
    'City': 'Brasilia',
    'Country': 'Brazil',
    'Latitude': -15.78,
    'Longitude': -70.66
  },
  {
    'City': 'Santiago',
    'Country': 'Chile ',
    'Latitude': -33.45,
    'Longitude': -70.66
  }
]

city2 = [
  {
    'City': 'Bogota',
    'Country': 'Colombia ',
    'Latitude': 4.6,
    'Longitude': -74.08
  },
  {
    'City': 'Caracas',
    'Country': 'Venezuela',
    'Latitude': 10.48,
    'Longitude': -66.86
  }
]


city1df = pd.DataFrame(city1)
city2df = pd.DataFrame(city2)

使用 geopandas 数据框进行交叉连接(Cross join)。

gcity1df = geopandas.GeoDataFrame(
    city1df, 
    geometry=geopandas.points_from_xy(city1df.Longitude, city1df.Latitude)
)
gcity2df = geopandas.GeoDataFrame(
    city2df, 
    geometry=geopandas.points_from_xy(city2df.Longitude, city2df.Latitude)
)

# cross join geopandas
gcity1df['key'] = 1
gcity2df['key'] = 1
merged = gcity1df.merge(gcity2df, on='key')

math函数和geopandas

# 6.64 ms ± 588 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%%timeit

# find distance
merged['dist'] = list(map(dist1, merged['geometry_x'], merged['geometry_y']))

mapping = {
    'City_x': 'City',
    'Country_x': 'Country',
    'Latitude_x': 'Latitude',
    'Longitude_x': 'Longitude',
    'geometry_x': 'geometry',
    'City_y': 'Nearest',
    'dist': 'Distance'
}

nearest = merged.loc[merged.groupby(['City_x', 'Country_x'])['dist'].idxmin()]
nearest.rename(columns=mapping)[list(mapping.values())]

           City    Country  Latitude  Longitude                     geometry  \
2      Brasilia     Brazil    -15.78     -70.66  POINT (-70.66000 -15.78000)   
0  Buenos Aires  Argentina    -34.58     -58.66  POINT (-58.66000 -34.58000)   
4      Santiago     Chile     -33.45     -70.66  POINT (-70.66000 -33.45000)   

  Nearest     Distance  
2  Bogota  2297.922808  
0  Bogota  4648.004515  
4  Bogota  4247.586882 
< p > geopygeopandas

# 9.99 ms ± 764 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%%timeit

# find distance
merged['dist'] = list(map(dist3, merged['geometry_x'], merged['geometry_y']))

mapping = {
    'City_x': 'City',
    'Country_x': 'Country',
    'Latitude_x': 'Latitude',
    'Longitude_x': 'Longitude',
    'geometry_x': 'geometry',
    'City_y': 'Nearest',
    'dist': 'Distance'
}

nearest = merged.loc[merged.groupby(['City_x', 'Country_x'])['dist'].idxmin()]
nearest.rename(columns=mapping)[list(mapping.values())]

           City    Country  Latitude  Longitude                     geometry  \
2      Brasilia     Brazil    -15.78     -70.66  POINT (-70.66000 -15.78000)   
0  Buenos Aires  Argentina    -34.58     -58.66  POINT (-58.66000 -34.58000)   
4      Santiago     Chile     -33.45     -70.66  POINT (-70.66000 -33.45000)   

  Nearest     Distance  
2  Bogota  2285.239605  
0  Bogota  4628.641817  
4  Bogota  4226.710978 

如果你想使用 pandas 而不是 geopandas

# cross join pandas
city1df['key'] = 1
city2df['key'] = 1
merged = city1df.merge(city2df, on='key')

使用math函数,

# 8.65 ms ± 2.21 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)
%%timeit

# find distance
merged['dist'] = list(
    map(
        dist2, 
        merged[['Longitude_x', 'Latitude_x']].values, 
        merged[['Longitude_y', 'Latitude_y']].values
    )
)

mapping = {
    'City_x': 'City',
    'Country_x': 'Country',
    'Latitude_x': 'Latitude',
    'Longitude_x': 'Longitude',
    'City_y': 'Nearest',
    'dist': 'Distance'
}

nearest = merged.loc[merged.groupby(['City_x', 'Country_x'])['dist'].idxmin()]
nearest.rename(columns=mapping)[list(mapping.values())]

           City    Country  Latitude  Longitude Nearest     Distance
2      Brasilia     Brazil    -15.78     -70.66  Bogota  2297.922808
0  Buenos Aires  Argentina    -34.58     -58.66  Bogota  4648.004515
4      Santiago     Chile     -33.45     -70.66  Bogota  4247.586882

使用 geopy

# 9.8 ms ± 807 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%%timeit

# find distance
merged['dist'] = list(
    map(
        dist4, 
        merged[['Longitude_x', 'Latitude_x']].values, 
        merged[['Longitude_y', 'Latitude_y']].values
    )
)

mapping = {
    'City_x': 'City',
    'Country_x': 'Country',
    'Latitude_x': 'Latitude',
    'Longitude_x': 'Longitude',
    'City_y': 'Nearest',
    'dist': 'Distance'
}

nearest = merged.loc[merged.groupby(['City_x', 'Country_x'])['dist'].idxmin()]
nearest.rename(columns=mapping)[list(mapping.values())]

           City    Country  Latitude  Longitude Nearest     Distance
2      Brasilia     Brazil    -15.78     -70.66  Bogota  2285.239605
0  Buenos Aires  Argentina    -34.58     -58.66  Bogota  4628.641817
4      Santiago     Chile     -33.45     -70.66  Bogota  4226.710978

1
这些距离是使用近似公式计算的,该公式没有考虑地球的扁平化。使用 geopy.distance.distance() 计算出的相同三个距离为(四舍五入)228546294227 公里。 - Walter Tross
我用链接 https://www.distance.to/-33.45,-70.66/4.6,-74.08 检查这些值,有什么问题吗? - E. Zeytinci
1
除了我更信任geopy之外,作为一个网站,我更信任http://edwilliams.org/gccalc.htm,它与`geopy`一致。NOAA的网站https://www.nhc.noaa.gov/gccalc.shtml说它是基于前者的,但是却得出不同的结果。它可能是基于前者的旧版本。 - Walter Tross

5
我认为很难找到一个时间复杂度O(m·n)更好的解决方案,其中m和n是city1city2的大小。保持距离比较(唯一的O(m·n)操作)简单,并利用numpy和pandas提供的向量化操作,速度对于任何合理的输入大小都不应该是问题。
这个想法是,在球面上比较距离,可以比较3D中点之间的距离。最近的城市也是通过球体的最近城市。此外,通常需要取平方根来计算距离,但如果只需要比较它们,可以避免平方根。
from geopy.distance import distance as dist
import numpy as np
import pandas as pd

def find_closest(lat1, lng1, lat2, lng2):
    def x_y_z_of_lat_lng_on_unit_sphere(lat, lng):
        rad_lat, rad_lng = np.radians(lat), np.radians(lng)
        sin_lat, sin_lng = np.sin(rad_lat), np.sin(rad_lng)
        cos_lat, cos_lng = np.cos(rad_lat), np.cos(rad_lng)
        return cos_lat * cos_lng, cos_lat * sin_lng, sin_lat
    x1, y1, z1 = x_y_z_of_lat_lng_on_unit_sphere(lat1, lng1)
    x2, y2, z2 = x_y_z_of_lat_lng_on_unit_sphere(lat2, lng2)
    return pd.Series(map(lambda x, y, z:
                         ((x2-x)**2 + (y2-y)**2 + (z2-z)**2).idxmin(),
                         x1, y1, z1))

city1 = [{"City":"Tokyo",    "Ctry":"JP", "Latitude": 35.68972, "Longitude": 139.69222},
         {"City":"Pretoria", "Ctry":"ZA", "Latitude":-25.71667, "Longitude": 28.28333},
         {"City":"London",   "Ctry":"GB", "Latitude": 51.50722, "Longitude": -0.12574}]
city2 = [{"City":"Seattle",  "Ctry":"US", "Latitude": 47.60972, "Longitude":-122.33306},
         {"City":"Auckland", "Ctry":"NZ", "Latitude":-36.84446, "Longitude": 174.76364}]
city1df = pd.DataFrame(city1)
city2df = pd.DataFrame(city2)

closest = find_closest(city1df.Latitude, city1df.Longitude, city2df.Latitude, city2df.Longitude)

resultdf = city1df.join(city2df, on=closest, rsuffix='2')
km = pd.Series(map(lambda latlng1, latlng2: round(dist(latlng1, latlng2).km),
                   resultdf[['Latitude',  'Longitude' ]].to_numpy(),
                   resultdf[['Latitude2', 'Longitude2']].to_numpy()))
resultdf['Distance'] = km
print(resultdf.to_string())
#        City Ctry  Latitude  Longitude     City2 Ctry2  Latitude2  Longitude2  Distance
# 0     Tokyo   JP  35.68972  139.69222   Seattle    US   47.60972  -122.33306      7715
# 1  Pretoria   ZA -25.71667   28.28333  Auckland    NZ  -36.84446   174.76364     12245
# 2    London   GB  51.50722   -0.12574   Seattle    US   47.60972  -122.33306      7723

请注意,任何将纬度和经度视为笛卡尔坐标的解决方案都是错误的,因为向极地移动时,经线(等经度线)会相互靠近。

3

这个解决方案可能不是解决问题的最快方法,但我相信它可以解决问题。

#New dataframe is basicly a copy of first but with more columns
gcity3df = gcity1df.copy()
gcity3df["Nearest"] = None
gcity3df["Distance"] = None

#For each city (row in gcity3df) we will calculate the nearest city from gcity2df and 
fill the Nones with results

for index, row in gcity3df.iterrows():
    #Setting neareast and distance to None, 
    #we will be filling those variables with results

    nearest = None
    distance = None
    for df2index, df2row in gcity2df.iterrows():
        d = row.geometry.distance(df2row.geometry)
        #If df2index city is closer than previous ones, replace nearest with it
        if distance is None or d < distance:
            distance = d
            nearest = df2row.City 
    #In the end we appends the closest city to gdf
    gcity3df.at[index, "Nearest"] = nearest
    gcity3df.at[index, "Distance"] = distance

如果您需要使用米而不是度来工作,您可以重新投影您的图层(这也会消除Walter所提到的错误)。您可以通过 gcity3df = gcity3df.to_crs({'init': 'epsg:XXXX'}) 来实现,其中XXXX是您所在世界地区使用的crs的EPSG代码。


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接