Spatial Join Returns Empty Results in GeoPandas: How to Fix It

Problem statement

Join Rutland's food businesses to the wards they sit in, and geopandas.sjoin() returns 0 rows with a single warning: the Food Standards Agency's register arrives in longitude and latitude, the ONS wards in British National Grid. Relabel the points with set_crs and the warning goes away while the result stays empty. Only a reprojection matches them: every business with a location gets its ward, and the 124 of 462 left over have no coordinates at all. The same empty result also comes from layers that do not overlap, a predicate stricter than the question, unusable geometries, and an inner join that drops what it could not match; each is checked below on the same two layers.

Quick answer

If geopandas.sjoin() returns no matches:

1. check that both GeoDataFrames have the same CRS with .crs
2. reproject with .to_crs(), not .set_crs()
3. compare .total_bounds to confirm the layers overlap
4. test a different predicate such as intersects
5. inspect null, empty, and invalid geometries before joining

Diagnosis at a glance

Working through an empty join

Check whether both GeoDataFrames use the same CRS

A CRS mismatch is one of the most common reasons a spatial join returns no rows.

Compare gdf.crs values before running sjoin

print("Points CRS:  ", points.crs)
print("Polygons CRS:", polygons.crs)
with warnings.catch_warnings(record=True) as caught:
    warnings.simplefilter("always")
    print("sjoin as loaded:", len(gpd.sjoin(points, polygons, predicate="within")), "rows |",
          str(caught[0].message).split(".")[0] if caught else "no warning")

Points CRS:   EPSG:4326
Polygons CRS: EPSG:27700
sjoin as loaded: 0 rows | CRS mismatch between the CRS of left geometries and the CRS of right geometries

One layer is in EPSG:4326 and the other in EPSG:27700, so the coordinates are not in the same system: the join ran, warned once, and returned nothing.

Reproject one layer with to_crs() instead of assigning CRS manually

Use to_crs() when the layer already has the correct CRS definition but needs to be converted.

if points.crs != polygons.crs:
    points_bng = points.to_crs(polygons.crs)

joined = gpd.sjoin(points_bng, polygons, how="inner", predicate="intersects")
print(len(joined), "matches;", joined["WD24NM"].value_counts().head(3).to_dict())

338 matches; {'Oakham North East': 73, 'Uppingham': 66, 'Normanton': 28}

Common mistake: confusing set_crs() with to_crs()

set_crs() only labels the data. It does not transform coordinates.

relabelled = points.set_crs(polygons.crs, allow_override=True)      # wrong: labels degrees as metres
print("set_crs:", len(gpd.sjoin(relabelled, polygons, predicate="within")), "matches |",
      "to_crs:", len(gpd.sjoin(points.to_crs(polygons.crs), polygons, predicate="within")), "matches")

set_crs: 0 matches | to_crs: 338 matches

set_crs() made the warning go away and the result stayed empty: the points now claim to be a few metres from the grid's origin, in the sea off the Isles of Scilly.

Use set_crs() only when the data has no CRS assigned and you know what the original CRS should be.

Confirm the layers actually overlap in space

Even if both layers have the same CRS, they may still be in different places.

Inspect total_bounds for both layers

print("Points bounds:  ", points.total_bounds.round(3))
print("Polygons bounds:", polygons.total_bounds.round(0))

Points bounds:   [-0.796 52.533 -0.454 52.75 ]
Polygons bounds: [479688. 292648. 506296. 318928.]

total_bounds returns [minx, miny, maxx, maxy]. Values between −1 and 53 against values in the hundreds of thousands are degrees against metres: two coordinate systems, whatever the metadata says.

Plot both layers together to verify overlap

A visual check is often the fastest way to confirm whether a join should produce matches.

ax = polygons.plot(facecolor="none", edgecolor="black", figsize=(8, 8))
points_bng.plot(ax=ax, color="red", markersize=10)
plt.close("all")

If the points are not visible or appear far away from the polygons, the join will return no matches.

Watch for wrong CRS metadata and bad source files

This often happens when one file is longitude/latitude and the other is projected in meters. It can also happen when a shapefile has a wrong .prj file or missing CRS metadata.

Use the correct spatial predicate

The predicate controls what counts as a match. If you choose the wrong one, the join can return zero rows even when the layers are close or overlapping.

Start with intersects when debugging

intersects is usually the safest first test because it matches geometries that overlap, cross, touch, or fall within one another.

joined = gpd.sjoin(points_bng, polygons, how="inner", predicate="intersects")
print(len(joined))

338

Be careful with within, contains, and direction

For a point-in-polygon join:

• points are usually within polygons
• polygons usually contain points

Direction matters. This returns no matches:

