Spatial Predicates Explained: Intersects, Within, Contains and the Rest

Problem statement

sjoin takes a predicate argument. Most people leave it at the default, discover the row count is wrong, and then try each option until the number looks plausible.

for predicate in ("intersects", "within", "contains"):
    print(f"{predicate:10}", len(gpd.sjoin(points, wards, predicate=predicate)), "rows")

intersects 432 rows
within     0 rows
contains   0 rows

These are points on Rutland's ward boundaries — the vertices of the Oakham North East ward, the kind of points a GPS track along a boundary road produces. Three predicates, two answers — every point, or none — and no obvious way to tell which one is correct — because correctness depends on a question nobody wrote down: what should happen to a point that falls exactly on a ward boundary?

Predicates are not interchangeable filters with different strictness. Each one asks a precise question about how two geometries relate, and the differences between them are almost entirely about boundaries — the edges that two shapes share.

Quick answer

Every predicate is a question about interiors and boundaries:

from shapely.geometry import Point, Polygon

square = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])
inside = Point(1, 1)
on_edge = Point(0, 1)      # exactly on the left boundary
outside = Point(3, 3)

for name, p in [("inside", inside), ("on_edge", on_edge), ("outside", outside)]:
    print(f"{name:8s} intersects={square.intersects(p)!s:5}  contains={square.contains(p)!s:5}  "
          f"covers={square.covers(p)!s:5}  within={p.within(square)!s:5}  touches={p.touches(square)}")

inside   intersects=True   contains=True   covers=True   within=True   touches=False
on_edge  intersects=True   contains=False  covers=True   within=False  touches=True
outside  intersects=False  contains=False  covers=False  within=False  touches=False

That middle row is the whole subject. contains is False for a point on the boundary; covers is True. Every "why did my count change" question about predicates traces back to it.

The relationships, drawn

Interior, boundary, exterior, and each predicate

The mental model: interior, boundary, exterior

Every geometry divides space into three parts. A polygon has an interior (the area inside), a boundary (its ring), and an exterior (everything else). A line has an interior (the line minus its endpoints), a boundary (the two endpoints), and an exterior. A point has an interior (itself), an empty boundary, and an exterior.

Every predicate is a statement about which of those parts intersect. That is not an analogy — it is literally how they are defined, in a 3×3 grid called the DE-9IM matrix:

print(square.relate(on_edge), "|", on_edge.relate(square))

FF20F1FF2 | F0FFFF212

Nine characters, one per interior/boundary/exterior pairing. You will rarely read one directly, but knowing they exist explains why the predicates behave so precisely: they are not heuristics, they are named patterns over that matrix.

intersects is the one to use unless you have a reason

intersects is true if the shapes share any point — interior, boundary, or a single corner. It is the most permissive, and the right default for "find things near/on/in each other".

print(len(gpd.sjoin(businesses, wards, predicate="intersects")), "of", len(businesses), "business points matched")

338 of 338 business points matched

Use something stricter only when you can say what you are excluding and why.

within and contains are the same question, reversed

a, b = businesses.geometry.iloc[0], wards.geometry[wards.contains(businesses.geometry.iloc[0])].iloc[0]
print(a.within(b), b.contains(a))      # is A inside B? is B holding A? — identical truth value

True True

In a spatial join, the difference is which frame is left and which is right:

per_point = gpd.sjoin(businesses, wards, predicate="within")       # each point tagged with its ward
per_ward = gpd.sjoin(wards, businesses, predicate="contains")      # each ward tagged with its points
print(len(per_point), "rows with point geometry |", len(per_ward), "rows with ward geometry")

338 rows with point geometry | 338 rows with ward geometry

Same relationship, different row shape. Pick based on what you want one row per.

The catch: a polygon does not contain its own boundary, so a point exactly on the edge is not "within" the ward — the problem statement's zero. Why the definition is written that way is in the explanation below.

covers is the fix for boundary cases

print(square.contains(on_edge), square.covers(on_edge))   # no interior intersection | no point of A outside B

