logo标签列表

最近邻连接并加入距离条件。

pythonpandasscikit-learngeopandas
5
在这个问题中，我指的是这个项目：

https://automating-gis-processes.github.io/site/master/notebooks/L3/nearest-neighbor-faster.html

我们有两个GeoDataFrame：
建筑物（Buildings）:
             name                   geometry
0            None  POINT (24.85584 60.20727)
1     Uimastadion  POINT (24.93045 60.18882)
2            None  POINT (24.95113 60.16994)
3  Hartwall Arena  POINT (24.92918 60.20570)

和公交车站：

     stop_name   stop_lat   stop_lon  stop_id                   geometry
0  Ritarihuone  60.169460  24.956670  1010102  POINT (24.95667 60.16946)
1   Kirkkokatu  60.171270  24.956570  1010103  POINT (24.95657 60.17127)
2   Kirkkokatu  60.170293  24.956721  1010104  POINT (24.95672 60.17029)
3    Vironkatu  60.172580  24.956554  1010105  POINT (24.95655 60.17258)

应用以下代码：

从sklearn.neighbors导入BallTree

from sklearn.neighbors import BallTree
import numpy as np

def get_nearest(src_points, candidates, k_neighbors=1):
    """Find nearest neighbors for all source points from a set of candidate points"""

    # Create tree from the candidate points
    tree = BallTree(candidates, leaf_size=15, metric='haversine')

    # Find closest points and distances
    distances, indices = tree.query(src_points, k=k_neighbors)

    # Transpose to get distances and indices into arrays
    distances = distances.transpose()
    indices = indices.transpose()

    # Get closest indices and distances (i.e. array at index 0)
    # note: for the second closest points, you would take index 1, etc.
    closest = indices[0]
    closest_dist = distances[0]

    # Return indices and distances
    return (closest, closest_dist)


def nearest_neighbor(left_gdf, right_gdf, return_dist=False):
    """
    For each point in left_gdf, find closest point in right GeoDataFrame and return them.

    NOTICE: Assumes that the input Points are in WGS84 projection (lat/lon).
    """

    left_geom_col = left_gdf.geometry.name
    right_geom_col = right_gdf.geometry.name

    # Ensure that index in right gdf is formed of sequential numbers
    right = right_gdf.copy().reset_index(drop=True)

    # Parse coordinates from points and insert them into a numpy array as RADIANS
    left_radians = np.array(left_gdf[left_geom_col].apply(lambda geom: (geom.x * np.pi / 180, geom.y * np.pi / 180)).to_list())
    right_radians = np.array(right[right_geom_col].apply(lambda geom: (geom.x * np.pi / 180, geom.y * np.pi / 180)).to_list())

    # Find the nearest points
    # -----------------------
    # closest ==> index in right_gdf that corresponds to the closest point
    # dist ==> distance between the nearest neighbors (in meters)

    closest, dist = get_nearest(src_points=left_radians, candidates=right_radians)

    # Return points from right GeoDataFrame that are closest to points in left GeoDataFrame
    closest_points = right.loc[closest]

    # Ensure that the index corresponds the one in left_gdf
    closest_points = closest_points.reset_index(drop=True)

    # Add distance if requested
    if return_dist:
        # Convert to meters from radians
        earth_radius = 6371000  # meters
        closest_points['distance'] = dist * earth_radius

            return closest_points


closest_stops = nearest_neighbor(buildings, stops, return_dist=True)

我们为每个建筑物索引获取到最近公交车站的距离：
    stop_name    stop_lat   stop_lon    stop_id                 geometry      distance
0   Muusantori   60.207490  24.857450   1304138 POINT (24.85745 60.20749)   180.521584
1   Eläintarha   60.192490  24.930840   1171120 POINT (24.93084 60.19249)   372.665221
2   Senaatintori 60.169010  24.950460   1020450 POINT (24.95046 60.16901)   119.425777
3   Veturitie    60.206610  24.929680   1174112 POINT (24.92968 60.20661)   106.762619

我正在寻找一种解决方案，以获取每个建筑物附近距离不超过250米的所有公交车站（可以多个）。谢谢您的帮助。
@s.k 是的和不是的，因为GIS工具无法处理大型数据集，这就是我正在寻找严格的Python解决方案的原因。 - cincin21
1
GIS工具是最早处理大数据的工具之一，所以我不确定你在评论中想表达什么。 - Paul H
我投票关闭此问题，因为它应该在gis.se上发布。 - Paul H
嗨@PaulH，我正在提到像scikit-learn、numpy这样的Python工具，它们可以帮助解决这个问题（以及其他类似的问题，不仅仅是空间问题）。如果没有人能够提供帮助，我可以将其转移到gis.se。关闭困难的问题并不是万能的答案... - cincin21
我不知道答案，但我认为你在GIS堆栈交换上会更有好运。这就是全部。这是免费的建议。按其价值接受它。如果你真的想让我考虑一下，我目前没有工作，收费$75/小时提供空间和统计分析咨询服务。 - Paul H
显示剩余3条评论
2个回答
5