print(len(gpd.sjoin(points_bng, polygons, how="inner", predicate="contains")))

0

That asks whether each point contains a polygon, which is false.

Test a few predicates quickly

for pred in ["intersects", "within", "contains", "touches"]:
    try:
        result = gpd.sjoin(points_bng, polygons, how="inner", predicate=pred)
        print(pred, len(result))
    except Exception as e:
        print(pred, "not available:", e)

intersects 338
within 338
contains 0
touches 0

Check available predicates in your environment

Available predicates can vary by GeoPandas, Shapely, and the spatial index backend.

print(sorted(p for p in points_bng.sindex.valid_query_predicates if p))

['contains', 'contains_properly', 'covered_by', 'covers', 'crosses', 'dwithin', 'intersects', 'overlaps', 'touches', 'within']

If a predicate is unsupported in your environment, GeoPandas raises an error. Check the available predicates before running the join.

Check for empty, null, or invalid geometries

Bad geometry can prevent expected matches.

Find null geometries with isna()

print("Null point geometries:", points.geometry.isna().sum())
print("Null polygon geometries:", polygons.geometry.isna().sum())

Null point geometries: 124
Null polygon geometries: 0

The register has businesses with no location, and a null geometry can never match: an inner join silently drops them.

Find empty geometries with is_empty

print("Empty point geometries:", points.geometry.is_empty.sum())
print("Empty polygon geometries:", polygons.geometry.is_empty.sum())

Empty point geometries: 0
Empty polygon geometries: 0

Detect invalid geometries with is_valid

print("Invalid polygons:", (~polygons.geometry.is_valid).sum())

Invalid polygons: 0

Repair invalid polygons before joining

The ONS wards are valid. If yours are not, repair them with make_valid() and rerun the join:

repaired = polygons.copy()
repaired["geometry"] = repaired.geometry.make_valid()
print("Invalid polygons after repair:", (~repaired.geometry.is_valid).sum(), "|",
      len(gpd.sjoin(points_bng, repaired, predicate="intersects")), "matches")

Invalid polygons after repair: 0 | 338 matches

buffer(0) also repairs some self-intersections, but it can silently drop part of a figure-eight polygon; prefer make_valid().

Make sure the geometry types match the task

The geometry types in each layer affect which predicates make sense.

Point-in-polygon joins

If your goal is assigning points to polygons, the usual pattern is:

joined = gpd.sjoin(points_bng, polygons, how="left", predicate="within")
print(len(points_bng), "points ->", len(joined), "rows")

462 points -> 462 rows

Polygon-to-polygon joins

For polygon overlays, intersects is often more useful than within or contains, especially when boundaries only partially overlap.

Line and polygon cases

If one layer contains lines and the other polygons, within may be too strict. If geometries only touch at boundaries, contains may return no matches while intersects or touches returns some.

Verify the join direction and join type

The join may not be empty. It may only be hiding unmatched rows.

Difference between left, right, and inner

• inner: returns only matched rows
• left: keeps all rows from the left GeoDataFrame
• right: keeps all rows from the right GeoDataFrame

If you use inner, unmatched features disappear completely, which makes debugging harder.

Keep the layer you want to inspect on the left

Use the dataset you want to preserve as the left table — joined from the point-in-polygon section above.

Inspect unmatched rows after the join

unmatched = joined[joined["index_right"].isna()]
print("Unmatched rows:", len(unmatched), "| of which have no geometry:", int(unmatched.geometry.isna().sum()))
print(unmatched["BusinessType"].value_counts().head(3).to_dict())

Unmatched rows: 124 | of which have no geometry: 124
{'Other catering premises': 53, 'Mobile caterer': 18, 'Retailers - other': 17}

Every unmatched business is one without a location — catering premises, mobile caterers and others the register does not place on a map — not a join failure. That is what a left join shows and an inner join hides.

Run a minimal test to isolate the issue

A small test can help confirm whether the problem is with the data or with your join logic.

Test with a small subset of rows

points_sample = points_bng.dropna(subset=["geometry"]).iloc[:5].copy()
polygons_sample = polygons.iloc[:5].copy()
test_join = gpd.sjoin(points_sample, polygons_sample, how="inner", predicate="intersects")
print(len(test_join), "matches from the first 5 of each")

1 matches from the first 5 of each

This only helps if the sample features are expected to overlap. A result from arbitrary first rows proves nothing either way: five businesses and five wards chosen by row order may simply be in different parts of the county.

Compare one geometry against another directly

pt = points_bng.geometry.dropna().iloc[0]
poly = polygons.geometry[polygons.contains(pt)].iloc[0]
print("Intersects:", pt.intersects(poly), "| Within:", pt.within(poly), "| Touches:", pt.touches(poly))