False True

covers is contains without the interior requirement. For administrative geography — points on a road that is the boundary, parcels sharing an edge with the ward outline — covers is usually what people mean when they say "inside".

GeoPandas supports it in joins from version 0.10:

print(len(gpd.sjoin(wards, points, predicate="covers")), "ward-point pairs with covers, against",
      len(gpd.sjoin(wards, points, predicate="contains")), "with contains")

432 ward-point pairs with covers, against 0 with contains

touches and overlaps are mutually exclusive

ward_pairs = gpd.sjoin(wards, wards, predicate="intersects")
ward_pairs = ward_pairs[ward_pairs.index < ward_pairs["index_right"]]
left, right = wards.geometry.loc[ward_pairs.index].values, wards.geometry.loc[ward_pairs["index_right"]].values
print(f"neighbouring ward pairs: {len(ward_pairs)}, touches: {int(shapely.touches(left, right).sum())}, "
      f"overlaps: {int(shapely.overlaps(left, right).sum())}")

parcel_pairs = gpd.sjoin(parcels, parcels, predicate="overlaps")
print(f"Land Registry parcel pairs that overlap: {int((parcel_pairs.index < parcel_pairs['index_right']).sum()):,}")

neighbouring ward pairs: 31, touches: 31, overlaps: 0
Land Registry parcel pairs that overlap: 11,453

Wards are drawn as a single coverage by one organisation: every neighbouring pair touches and none overlaps. Two adjacent parcels sharing a fence line should behave the same way; two with a digitising error creating a 2 cm sliver of overlap flip from touches to overlaps. That flip is a useful topology check — see topological coverage explained. Land Registry title extents overlap legitimately in places (a lease inside a freehold), so thousands of overlapping pairs there are a feature of the data model, not all errors.

crosses is for mixed dimensions

crossing = gpd.sjoin(roads, wards, predicate="crosses")
inside = gpd.sjoin(roads, wards, predicate="within")
print(f"{crossing.index.nunique()} road segments cross a ward boundary; {inside.index.nunique()} lie within one ward")
print("do two wards ever cross?", bool(shapely.crosses(left, right).any()))

390 road segments cross a ward boundary; 6559 lie within one ward
do two wards ever cross? False

crosses requires the intersection to have a lower dimension than at least one input. A line crossing a polygon produces a line (dimension 1) where the polygon is dimension 2 — so it crosses. Two polygons overlapping produce a polygon, so polygons never cross each other; they overlap.

Ward-boundary vertices that test the edge cases

The data used on this page

Rutland's 15 wards from the ONS December 2024 full-resolution boundaries (OGL); the vertices of the Oakham North East ward as points — real coordinates lying exactly on ward boundaries; the Food Standards Agency's register of Rutland's food businesses (OGL); HM Land Registry's INSPIRE parcels (© Crown copyright and database rights; OGL); and roads from Geofabrik's OpenStreetMap extract of 1 January 2026 (© OpenStreetMap contributors, ODbL). All in British National Grid.

import json
import time
from pathlib import Path

import geopandas as gpd
import pandas as pd
import requests
import shapely
from shapely.geometry import Point, Polygon

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)
wards = gpd.read_file("data/wards.geojson")
edge = wards.loc[wards["WD24NM"] == "Oakham North East", "geometry"].iloc[0].exterior
points = gpd.GeoDataFrame(geometry=gpd.points_from_xy(*shapely.get_coordinates(edge)[:-1].T), crs=27700)

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[["BusinessName"]], geometry=gpd.points_from_xy(x, y), crs=4326)[x.notna().values].to_crs(27700)
parcels = gpd.read_file("/vsizip/data/Rutland.zip/Land_Registry_Cadastral_Parcels.gml", columns=["INSPIREID"], DOWNLOAD_SCHEMA="NO")
roads = gpd.read_file("/vsizip/data/rutland.shp.zip/gis_osm_roads_free_1.shp", columns=["name"]).to_crs(27700)
print(len(wards), "wards |", len(points), "points on the Oakham North East boundary |", len(businesses), "businesses")

