Finding nearest point in other GeoDataFrame using GeoPandas

Question

I've got two 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'])

and I want to find the name of the nearest point in gpd2 for each row in 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

I've been trying to get this working using a lambda function:

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

with

def min_dist(point, gpd2):
    geoseries = some_function(...)
    return geoseries

This method worked for me: stackoverflow.com/questions/37402046/… look at the link — Johnny Cheesecutter, Commented Mar 25, 2019 at 20:07

Ric S · Accepted Answer · 2021-02-25 08:12:30Z

If you have large dataframes, I've found that scipy's cKDTree spatial index .query method returns very fast results for nearest neighbor searches. As it uses a spatial index it's orders of magnitude faster than looping though the dataframe and then finding the minimum of all distances. It is also faster than using shapely's nearest_points with RTree (the spatial index method available via geopandas) because cKDTree allows you to vectorize your search whereas the other method does not.

Here is a helper function that will return the distance and 'Name' of the nearest neighbor in gpd2 from each point in gpd1. It assumes both gdfs have a geometry column (of points).


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(gdA.geometry.apply(lambda x: (x.x, x.y))))
    nB = np.array(list(gdB.geometry.apply(lambda x: (x.x, x.y))))
    btree = cKDTree(nB)
    dist, idx = btree.query(nA, k=1)
    gdB_nearest = gdB.iloc[idx].drop(columns="geometry").reset_index(drop=True)
    gdf = pd.concat(
        [
            gdA.reset_index(drop=True),
            gdB_nearest,
            pd.Series(dist, name='dist')
        ], 
        axis=1)

    return gdf

ckdnearest(gpd1, gpd2)

And if you want to find the closest point to a LineString, here is a full working example:

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)

Is it possible to give the nearest point on the line as well, using this method? For example to snap a GPS location to the nearest street. — hyperknot, Commented Dec 5, 2018 at 16:46
This answer is amazing! However, the code for nearest points to line produces a bug for me. It seems the correct distance from the closest line is returned for each point, but the line id that that is returned is wrong. I think its the idx calculation, but I'm pretty new to Python, so I can't manage to wrap my head around it. — Shakedk, Commented Nov 19, 2019 at 19:45
What is the gdfB_cols parameter for? It will select the columns of LineString to be concatenated with the Point? RecursionError: maximum recursion depth exceeded while calling a Python object and import sys; sys.setrecursionlimit(10000), crashing. Are there something to optimise, improve in the nearest point from Point to LineString? — hhh, Commented Nov 19, 2020 at 1:34
This extracts the points from the paths right, and uses them as "proxies" for the line? But seems like that wont work if the points on the line is far away from each other, but you are close to the line. — epa095, Commented Apr 26, 2021 at 17:58

gene · Accepted Answer · 2016-12-22 20:23:43Z

32

You can directly use the Shapely function Nearest points (the geometries of the GeoSeries are Shapely geometries):

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

Explication

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)

edited Dec 22, 2016 at 20:23

answered Dec 22, 2016 at 20:16

gene

55.4k3 gold badges113 silver badges191 bronze badges

Something isn't working for me and I can't figure it out. The function returns an empty GeoSeries even though the geometry is solid. For example: sample_point = gpd2.geometry.unary_union[400] / sample_point in gpd2.geometry This returns True. gpd2.geometry == sample_point This comes out all False.
– robroc
Commented Aug 21, 2018 at 15:47
1

Addition to above: gpd2.geometry.geom_equals(sample_point) works.
– robroc
Commented Aug 21, 2018 at 16:26

Add a comment |

do-me · Accepted Answer · 2023-10-13 13:59:30Z

21

As of v0.10.0 geopandas supports sjoin_nearest natively - see here.

Example: Get the nearest distances between points of two different GeoDataFrames or within one GeoDataFrame. Note that you should use the right EPSG for your region.

import geopandas as gpd 
gdf1, gdf2 = # load here and convert to the right EPSG
gdf1 = gdf1.to_crs("32633") # transverse mercator
gdf2 = gdf2.to_crs("32633")

gpd.sjoin_nearest(gdf1, gdf2, distance_col="distances",
    lsuffix="left", rsuffix="right", exclusive=True)

Use exclusive=True to exclude 0 distances between same points.

edited Oct 13, 2023 at 13:59

do-me

1033 bronze badges

answered Oct 8, 2021 at 12:54

rdmolony

3332 silver badges5 bronze badges

1

That looks pretty good. I can't believe this question still gets hits though.
– RedM
Commented Oct 8, 2021 at 13:17
absolutely the best answer in 2023
– wfgeo
Commented Oct 12, 2023 at 11:27

