Finding nearest point in other GeoDataFrame using GeoPandas

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 g`pd1`:

``````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
``````

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. 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. 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. Commented Apr 26, 2021 at 17:58

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)
``````
• 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. Commented Aug 21, 2018 at 15:47
• Addition to above: `gpd2.geometry.geom_equals(sample_point)` works. Commented Aug 21, 2018 at 16:26

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.

• 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 Commented Oct 12, 2023 at 11:27

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...

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? 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/… Commented Mar 8, 2021 at 13:50

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)
``````

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

``````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
``````

``````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)

# 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)

if return_dist:
# Convert to meters from radians

return closest_points
``````

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"
)
``````

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
``````