GeoPandy: znajdź najbliższy punkt w innej ramce danych

20

Mam 2 ramki geodata:

import geopandas as gpd
from shapely.geometry import Point
gpd1 = gpd.GeoDataFrame([['John',1,Point(1,1)],['Smith',1,Point(2,2)],['Soap',1,Point(0,2)]],columns=['Name','ID','geometry'])
gpd2 = gpd.GeoDataFrame([['Work',Point(0,1.1)],['Shops',Point(2.5,2)],['Home',Point(1,1.1)]],columns=['Place','geometry'])

i chcę znaleźć nazwę najbliższego punktu w gpd2 dla każdego wiersza w gpd1:

desired_output = 

    Name  ID     geometry  Nearest
0   John   1  POINT (1 1)     Home
1  Smith   1  POINT (2 2)    Shops
2   Soap   1  POINT (0 2)     Work

Próbowałem uruchomić to za pomocą funkcji lambda:

gpd1['Nearest'] = gpd1.apply(lambda row: min_dist(row.geometry,gpd2)['Place'] , axis=1)

z

def min_dist(point, gpd2):

    geoseries = some_function()
    return geoseries
RedM
źródło
Ta metoda działała dla mnie: stackoverflow.com/questions/37402046/... spójrz na link
Johnny Cheesecutter

Odpowiedzi:

16

Możesz bezpośrednio użyć funkcji Shapely Najbliższe punkty (geometrie GeoSeries są geometriami Shapely):

from shapely.ops import nearest_points
# unary union of the gpd2 geomtries 
pts3 = gpd2.geometry.unary_union
def near(point, pts=pts3):
     # find the nearest point and return the corresponding Place value
     nearest = gpd2.geometry == nearest_points(point, pts)[1]
     return gpd2[nearest].Place.get_values()[0]
gpd1['Nearest'] = gpd1.apply(lambda row: near(row.geometry), axis=1)
gpd1
    Name  ID     geometry  Nearest
0   John   1  POINT (1 1)     Home
1  Smith   1  POINT (2 2)    Shops
2   Soap   1  POINT (0 2)     Work

Rozwinięcie

for i, row in gpd1.iterrows():
    print nearest_points(row.geometry, pts3)[0], nearest_points(row.geometry, pts3)[1]
 POINT (1 1) POINT (1 1.1)
 POINT (2 2) POINT (2.5 2)
 POINT (0 2) POINT (0 1.1)
gen
źródło
Coś mi nie działa i nie mogę tego rozgryźć. Funkcja zwraca pusty GeoSeries, mimo że geometria jest solidna. Na przykład: sample_point = gpd2.geometry.unary_union[400] / sample_point in gpd2.geometry Zwraca wartość True. gpd2.geometry == sample_point To wychodzi na wszystkie fałszywe.
robroc
Dodatek do powyższego: gpd2.geometry.geom_equals(sample_point)działa.
robroc
13

Jeśli masz duże ramki danych, zauważyłem, że metoda scipyindeksu przestrzennego cKDTree .queryzwraca bardzo szybkie wyniki wyszukiwania najbliższego sąsiada. Ponieważ wykorzystuje indeks przestrzenny, jego rzędy wielkości są szybsze niż zapętlanie przez ramkę danych, a następnie znajdowanie minimum wszystkich odległości. Jest także szybszy niż używanie foremnych nearest_pointsz RTree (metoda indeksu przestrzennego dostępna za pośrednictwem geopandas), ponieważ cKDTree pozwala na wektoryzację wyszukiwania, podczas gdy druga metoda tego nie robi.

Oto funkcja pomocnika, która zwróci odległość i „imię” najbliższego sąsiada gpd2z każdego punktu w gpd1. Zakłada się, że oba pliki gdf mają geometrykolumnę (punktów).

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

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

gpd1 = gpd.GeoDataFrame([['John', 1, Point(1, 1)], ['Smith', 1, Point(2, 2)],
                         ['Soap', 1, Point(0, 2)]],
                        columns=['Name', 'ID', 'geometry'])
gpd2 = gpd.GeoDataFrame([['Work', Point(0, 1.1)], ['Shops', Point(2.5, 2)],
                         ['Home', Point(1, 1.1)]],
                        columns=['Place', 'geometry'])