15 wards | 214 points on the Oakham North East boundary | 338 businesses

Example 1: choosing a predicate by asking the right question

QUESTIONS = {
    "which ward is this business in?":          (businesses, wards, "within"),
    "which businesses are in this ward?":       (wards, businesses, "contains"),
    "which wards does this road pass through?": (roads, wards, "intersects"),
    "which wards share a border with this one?": (wards, wards, "touches"),
    "which boundary points are on or in a ward?": (wards, points, "covers"),
}
for question, (left_frame, right_frame, predicate) in QUESTIONS.items():
    print(f"{question:45} {predicate:10} {len(gpd.sjoin(left_frame, right_frame, predicate=predicate)):>6,} rows")

which ward is this business in?               within        338 rows
which businesses are in this ward?            contains      338 rows
which wards does this road pass through?      intersects  7,222 rows
which wards share a border with this one?     touches        62 rows
which boundary points are on or in a ward?    covers        432 rows

Writing the question in words first is not a formality — it is the only way to notice that "in the ward" and "on or in the ward" are different questions with different answers.

Example 2: counting the boundary cases before they surprise you

def boundary_report(points, polys):
    """How many points sit exactly on a boundary? That is the predicate risk."""
    within = gpd.sjoin(points, polys, predicate="within")
    covered = gpd.sjoin(points, polys, predicate="covered_by")
    return {"points": len(points), "matched_within": within.index.nunique(),
            "matched_covered_by": covered.index.nunique(),
            "on_boundary": covered.index.nunique() - within.index.nunique(),
            "rows_covered_by": len(covered)}

print("businesses:     ", boundary_report(businesses.dropna(subset=["geometry"]), wards))
print("boundary points:", boundary_report(points, wards))

businesses:      {'points': 338, 'matched_within': 338, 'matched_covered_by': 338, 'on_boundary': 0, 'rows_covered_by': 338}
boundary points: {'points': 214, 'matched_within': 0, 'matched_covered_by': 214, 'on_boundary': 214, 'rows_covered_by': 432}

The businesses have no boundary cases at all; the boundary points are nothing else, and covered_by matches most of them to two wards. Whether a point on a line belongs to one ward, the other, or both is a decision for you — but now you know the decision exists.

Example 3: a duplicate-safe point-in-polygon join

A point on a shared boundary is matched by both neighbouring polygons under covered_by (the left-hand form of covers), producing two rows. That is correct behaviour and usually not what you want.

joined = gpd.sjoin(points, wards, predicate="covered_by", how="left")
dupes = joined.index.duplicated(keep=False)
print(f"{dupes.sum()} rows from points on a shared boundary")

# keep one deterministically — the alphabetically first ward, so re-runs agree
ordered = joined.sort_values("WD24NM", kind="stable")
one_each = ordered[~ordered.index.duplicated(keep="first")].sort_index()
print(len(points), "points ->", len(joined), "rows ->", len(one_each), "after tie-breaking;",
      one_each["WD24NM"].value_counts().to_dict())

432 rows from points on a shared boundary
214 points -> 432 rows -> 214 after tie-breaking; {'Oakham North East': 161, 'Exton': 31, 'Barleythorpe': 22}

Deterministic tie-breaking matters more than which side wins. A join that picks differently each run makes every downstream total unreproducible — see idempotency explained.

Predicates are defined by DE-9IM

Spatial predicates come from the OGC Simple Features specification, and their precision is deliberate: they are defined against the DE-9IM matrix so that any two implementations agree exactly. GEOS, PostGIS, JTS and Shapely all give the same answer for the same pair of shapes, which is why the same predicate names appear in SQL and in Python.

The design decision that surprises people is that contains requires interior intersection. A.contains(B) is defined as "no point of B is in the exterior of A, and at least one point of B is in the interior of A". For a point exactly on the ring, the second clause fails — the point is in A's boundary, not its interior. covers drops that second clause, which is why it is the more forgiving of the two.

