Поиск ближайшей станции к каждому магазину с помощью BallTree - PullRequest
0 голосов
/ 12 июля 2020

У меня есть 2 набора данных, список магазинов с координатами в Великобритании и железнодорожный вокзал с координатами.

Я использую BallTree, чтобы найти ближайшую станцию ​​к каждому магазину с указанием расстояния, используя aa код с этого веб-сайта, и я соответствующим образом поменял свои фреймы данных.

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

Код:

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

from sklearn.neighbors import BallTree

df_pocs = pd.read_csv(r'C:\Users\FLETCHWI\Desktop\XX\shops.csv', encoding = "ISO-8859-1", engine='python')

df_stations = pd.read_csv(r'C:\Users\FLETCHWI\Desktop\xx\uk_stations.csv', encoding = "ISO-8859-1", engine='python')


gdf_pocs = gpd.GeoDataFrame(
    df_pocs, geometry=gpd.points_from_xy(df_pocs.longitude, df_pocs.latitude))

gdf_stations = gpd.GeoDataFrame(
    df_stations, geometry=gpd.points_from_xy(df_stations.longitude, df_stations.latitude))



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

# Find closest public transport stop for each building and get also the distance based on haversine distance
# Note: haversine distance which is implemented here is a bit slower than using e.g. 'euclidean' metric
# but useful as we get the distance between points in meters
closest_stations = nearest_neighbor(gdf_pocs, gdf_stations, return_dist=True)

После запуска кода , он возвращает одну и ту же станцию ​​для каждого магазина, который у меня есть. Однако я бы хотел, чтобы он находил ближайшую станцию ​​для каждого магазина и расстояние до нее.

Любая помощь приветствуется, спасибо!

введите описание изображения здесь

1 Ответ

0 голосов
/ 13 июля 2020

Я провел некоторое тестирование функций, и действительно, чтобы оно работало, необходимо поменять местами lat / long.

Обратите внимание на предупреждение:

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

Следовательно, при определении точки simple изменить

gdf_pocs = gpd.GeoDataFrame(
    df_pocs, geometry=gpd.points_from_xy(df_pocs.longitude, df_pocs.latitude))

gdf_stations = gpd.GeoDataFrame(
    df_stations, geometry=gpd.points_from_xy(df_stations.longitude, df_stations.latitude))

на

gdf_pocs = gpd.GeoDataFrame(
    df_pocs, geometry=gpd.points_from_xy(df_pocs.latitude, df_pocs.longitude))

gdf_stations = gpd.GeoDataFrame(
    df_stations, geometry=gpd.points_from_xy(df_stations.latitude, df_stations.longitude))

...