def ckdnearest(gdA, gdB):
    nA = np.array(list(zip(gdA.geometry.x, gdA.geometry.y)) )
    nB = np.array(list(zip(gdB.geometry.x, gdB.geometry.y)) )
    btree = cKDTree(nB)
    dist, idx = btree.query(nA, k=1)
    gdf = pd.concat(
        [gdA, gdB.loc[idx, gdB.columns != 'geometry'].reset_index(),
         pd.Series(dist, name='dist')], axis=1)
    return gdf

ckdnearest(gpd1, gpd2)

A jeśli chcesz znaleźć najbliższy punkt LineString, oto pełny działający przykład:

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

gpd1 = gpd.GeoDataFrame([['John', 1, Point(1, 1)],
                         ['Smith', 1, Point(2, 2)],
                         ['Soap', 1, Point(0, 2)]],
                        columns=['Name', 'ID', 'geometry'])
gpd2 = gpd.GeoDataFrame([['Work', LineString([Point(100, 0), Point(100, 1)])],
                         ['Shops', LineString([Point(101, 0), Point(101, 1), Point(102, 3)])],
                         ['Home',  LineString([Point(101, 0), Point(102, 1)])]],
                        columns=['Place', 'geometry'])


def ckdnearest(gdfA, gdfB, gdfB_cols=['Place']):
    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(gpd1, gpd2)
JHuw
źródło
Czy za pomocą tej metody można również podać najbliższy punkt na linii? Na przykład, aby przyciągnąć pozycję GPS do najbliższej ulicy.
hyperknot
Ta odpowiedź jest niesamowita! Jednak kod najbliższych punktów do linii powoduje błąd. Wygląda na to, że dla każdego punktu zwracana jest poprawna odległość od najbliższej linii, ale zwracany identyfikator linii jest nieprawidłowy. Myślę, że to obliczenia idx, ale jestem całkiem nowy w Pythonie, więc nie mogę sobie tego poradzić.
Shakedk
1

Domyśliłam się:

def min_dist(point, gpd2):
    gpd2['Dist'] = gpd2.apply(lambda row:  point.distance(row.geometry),axis=1)
    geoseries = gpd2.iloc[gpd2['Dist'].argmin()]
    return geoseries

Oczywiście krytyka jest mile widziana. Nie jestem fanem przeliczania gpd2 ['Dist'] dla każdego wiersza gpd1 ...

RedM
źródło
1

Odpowiedź Gene'a nie zadziałała dla mnie. W końcu odkryłem, że gpd2.geometry.unary_union zaowocowało geometrią, która zawierała tylko około 30 000 z moich ogółem około 150 000 punktów. Dla każdego, kto napotka ten sam problem, oto jak go rozwiązałem:

    from shapely.ops import nearest_points
    from shapely.geometry import MultiPoint

    gpd2_pts_list = gpd2.geometry.tolist()
    gpd2_pts = MultiPoint(gpd2_pts_list)
    def nearest(point, gpd2_pts, gpd2=gpd2, geom_col='geometry', src_col='Place'):
         # find the nearest point
         nearest_point = nearest_points(point, gpd2_pts)[1]
         # return the corresponding value of the src_col of the nearest point
         value = gpd2[gpd2[geom_col] == nearest_point][src_col].get_values()[0]
         return value

    gpd1['Nearest'] = gpd1.apply(lambda x: nearest(x.geometry, gpd2_pts), axis=1)
Inske
źródło
0

Dla każdego, kto ma błędy indeksowania z własnymi danymi podczas korzystania z doskonałej odpowiedzi z @ JHuw , moim problemem było to, że moje indeksy nie były wyrównane. Zresetowanie indeksu gdfA i gdfB rozwiązało moje problemy, być może może to również pomóc @ Shakedk .

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

gpd1 = gpd.GeoDataFrame([['John', 1, Point(1, 1)],
                         ['Smith', 1, Point(2, 2)],
                         ['Soap', 1, Point(0, 2)]],
                        columns=['Name', 'ID', 'geometry'])
gpd2 = gpd.GeoDataFrame([['Work', LineString([Point(100, 0), Point(100, 1)])],
                         ['Shops', LineString([Point(101, 0), Point(101, 1), Point(102, 3)])],
                         ['Home',  LineString([Point(101, 0), Point(102, 1)])]],
                        columns=['Place', 'geometry'])


def ckdnearest(gdfA, gdfB, gdfB_cols=['Place']):
    # 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(gpd1, gpd2)
Markus Rosenfelder
źródło