How to Perform a Spatial Join in Python (GeoPandas)
On this page
Problem statement
A spatial join lets you transfer attributes from one spatial layer to another based on location instead of a shared ID field.
Typical GIS tasks include:
- assigning district names to customer points
- joining addresses to census tracts
- attaching flood zone attributes to parcel polygons
- matching building footprints to administrative boundaries
In GeoPandas, this means taking two vector layers and joining them using a spatial relationship such as intersects, within, or contains. This page does it with real layers for Rutland: the Food Standards Agency's register of food businesses as points, the county's 15 wards from the ONS, its Land Registry parcels, and OpenStreetMap roads, pubs and land use.
Quick answer
Use geopandas.sjoin() with two GeoDataFrames that use the same CRS, then choose the correct spatial predicate for your workflow:
points = businesses
polygons = wards
if points.crs != polygons.crs:
polygons = polygons.to_crs(points.crs)
joined = gpd.sjoin(points, polygons, how="left", predicate="within")
print(joined[["BusinessName", "WD24NM"]].head(3).to_string())
print(len(points), "points in,", len(joined), "rows out") BusinessName WD24NM
0 10 Dining Street Oakham North East
1 5 Elements Cafe and Lounge Uppingham
2 A Simpler Life Uppingham
462 points in, 462 rows out
This keeps all points and adds polygon attributes where each point falls within a polygon.
Geometry types, one CRS, a predicate, the row count
Load the two spatial layers and check geometry types
Geometry type affects which spatial predicate makes sense:
for name, layer in (("businesses", businesses), ("wards", wards), ("parcels", parcels), ("roads", roads)):
print(f"{name:10} {type(layer).__name__:12} {layer.geom_type.value_counts().to_dict()}")businesses GeoDataFrame {'Point': 338}
wards GeoDataFrame {'Polygon': 15}
parcels GeoDataFrame {'Polygon': 21253}
roads GeoDataFrame {'LineString': 7292}
Typical join patterns are points to polygons, polygons to polygons, and lines to polygons. If you want to assign polygon attributes to points, within is usually the clearest choice.
Make sure both layers use the same CRS
Spatial joins require both layers to use the same coordinate reference system:
print(businesses.crs.to_epsg(), wards.crs.to_epsg())
wards_wgs84 = wards.to_crs(businesses.crs)
print("after to_crs:", wards_wgs84.crs.to_epsg())4326 27700
after to_crs: 4326
If you skip this step, GeoPandas warns about the mismatch and the join returns no matches — see spatial join returns empty results.
Choose the right spatial predicate
Common predicates include:
intersects: features touch or overlap in any waywithin: left geometry is completely inside right geometrycontains: left geometry completely contains right geometry
The ones your installation supports come from the spatial index:
print(sorted(p for p in businesses.sindex.valid_query_predicates if p))['contains', 'contains_properly', 'covered_by', 'covers', 'crosses', 'dwithin', 'intersects', 'overlaps', 'touches', 'within']
Common choices: points inside polygons, usually within; overlapping polygons, usually intersects; lines crossing polygons, usually intersects.
Run the spatial join, and read the row count
for how in ("inner", "left", "right"):
result = gpd.sjoin(businesses, wards_wgs84, how=how, predicate="within")
print(f"how={how:5}: {len(result)} rows, geometry from the {'right' if how == 'right' else 'left'} layer")how=inner: 338 rows, geometry from the left layer
how=left : 462 rows, geometry from the left layer
how=right: 338 rows, geometry from the right layer
inner: keep only matched rows — here the 124 businesses without a location disappearleft: keep all rows from the left GeoDataFrame, unmatched ones withNaNin the joined columnsright: keep all rows from the right one, and take the right layer's geometry — a row per match, plus one for any ward with no business (every Rutland ward has at least one)
In most point-in-polygon workflows, the left layer is the layer you want to preserve.
Review the output columns
The result includes the left layer's columns, the right layer's columns (unless you subset the right GeoDataFrame first) and an index_right column. When both layers have a column with the same name, GeoPandas adds suffixes:
print(joined.columns.tolist()[-4:])
clash = gpd.sjoin(pubs, landuse[["name", "fclass", "geometry"]], how="inner", predicate="within")
print([c for c in clash.columns if c.startswith(("name", "fclass"))], "|", clash[["name_left", "fclass_right"]].head(2).values.tolist())['BusinessType', 'geometry', 'index_right', 'WD24NM']
['fclass_left', 'name_left', 'name_right', 'fclass_right'] | [['George and Dragon', 'residential'], ['Exeter Arms', 'residential']]
Save the joined result
Path("output").mkdir(exist_ok=True)
joined.to_file("output/businesses_with_wards.gpkg", driver="GPKG")
print(gpd.read_file("output/businesses_with_wards.gpkg").columns.tolist()[:4], "…")['FHRSID', 'BusinessName', 'BusinessType', 'index_right'] …
Write GeoPackage rather than shapefile: the joined column names here are long, and a shapefile would cut them to ten characters.
Points, overlapping polygons and lines joined to Rutland wards
The data used on this page
The Food Standards Agency's food hygiene rating data for Rutland (OGL); Rutland's wards from the ONS December 2024 boundaries (OGL); HM Land Registry's INSPIRE parcels (© Crown copyright and database rights; OGL); and roads, pubs and land use from Geofabrik's OpenStreetMap extract of 1 January 2026 (© OpenStreetMap contributors, ODbL). One trap from writing this setup: a where= filter only sees the columns you select, so where="fclass = 'pub'" with columns=["name"] returned no rows at all, and no error — include the filtered column.
import json
from pathlib import Path
import geopandas as gpd
import pandas as pd
import requests
Path("data").mkdir(exist_ok=True)
downloads = {
# the FSA register file for Rutland as saved on 9 October 2026 (updated daily); live:
# ratings.food.gov.uk/api/open-data-files/FHRS879en-GB.json
"data/food.json": "https://data.spatialworkflow.io/fhrs-879-2026-10-09.json",
"data/rutland.shp.zip": "https://download.geofabrik.de/europe/united-kingdom/england/rutland-260101-free.shp.zip",
"data/Rutland.zip": "https://use-land-property-data.service.gov.uk/datasets/inspire/download/Rutland.zip",
}
for target, url in downloads.items():
if not Path(target).exists():
Path(target).write_bytes(requests.get(url, timeout=300).content)
if not Path("data/wards.geojson").exists():
Path("data/wards.geojson").write_bytes(requests.get(
"https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
"Wards_December_2024_Boundaries_UK_BFC/FeatureServer/0/query",
params={"where": "LAD24CD='E06000017'", "outFields": "WD24NM", "outSR": 27700, "f": "geojson"}, timeout=120).content)
found = json.loads(Path("data/food.json").read_text())["FHRSEstablishment"]["EstablishmentCollection"]
records = pd.json_normalize(found if isinstance(found, list) else found["EstablishmentDetail"])
x = pd.to_numeric(records["Geocode.Longitude"], errors="coerce")
y = pd.to_numeric(records["Geocode.Latitude"], errors="coerce")
businesses = gpd.GeoDataFrame(records[["FHRSID", "BusinessName", "BusinessType"]],
geometry=gpd.GeoSeries(gpd.points_from_xy(x, y)).where(x.notna()), crs=4326)
wards = gpd.read_file("data/wards.geojson")
parcels = gpd.read_file("/vsizip/data/Rutland.zip/Land_Registry_Cadastral_Parcels.gml", columns=["INSPIREID"], DOWNLOAD_SCHEMA="NO")
OSM = "/vsizip/data/rutland.shp.zip/"
roads = gpd.read_file(OSM + "gis_osm_roads_free_1.shp", columns=["name", "fclass"]).to_crs(27700)
pubs = gpd.read_file(OSM + "gis_osm_pois_free_1.shp", where="fclass = 'pub'", columns=["name", "fclass"])
landuse = gpd.read_file(OSM + "gis_osm_landuse_a_free_1.shp", columns=["name", "fclass"])
print(len(businesses), "businesses,", len(wards), "wards,", len(parcels), "parcels,", len(roads), "road segments")462 businesses, 15 wards, 21253 parcels, 7292 road segments
Example 1: Count points per polygon
per_ward = (gpd.sjoin(businesses, wards_wgs84[["WD24NM", "geometry"]], how="inner", predicate="within")
.groupby("WD24NM").size().sort_values(ascending=False))
print(per_ward.head(5).to_string())
print("wards with no located business:", sorted(set(wards["WD24NM"]) - set(per_ward.index)))WD24NM
Oakham North East 73
Uppingham 66
Normanton 28
Exton 24
Braunston & Martinsthorpe 23
wards with no located business: []
Example 2: Join overlapping polygons — and the duplicate rows
Attach each parcel to the ward it lies in:
parcel_ward = gpd.sjoin(parcels, wards[["WD24NM", "geometry"]], how="inner", predicate="intersects")
repeats = parcel_ward.index.value_counts()
print(f"{len(parcels):,} parcels -> {len(parcel_ward):,} rows; {int((repeats > 1).sum())} parcels matched more than one ward")
print("most wards for one parcel:", int(repeats.max()))21,253 parcels -> 21,843 rows; 569 parcels matched more than one ward
most wards for one parcel: 4
One parcel produces a row per ward it touches. If you need exactly one ward per parcel, decide how: by the ward containing the parcel's representative point, or by the largest overlap:
points_in = parcels.assign(geometry=parcels.representative_point())
one_each = gpd.sjoin(points_in, wards[["WD24NM", "geometry"]], how="left", predicate="within")
print(f"representative point: {len(one_each):,} rows for {len(parcels):,} parcels; unmatched {int(one_each['WD24NM'].isna().sum())}")representative point: 21,253 rows for 21,253 parcels; unmatched 246
One row per parcel. The unmatched ones lie outside every Rutland ward: the Land Registry's county file includes parcels that straddle or sit just beyond the boundary.
Example 3: Lines to polygons
road_ward = gpd.sjoin(roads, wards[["WD24NM", "geometry"]], how="inner", predicate="intersects")
crossing = road_ward.index.value_counts()
print(f"{len(roads):,} segments -> {len(road_ward):,} rows; {int((crossing > 1).sum())} segments cross a ward boundary")7,292 segments -> 7,222 rows; 258 segments cross a ward boundary
A road segment that crosses a boundary belongs to both wards. For lengths per ward, clip the roads with gpd.overlay(roads, wards, how="intersection") rather than joining them, so each piece is counted once.
Example 4: Reproject before joining
sites = pubs # EPSG:4326
boundaries = wards # EPSG:27700
print(sites.crs.to_epsg(), boundaries.crs.to_epsg())
if sites.crs != boundaries.crs:
boundaries = boundaries.to_crs(sites.crs)
joined_pubs = gpd.sjoin(sites, boundaries[["WD24NM", "geometry"]], how="left", predicate="within")
print(joined_pubs[["name", "WD24NM"]].dropna().head(3).to_string(index=False))4326 27700
name WD24NM
George and Dragon Lyddington
Exeter Arms Ketton
The Northwick Arms Ketton
Geometry relationships instead of a key
A spatial join is different from a normal attribute join.
- An attribute join uses a shared key such as
parcel_id - A spatial join uses geometry relationships such as overlap, containment, or intersection
In GeoPandas, sjoin() compares the geometry in the left GeoDataFrame with the geometry in the right GeoDataFrame using the predicate you choose, using a spatial index to cut the number of comparisons.
A spatial join can also increase row counts. One feature may match multiple features, especially when polygons overlap a boundary, a line crosses several polygons, or the right-hand layer is not cleanly non-overlapping. That is why the output may contain more rows than the left input layer — as the parcels and roads above did.
Duplicate matches and boundary behaviour
Duplicate matches
One feature can match more than one feature.
This often happens with:
- overlapping polygon layers
- polygon-to-polygon joins
- line-to-polygon joins
- boundary cases depending on the predicate
If you need exactly one match per feature, you will need additional cleanup logic after the join.
CRS mismatches
Always check CRS before every spatial join.
print(left_gdf.crs)
print(right_gdf.crs)If needed:
right_gdf = right_gdf.to_crs(left_gdf.crs)Missing or invalid geometries
Null or invalid geometries can cause missing matches or errors.
left_gdf = left_gdf[left_gdf.geometry.notna()]
right_gdf = right_gdf[right_gdf.geometry.notna()]
right_gdf = right_gdf[right_gdf.is_valid]In some workflows, you may need to repair invalid geometry before joining.
Boundary behavior
Points or lines on polygon boundaries can produce unexpected results.
withinrequires the left geometry to be in the interior of the right geometrycontainsis the mirror image: the right geometry must lie in the left one's interior, with boundary-only contact not enoughintersectsis more inclusive and matches touching geometries
If features on edges are not matching as expected, test intersects and compare the output.
Performance on large datasets
Large joins can be slow.
To reduce overhead:
- keep only the columns you need before joining
- make sure a spatial index is available
- filter either dataset first if you can reduce the search area
- avoid unnecessary reprojection or repeated file reads inside loops
Internal links
- How to Read a Shapefile in Python with GeoPandas
- How to Export GeoJSON in Python with GeoPandas
- How to Reproject Spatial Data in Python (GeoPandas)
FAQ
How do I do a spatial join in GeoPandas?
Use gpd.sjoin(left_gdf, right_gdf, how=..., predicate=...). Both inputs must be GeoDataFrames and both must use the same CRS.
joined = gpd.sjoin(points, polygons, how="left", predicate="within")What is the difference between intersects, within, and contains in a spatial join?
intersects: geometries touch or overlap in any waywithin: left geometry is inside the right geometrycontains: left geometry contains the right geometry
For joining points to polygons, within is usually the best starting point.
Why does my spatial join return empty results?
Common causes include:
- CRS mismatch
- wrong predicate
- invalid or null geometries
- empty input layers
Check CRS first, then inspect geometry validity and test a different predicate.
Why do I get multiple rows for one feature after a spatial join?
Because one feature can match multiple features in the other layer. This is common with overlapping polygons and polygon intersection joins.
How do I count how many points fall inside each polygon with a spatial join?
Run the join with points on the left, then group by the polygon's index column and count:
joined = gpd.sjoin(points, polygons, how="inner", predicate="within")
counts = joined.groupby("index_right").size()You can also use gpd.sjoin together with a follow-up groupby, or reach for gdf.sjoin(...).value_counts() on the matched key.
What is the difference between sjoin and sjoin_nearest?
sjoin matches features by a topological predicate such as within or intersects, so a point that falls outside every polygon gets no match. sjoin_nearest instead matches each left feature to the closest right feature by distance, which is useful for snapping points to the nearest road or facility.
How do I keep only one match per feature after a spatial join produces duplicates?
After the join, drop duplicates on the left layer's identifier, for example joined.drop_duplicates(subset="customer_id"). To pick a specific match rather than an arbitrary one, sort first by a distance or priority column, then drop duplicates so the preferred row is kept.
Can I run a spatial join on geographic coordinates in EPSG:4326?
Yes. Predicate-based joins like within and intersects are topological, so they work correctly in EPSG:4326 as long as both layers share that CRS. Reproject to a projected CRS only when you need distance-based work such as sjoin_nearest with a max_distance in meters.