Add a comment |

RedM · Accepted Answer · 2016-12-22 09:15:16Z

9

Figured it out:

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

Of course some criticism is welcome. I'm not a fan of recalculating gpd2['Dist'] for every row of gpd1...

answered Dec 22, 2016 at 9:15

RedM

1,6531 gold badge11 silver badges17 bronze badges

Add a comment |

Markus Rosenfelder · Accepted Answer · 2019-11-30 20:15:02Z

For anyone having indexing errors with their own data while using the excellent answer from @JHuw, my problem was that my indexes did not align. Resetting the index of gdfA and gdfB solved my issues, maybe this can help you as well @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)

in case of comparing the distance between two pairs of lat and long how to read the distance? — Apricot, Commented Mar 4, 2021 at 15:36
You can use GeoSeries.distance() to get the distance between two points. See geopandas.org/docs/reference/api/… — Markus Rosenfelder, Commented Mar 8, 2021 at 13:50

Kadir Şahbaz · Accepted Answer · 2020-06-13 13:54:33Z

The answer by Gene didn't work for me. Finally I discovered that gpd2.geometry.unary_union resulted in a geometry that only contained about 30.000 of my total of roughly 150.000 points. For anyone else running into the same problem, here's how I solved it:

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)

The line has to be value = gpd2[gpd2[geom_col] == nearest_point][src_col].values[0], as the get_values() is deprecated pandas.pydata.org/pandas-docs/version/0.25.3/reference/api/… — nish, Commented Mar 18, 2022 at 13:04

rdmolony · Accepted Answer · 2021-03-08 13:53:31Z

If interested in an interactive tutorial on using either shapely nearest_points or scikit-learn BallTree larger data sets (instead of scipy KDTree) I found Automating GIS processes by the University of Helsinki, Finland) to be really useful

Their Github (the course notes update every year): https://github.com/Automating-GIS-processes/site

nearest_points (2020): https://autogis-site.readthedocs.io/en/latest/notebooks/L3/04_nearest-neighbour.html

from shapely.ops import nearest_points

def get_nearest_values(row, other_gdf, point_column='geometry', value_column="geometry"):
    """Find the nearest point and return the corresponding value from specified value column."""
    
    # Create an union of the other GeoDataFrame's geometries:
    other_points = other_gdf["geometry"].unary_union
    
    # Find the nearest points
    nearest_geoms = nearest_points(row[point_column], other_points)
    
    # Get corresponding values from the other df
    nearest_data = other_gdf.loc[other_gdf["geometry"] == nearest_geoms[1]]
    
    nearest_value = nearest_data[value_column].values[0]
    
    return nearest_value

BallTree (2020): https://autogis-site.readthedocs.io/en/latest/notebooks/L3/06_nearest-neighbor-faster.html

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

wfgeo · Accepted Answer · 2020-11-27 10:32:33Z

This solution is extremely inefficient but it should work for any and all geometry types (including mixed geometry type gdfs). I would only try this if your gdfs are small (my use case was a gdf with about 2000 rows for which I wanted to find the nearest feature from another gdf with about 15 rows, and it took a few seconds on a typical office laptop). Intersecting features may cause it to spaz out though, so be warned. It was originally based on @RedM's solution but will instead assign the index of the feature in gdf2 that is nearest to gdf1

gdf1["gdf2_idx"] = gdf1.apply(
    lambda row1: gdf2.apply(
        lambda row2: row1.geometry.distance(row2.geometry), axis="columns").idxmin(),
    axis="columns"
)

Vincent · Accepted Answer · 2021-04-19 22:51:36Z

Note, Gene's answer (https://gis.stackexchange.com/a/222388/157928) didn't work for me. I kept getting errors because of the following:

AttributeError: 'Series' object has no attribute 'get_values'

This worked for me:

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])
    nearest = nearest.index[nearest == True].tolist()
    return gpd2.loc[nearest].Place.iloc[0]
gpd1['Nearest'] = gpd1.apply(lambda row: near(row.geometry), axis=1)
gpd1

Stack Exchange Network

Finding nearest point in other GeoDataFrame using GeoPandas

9 Answers 9

Your Answer

Not the answer you're looking for? Browse other questions tagged
python
distance
geopandas
proximity
pandas
or ask your own question.

Linked

Hot Network Questions

Finding nearest point in other GeoDataFrame using GeoPandas

9 Answers 9

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged pythondistancegeopandasproximitypandas or ask your own question.

Linked

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
python
distance
geopandas
proximity
pandas
or ask your own question.