How to Find the Nearest Point in GeoPandas
On this page
Problem statement
"Which supermarket is nearest to each village?" is one call in GeoPandas, sjoin_nearest, and it always returns an answer. Whether the answer is right depends on two things the call does not check: whether distance is measured in metres or in degrees, and whether every feature in the target layer is what the layer's name says.
This page matches Rutland's 64 towns, villages and hamlets to their nearest OpenStreetMap supermarket, and both things turn out to matter: measured in degrees, one place in six gets a different supermarket, and the nearest "supermarket" for a third of the county is a military base.
Quick answer
The most practical way to find the nearest point in GeoPandas is usually GeoDataFrame.sjoin_nearest().
Basic pattern, for the nearest supermarket to each town, village and hamlet in Rutland (OpenStreetMap data, set up under Code examples):
places = places.to_crs("EPSG:27700") # British National Grid, in metres
supermarkets = supermarkets.to_crs(places.crs)
nearest = places.sjoin_nearest(
supermarkets,
how="left",
distance_col="distance_m"
)
print(f"{len(places)} places -> {len(nearest)} rows; median distance {nearest['distance_m'].median():,.0f} m")64 places -> 64 rows; median distance 4,824 m
Use a projected CRS before measuring distance. If your data is still in latitude/longitude such as EPSG:4326, the distance values are in degrees — Example 2 shows what that looks like.
Align the layers, sjoin_nearest, cap the distance
Load the point datasets
Start by reading the two layers. Here the places are points, and the supermarkets are a mix: OpenStreetMap maps most shops as building outlines, so the points layer has one supermarket and the polygons layer six.
import pandas as pd
print(places.geom_type.unique(), supermarkets_raw.geom_type.value_counts().to_dict())
# one point per supermarket, inside its outline
supermarkets = supermarkets_raw.to_crs(27700)
supermarkets["geometry"] = supermarkets.representative_point()
print(supermarkets[["name", "geometry"]].to_string())<ArrowStringArray>
['Point']
Length: 1, dtype: str {'Polygon': 6, 'Point': 1}
name geometry
0 Kendrew Barracks Headquarters POINT (491254.402 315129.476)
1 Tesco POINT (485892.848 308613.211)
2 The Co-operative Food POINT (486934.372 309847.769)
3 Co-op Food POINT (486760.029 299753.082)
4 The Co-operative Food POINT (485453.287 308617.001)
5 Aldi POINT (485425.322 310161.945)
6 Lidl POINT (485549.802 309351.697)
Check what the target layer contains before trusting its nearest matches: one of Rutland's seven "supermarkets" is named Kendrew Barracks Headquarters. sjoin_nearest finds the nearest feature; whether that feature is what the name of the layer says is a question about the data.
Check and align the coordinate reference system
Both layers must use the same CRS.
print(places.crs, supermarkets.crs)EPSG:27700 EPSG:27700
If the CRS differs, reproject one layer to match the other. For nearest-distance analysis, use a projected CRS with metre-based units: in Great Britain, EPSG:27700.
places = places.to_crs("EPSG:27700")
supermarkets = supermarkets.to_crs(places.crs)Find the nearest point with sjoin_nearest()
Use sjoin_nearest() to attach the nearest supermarket to each place.
nearest = places.sjoin_nearest(
supermarkets[["name", "geometry"]],
how="left",
distance_col="distance_m"
)
print(len(places), "places,", len(nearest), "rows")64 places, 64 rows
One row per place, because no place had two supermarkets at exactly the same distance; a tie would have given that place two rows. To cap the search, max_distance takes a distance in the CRS's units:
nearest_5km = places.sjoin_nearest(
supermarkets[["name", "geometry"]],
how="left",
distance_col="distance_m",
max_distance=5000
)
print(nearest_5km["index_right"].isna().sum(), "of", len(places), "places have no supermarket within 5 km")31 of 64 places have no supermarket within 5 km
Review the output columns
Inspect the result:
print(list(nearest.columns))
print(nearest[["name_left", "fclass", "name_right", "distance_m"]].round(0).head())['osm_id', 'code', 'fclass', 'population', 'name_left', 'geometry', 'index_right', 'name_right', 'distance_m']
name_left fclass name_right distance_m
0 Ryhall village Kendrew Barracks Headquarters 12992.0
1 Pickworth hamlet Kendrew Barracks Headquarters 8167.0
2 Cottesmore village Kendrew Barracks Headquarters 1642.0
3 Seaton village Co-op Food 3619.0
4 Thorpe by Water village Co-op Food 4132.0
Both layers have a name, so the result has name_left and name_right; index_right is the matched supermarket's row. And the first three places already show the problem Example 1 measures.
Filter by maximum search distance if needed
If you do not want unrealistic matches, you can filter after the join.
nearby_only = nearest[nearest["distance_m"] <= 5000].copy()If you used max_distance inside sjoin_nearest(), rows with no match inside that threshold have missing joined values:
unmatched = nearest_5km[nearest_5km["index_right"].isna()].copy()
print(len(nearby_only), "within 5 km after the join;", len(unmatched), "unmatched with max_distance")33 within 5 km after the join; 31 unmatched with max_distance
The two agree — max_distance just does the filtering inside the join, and keeps the unmatched places as rows with blanks rather than dropping them.
Rutland places to their nearest supermarket
The data used on this page
Geofabrik's OpenStreetMap extract for Rutland, from its dated 1 January 2026 snapshot so the counts below stay reproducible: the town, village and hamlet points from the places layer, and every feature tagged as a supermarket from the points and polygons layers.
import io
import zipfile
from pathlib import Path
import geopandas as gpd
import pandas as pd
import requests
if not Path("osm").exists():
reply = requests.get("https://download.geofabrik.de/europe/united-kingdom/england/"
"rutland-260101-free.shp.zip", timeout=300)
reply.raise_for_status()
zipfile.ZipFile(io.BytesIO(reply.content)).extractall("osm")
places = gpd.read_file("osm/gis_osm_places_free_1.shp")
places = places[places["fclass"].isin(["town", "village", "hamlet"])].reset_index(drop=True)
pois = [gpd.read_file(f"osm/gis_osm_{layer}_free_1.shp") for layer in ("pois", "pois_a")]
supermarkets_raw = pd.concat([p[p["fclass"] == "supermarket"] for p in pois], ignore_index=True)
supermarkets = supermarkets_raw.to_crs(27700)
supermarkets["geometry"] = supermarkets.representative_point()
print(len(places), "places:", places["fclass"].value_counts().to_dict(), "|", len(supermarkets), "supermarkets")64 places: {'village': 52, 'hamlet': 10, 'town': 2} | 7 supermarkets
Example 1: Find the nearest supermarket for each place
places_bng = places.to_crs("EPSG:27700")
nearest = places_bng.sjoin_nearest(
supermarkets[["name", "geometry"]],
how="left",
distance_col="distance_m"
)
print(nearest.sort_values("distance_m", ascending=False)
[["name_left", "fclass", "name_right", "distance_m"]].round(0).head(5).to_string(index=False))
print("\nplaces served by each supermarket:", nearest["name_right"].value_counts().to_dict())
shops = supermarkets[supermarkets["name"] != "Kendrew Barracks Headquarters"]
without = places_bng.sjoin_nearest(shops[["name", "geometry"]], how="left", distance_col="distance_m")
moved = nearest["name_right"].eq("Kendrew Barracks Headquarters")
print(f"without it: {moved.sum()} places re-matched; their median distance "
f"{nearest.loc[moved, 'distance_m'].median():,.0f} m -> {without.loc[moved, 'distance_m'].median():,.0f} m") name_left fclass name_right distance_m
Belmesthorpe village Kendrew Barracks Headquarters 13943.0
Borderville hamlet Kendrew Barracks Headquarters 13883.0
Essendine village Kendrew Barracks Headquarters 13648.0
Ryhall village Kendrew Barracks Headquarters 12992.0
Tinwell village Kendrew Barracks Headquarters 12780.0
places served by each supermarket: {'Co-op Food': 25, 'Kendrew Barracks Headquarters': 21, 'The Co-operative Food': 9, 'Aldi': 5, 'Tesco': 4}
without it: 21 places re-matched; their median distance 8,167 m -> 12,143 m
Every one of the five places furthest from a supermarket is matched to Kendrew Barracks Headquarters, and so are 21 of the 64 places. Whether that feature is a shop the public can use is exactly the kind of question the join cannot answer; the last line shows how much the answer depends on it.
Example 2: The same join in degrees
import warnings
with warnings.catch_warnings(record=True) as caught:
warnings.simplefilter("always")
in_degrees = places.to_crs(4326).sjoin_nearest(supermarkets.to_crs(4326)[["name", "geometry"]],
distance_col="distance")
print(str(caught[0].message)[:75], "…")
print(f"median 'distance': {in_degrees['distance'].median():.4f} (degrees)")
same = (in_degrees["name_right"].values == nearest["name_right"].values).mean()
print(f"same nearest supermarket as in metres: {same:.0%} of places")Geometry is in a geographic CRS. Results from 'sjoin_nearest' are likely in …
median 'distance': 0.0550 (degrees)
same nearest supermarket as in metres: 84% of places
The distances are in degrees and meaningless as metres. At Rutland's latitude a degree of longitude is about 60% of a degree of latitude, so "nearest" in degrees can pick a different feature than nearest on the ground — and here it did, for about one place in six. The warning says the results are "likely incorrect"; the count says which ones.
Example 3: Find the nearest point to one single geometry
For a single point, direct distance calculation is simple.
site = places_bng.loc[places_bng["name"] == "Uppingham", "geometry"].iloc[0]
candidates = supermarkets.copy()
candidates["distance_m"] = candidates.geometry.distance(site)
nearest_one = candidates.loc[candidates["distance_m"].idxmin()]
print(nearest_one[["name", "distance_m"]].to_dict()){'name': 'Co-op Food', 'distance_m': 179.0182651093306}
This works well for one geometry or a small number of checks.
A nearest join trusts its inputs completely
sjoin_nearest builds a spatial index over the target layer and, for each source feature, asks it for the closest geometry. That is far faster than a loop of distance calculations, and it is exact — for the coordinates and the features it is given. Both of this page's surprises are about what it is given.
Distance is computed in the layer's units. In longitude and latitude a degree east is shorter than a degree north — at Rutland's latitude about 60% as long — so "nearest" in degrees favours targets to the east or west, and for 16% of Rutland's places it named a different supermarket from the one nearest on the ground. And the join has no idea what a feature is. Kendrew Barracks Headquarters is tagged as a supermarket in OpenStreetMap, so for 21 places it is the nearest supermarket; remove it and their median distance goes from 8.2 km to 12.1 km. The nearest feature is always found; whether it answers the question depends on the layer. How ties and search radius work in general is in nearest-neighbour joins explained.
Polygons as targets, ties and invalid input
- Polygon targets measure to the edge. Six of the seven supermarkets are building outlines; this page reduced each to a point inside it with
representative_point(), so distances run to a point inside each building rather than to its nearest wall. - Ties give extra rows. None here; with a regular grid of targets, check
nearest.index.duplicated().sum().
Check for null or invalid geometry values before running the join.
places_bng = places_bng[places_bng.geometry.notna() & places_bng.is_valid].copy()
supermarkets = supermarkets[supermarkets.geometry.notna() & supermarkets.is_valid].copy()
print(len(places_bng), "places,", len(supermarkets), "supermarkets after the check")64 places, 7 supermarkets after the check
Internal links
If you are not sure which CRS to use, see the broader concept in How to Reproject Spatial Data in Python (GeoPandas).
Related tasks:
If your output distances look wrong, first check your reprojection workflow in How to Reproject Spatial Data in Python (GeoPandas).
FAQ
How do I find the nearest point in GeoPandas?
Use sjoin_nearest() to match each feature in one GeoDataFrame to the closest feature in another:
result = source.sjoin_nearest(target, how="left", distance_col="distance_m")Should I use sjoin_nearest() or calculate distance manually?
Use sjoin_nearest() for layer-to-layer matching. Use manual distance calculation when you only need the nearest geometry for one feature or a very small number of features.
Why are my nearest distance values incorrect?
The most common reason is CRS. If your data is in latitude/longitude, distances are calculated in degrees. Reproject to a local projected CRS before running the analysis.
Can I set a maximum distance for nearest matches?
Yes. Use max_distance in sjoin_nearest():
result = source.sjoin_nearest(
target,
max_distance=5000,
distance_col="distance_m"
)How do I keep the geometry of the matched nearest point, not just its attributes?
sjoin_nearest keeps the left geometry, so rename and carry the target geometry before the join, for example stores["store_geom"] = stores.geometry. That column survives the join and lets you draw connecting lines or measure later.
Why does sjoin_nearest return more rows than my source layer?
This happens when two or more target features are tied at the same nearest distance, producing duplicate rows for one source point. Drop the extras with nearest = nearest[~nearest.index.duplicated(keep="first")] if you need exactly one match per source feature.
How do I exclude a point from matching itself when both layers are the same?
Run the join, then filter out self-matches where the indices line up, for example by comparing an ID column and keeping rows where left_id != right_id. A plain sjoin_nearest of a layer against itself returns each point matched to itself at distance zero.
How do I find the k nearest points instead of just the single closest one?
sjoin_nearest accepts exclusive and a max_distance, but for k-nearest use scipy.spatial.cKDTree on the coordinate arrays from points_from_xy, querying with k=3. Build the tree in a projected CRS so the returned distances are in metres.