GeoPandas: Encuentra el punto más cercano en otro marco de datos

20

Tengo 2 geodataframes:

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'])

y quiero encontrar el nombre del punto más cercano en gpd2 para cada fila en 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

He estado tratando de hacer que esto funcione usando una función lambda:

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

con

def min_dist(point, gpd2):

    geoseries = some_function()
    return geoseries
RedM
fuente
Este método funcionó para mí: stackoverflow.com/questions/37402046/… mira el enlace
Johnny Cheesecutter

Respuestas:

16

Puede usar directamente la función Shapely Puntos más cercanos (las geometrías de GeoSeries son geometrías 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

Explicación

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)
gene
fuente
Algo no funciona para mí y no puedo entenderlo. La función devuelve una GeoSeries vacía aunque la geometría sea sólida. Por ejemplo: sample_point = gpd2.geometry.unary_union[400] / sample_point in gpd2.geometry Esto devuelve True. gpd2.geometry == sample_point Esto sale todo falso.
robroc
Además de lo anterior: gpd2.geometry.geom_equals(sample_point)obras.
robroc
13

Si tiene grandes marcos de datos, he descubierto que scipyel .querymétodo de índice espacial cKDTree devuelve resultados muy rápidos para las búsquedas de vecinos más cercanos. Como utiliza un índice espacial, es de órdenes de magnitud más rápido que recorrer el marco de datos y luego encontrar el mínimo de todas las distancias. También es más rápido que usar nearest_pointsShapely 's con RTree (el método de índice espacial disponible a través de geopandas) porque cKDTree le permite vectorizar su búsqueda, mientras que el otro método no.

Aquí hay una función auxiliar que devolverá la distancia y el 'Nombre' del vecino más cercano gpd2desde cada punto gpd1. Se supone que ambos gdfs tienen una geometrycolumna (de puntos).

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)

Y si desea encontrar el punto más cercano a LineString, aquí hay un ejemplo de trabajo completo:

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
fuente
¿Es posible dar también el punto más cercano en la línea, usando este método? Por ejemplo, para ajustar una ubicación GPS a la calle más cercana.
hyperknot
¡Esta respuesta es asombrosa! Sin embargo, el código para los puntos más cercanos a la línea produce un error para mí. Parece que se devuelve la distancia correcta desde la línea más cercana para cada punto, pero la identificación de la línea que se devuelve es incorrecta. Creo que es el cálculo de idx, pero soy bastante nuevo en Python, por lo que no puedo entenderlo.
Shakedk
1

Lo averigué:

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

Por supuesto, algunas críticas son bienvenidas. No soy fanático de volver a calcular gpd2 ['Dist'] para cada fila de gpd1 ...

RedM
fuente
1

La respuesta de Gene no funcionó para mí. Finalmente descubrí que gpd2.geometry.unary_union resultó en una geometría que solo contenía aproximadamente 30,000 de mi total de aproximadamente 150,000 puntos. Para cualquier otra persona que tenga el mismo problema, así es como lo resolví:

    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
fuente
0

Para cualquiera que tenga errores de indexación con sus propios datos mientras usa la excelente respuesta de @ JHuw , mi problema fue que mis índices no se alinearon. Restablecer el índice de gdfA y gdfB resolvió mis problemas, tal vez esto también pueda ayudarlo @ 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
fuente