有一种方法可以重复使用BallTree所做的事情，就像问题中那样，但是使用query_radius而不是函数格式，但你仍然可以轻松更改它。

from sklearn.neighbors import BallTree
import numpy as np
import pandas as pd
## here I start with buildings and stops as loaded in the link provided

# variable in meter you can change
radius_max = 250 # meters
# another parameter, in case you want to do with Mars radius ^^
earth_radius = 6371000  # meters

# similar to the method with apply in the tutorial 
# to create left_radians and right_radians, but faster
candidates = np.vstack([stops['geometry'].x.to_numpy(), 
                        stops['geometry'].y.to_numpy()]).T*np.pi/180
src_points = np.vstack([buildings['geometry'].x.to_numpy(), 
                        buildings['geometry'].y.to_numpy()]).T*np.pi/180

# Create tree from the candidate points
tree = BallTree(candidates, leaf_size=15, metric='haversine')
# use query_radius instead
ind_radius, dist_radius = tree.query_radius(src_points, 
                                            r=radius_max/earth_radius, 
                                            return_distance=True)

现在你可以操纵结果以获得你想要的内容。

# create a dataframe build with
# index based on row position of the building in buildings
# column row_stop is the row position of the stop
# dist is the distance
closest_dist = pd.concat([pd.Series(ind_radius).explode().rename('row_stop'), 
                          pd.Series(dist_radius).explode().rename('dist')*earth_radius], 
                         axis=1)
print (closest_dist.head())
#  row_stop     dist
#0     1131  180.522
#1      NaN      NaN
#2       64  174.744
#2       61  119.426
#3      532  106.763

# merge the dataframe created above with the original data stops
# to get names, id, ... note: the index must be reset as in closest_dist
# it is position based
closest_stop = closest_dist.merge(stops.reset_index(drop=True), 
                                  left_on='row_stop', right_index=True, how='left')
print (closest_stop.head())
#  row_stop     dist     stop_name  stop_lat  stop_lon    stop_id  \
#0     1131  180.522    Muusantori  60.20749  24.85745  1304138.0   
#1      NaN      NaN           NaN       NaN       NaN        NaN   
#2       64  174.744  Senaatintori  60.16896  24.94983  1020455.0   
#2       61  119.426  Senaatintori  60.16901  24.95046  1020450.0   
#3      532  106.763     Veturitie  60.20661  24.92968  1174112.0   
#
#                    geometry  
#0  POINT (24.85745 60.20749)  
#1                       None  
#2  POINT (24.94983 60.16896)  
#2  POINT (24.95046 60.16901)  
#3  POINT (24.92968 60.20661)

最后回到建筑物

# join buildings with reset_index with 
# closest_stop as index in closest_stop are position based
final_df = buildings.reset_index(drop=True).join(closest_stop, rsuffix='_stop')
print (final_df.head(10))
#              name                   geometry row_stop     dist     stop_name  \
# 0            None  POINT (24.85584 60.20727)     1131  180.522    Muusantori   
# 1     Uimastadion  POINT (24.93045 60.18882)      NaN      NaN           NaN   
# 2            None  POINT (24.95113 60.16994)       64  174.744  Senaatintori   
# 2            None  POINT (24.95113 60.16994)       61  119.426  Senaatintori   
# 3  Hartwall Arena  POINT (24.92918 60.20570)      532  106.763     Veturitie   

#    stop_lat  stop_lon    stop_id              geometry_stop  
# 0  60.20749  24.85745  1304138.0  POINT (24.85745 60.20749)  
# 1       NaN       NaN        NaN                       None  
# 2  60.16896  24.94983  1020455.0  POINT (24.94983 60.16896)  
# 2  60.16901  24.95046  1020450.0  POINT (24.95046 60.16901)  
# 3  60.20661  24.92968  1174112.0  POINT (24.92968 60.20661)
在最后一个DF中，为什么索引＃0和2x＃2处有“None”建筑物？而索引＃1的建筑物半径250m内没有公交车站，对吧？ - cincin21
1
@cincin21，关于None的问题，我在打开建筑物数据框时发现有一些没有名称的行（None），这就是为什么最后你也得到了None。对于第1行，所有的nan都表示在限制范围内（在这种情况下为250m）没有停靠站。 - Ben.T
1
@cincin21，如果你看一下你在问题中打印closest_stops的那一行，在末尾处，第一行距离为372，这证实了在此找到的结果。 - Ben.T
0

获取距离250米以内的最近公交车站：

filtered_by_distance = closest_stops[closest_stops.distance < 250]
result = buildings.join(filtered_by_distance)

对于半径内的所有站点，您需要使用BallTree.query_radius。但是您需要将米转换为弧度。

谢谢，我会阅读BallTree.query_radius的相关内容。第一个解决方案不适用于我的情况，因为我想要每个公交站（可能不止一个）。 - cincin21
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接