Intersects: True | Within: True | Touches: False

If direct geometry tests also return False, the problem is usually CRS, bad geometry, the wrong predicate, or no real overlap.

Debugging Rutland hygiene ratings against wards

The data used on this page

The Food Standards Agency's food hygiene rating data for Rutland (OGL), with points built from its longitude and latitude, and Rutland's wards from the ONS December 2024 boundaries (OGL) in British National Grid.

import json
import warnings
from pathlib import Path

import geopandas as gpd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import pandas as pd
import requests

Path("data").mkdir(exist_ok=True)
if not Path("data/food.json").exists():
    # the FSA register file as saved on 9 October 2026 (updated daily); live:
    # ratings.food.gov.uk/api/open-data-files/FHRS879en-GB.json
    Path("data/food.json").write_bytes(requests.get("https://data.spatialworkflow.io/fhrs-879-2026-10-09.json", timeout=120).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")
points = gpd.GeoDataFrame(records[["FHRSID", "BusinessName", "BusinessType"]],
                          geometry=gpd.GeoSeries(gpd.points_from_xy(x, y)).where(x.notna()), crs=4326)
polygons = gpd.read_file("data/wards.geojson")
print(len(points), "businesses,", len(polygons), "wards")

462 businesses, 15 wards

Example: debug a point-to-polygon join

The whole sequence as one function, run on the layers exactly as they arrive.

def debug_join(points, polygons):
    print("CRS:", points.crs.to_epsg(), polygons.crs.to_epsg())                  # 1. check CRS
    if points.crs != polygons.crs:
        points = points.to_crs(polygons.crs)
    print("bounds:", points.total_bounds.round(0), polygons.total_bounds.round(0))  # 2. check extents
    print("null points:", int(points.geometry.isna().sum()), "| empty:", int(points.geometry.is_empty.sum()),
          "| invalid polygons:", int((~polygons.geometry.is_valid).sum()))       # 3. geometry quality
    joined = gpd.sjoin(points, polygons, how="left", predicate="intersects")      # 4. a safe first join
    unmatched = joined[joined["index_right"].isna()]                              # 5. inspect unmatched
    print(f"{len(joined)} rows, {len(unmatched)} unmatched")
    return joined

joined = debug_join(points, polygons)

CRS: 4326 27700
bounds: [481609. 293566. 504554. 317703.] [479688. 292648. 506296. 318928.]
null points: 124 | empty: 0 | invalid polygons: 0
462 rows, 124 unmatched

Four things a join needs to match

GeoPandas narrows the candidates with a spatial index and then tests each with the predicate. Degrees against metres fail at the first stage: the index finds no candidate pairs, so there is nothing for the predicate to reject, and no error to raise.

Border-only contact, precision and wrong CRS metadata

• Layers that only touch at borders: within and contains may return no matches, while intersects or touches may work.
• Multipart geometries: MultiPolygons can behave differently than expected when testing containment near holes or gaps.
• Sampling for tests: head() or iloc[:5] is only useful if those features are expected to overlap.

Internal links

If you need the underlying concept, read How Spatial Join Works in GeoPandas.

For the standard workflow, see Perform a Spatial Join in GeoPandas and Reproject a GeoDataFrame in GeoPandas.

If the problem looks projection-related, use GeoPandas CRS Mismatch: How to Fix Projection Problems.

You may also need How to Fix CRS Mismatch in GeoPandas or GeoPandas Not Reading Shapefile: Common Causes and Fixes if the issue starts earlier in the workflow.

FAQ

Why does sjoin() return an empty GeoDataFrame without an error?

Because the function can run correctly even when no features satisfy the spatial relationship. The most common causes are CRS mismatch, no spatial overlap, the wrong predicate, or bad geometries.

How do I check whether two GeoDataFrames actually overlap?

Compare .total_bounds for both layers and plot them together on the same axis. If the extents do not overlap visually or numerically, the join will return no matches.

Can a missing spatial index cause sjoin() to return nothing?

No. GeoPandas builds the spatial index automatically, and a missing index would not silently drop matches. Empty results come from CRS mismatch, no real overlap, the wrong predicate, or invalid geometries, not from the index itself.

Could a predicate="op" argument be why my join fails?

The old op= keyword was removed; current GeoPandas uses predicate=. If you see an error or unexpected behavior from an op= argument, update the call to gpd.sjoin(points, polygons, predicate="intersects").

After fixing the CRS, why are only some features matching?

Partial matches usually mean the geometries genuinely overlap only in part, or some rows have null, empty, or invalid geometry that fails the predicate. Check the unmatched rows with a left join and inspect their geometry with is_empty and is_valid before assuming the join logic is wrong.