This is not pedantry with a practical cost of zero. Administrative boundaries are frequently digitised from the same source as the features inside them, so vertices coincide exactly. Points snapped to a road network sit precisely on the road that forms a ward edge. In those datasets, the within and covers counts differ by a real number of rows, every time.

The second thing worth internalising: Speed is not a reason to prefer intersects. GeoPandas uses the spatial index to find bounding-box candidates for every predicate, then runs the exact test, and which exact test is cheapest depends on the data. On Rutland's parcels against its wards the stricter predicates were not slower — within was the fastest by far — so choose by the question, not by an assumption about speed:

for predicate in ("intersects", "within", "covered_by"):
    t0 = time.perf_counter()
    rows = len(gpd.sjoin(parcels, wards, predicate=predicate))
    print(f"{predicate:10} {rows:>6,} rows in {time.perf_counter() - t0:.2f} s")

intersects 21,843 rows in 0.52 s
within     20,275 rows in 0.03 s
covered_by 20,275 rows in 0.43 s

Identical geometries, precision and dwithin

• within on identical geometries is True. A shape is within itself, and contains itself. Self-joins need index != index_right filtering.
• overlaps is False for identical geometries — neither contains the other is violated, since each contains the other.
• Invalid geometries make every predicate unreliable. A self-intersecting polygon has an ambiguous interior. Run make_valid first.
• touches is very sensitive to coordinate precision. Two parcels that "should" share an edge may miss by 10⁻⁹ and report False. See coordinate precision.
• Proximity is not a topological relationship. None of the DE-9IM predicates says "within 500 m"; GeoPandas 1.0+ adds dwithin for that, mirroring PostGIS's ST_DWithin, and sjoin_nearest answers "closest".
• Empty geometries return False for everything including intersects, and never raise. Filter them before joining.

ward = wards.geometry.iloc[0]
print("identical: within", ward.within(ward), "| overlaps", ward.overlaps(ward),
      "| empty intersects", shapely.Polygon().intersects(ward))

identical: within True | overlaps False | empty intersects False

Internal links

• How to perform a spatial join in Python (GeoPandas) — where the predicate argument lives
• Spatial join returns empty results in GeoPandas — when the predicate is right and the CRS is not
• GeoPandas spatial join returns duplicate rows — the boundary-point duplication in example 3
• How to select features by location in GeoPandas — predicates outside a join
• Topological Coverage: What "Clean" Means for Adjacent Polygons
• What makes a geometry valid? — why invalid input breaks predicates
• Spatial indexes explained — the candidate filter every predicate runs first
• How to run spatial SQL queries from Python with PostGIS — the same predicates as ST_ functions

FAQ

Which predicate should I use for point-in-polygon?

within (or contains with the frames swapped) for strict containment, covers if a point sitting exactly on the boundary should count. Compare the two counts once on your data to see whether it matters.

Why is contains False for a point on the edge?

contains requires at least one point of the inner shape to be in the outer shape's interior. A point on the ring is in the boundary, not the interior. Use covers, which drops that requirement.

What is the difference between intersects and overlaps?

intersects is true for any shared point, including a boundary-only touch or full containment. overlaps requires the interiors to intersect and neither shape to contain the other.

Do two polygons ever cross?

No. crosses requires the intersection to have a lower dimension than at least one input; two overlapping polygons intersect in a polygon. Use overlaps.

Is intersects faster than the others?

Not reliably. On Rutland's parcels against its wards, within ran in a tenth of the time intersects took. All predicates use the spatial index for candidate selection first; after that, the cost depends on the geometries and on GEOS's prepared-geometry shortcuts more than on the predicate's name.

How do I find features within a distance?

No predicate does distance. Buffer one layer and use intersects, use sjoin_nearest with max_distance, or use dwithin in GeoPandas 1.0+.

Are these the same as PostGIS ST_ functions?

Yes — ST_Intersects, ST_Within, ST_Covers and the rest implement the same OGC definitions, so the answers match. Only the performance characteristics differ.