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.

gpd.sjoin(points, wards, predicate="intersects")   # 4,118 rows
gpd.sjoin(points, wards, predicate="within")       # 4,096 rows
gpd.sjoin(points, wards, predicate="contains")     # 0 rows

Three predicates, three answers, 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. The three that matter most:

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)}  "
          f"contains={square.contains(p)}  covers={square.covers(p)}")

# inside   intersects=True   contains=True   covers=True
# on_edge  intersects=True   contains=False  covers=True
# outside  intersects=False  contains=False  covers=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 #

Step-by-step solution #

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:

square.relate(on_edge)      # '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, the fastest, and the right default for "find things near/on/in each other".

gpd.sjoin(points, wards, predicate="intersects")

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

within and contains are the same question, reversed #

a.within(b)      # is A inside B?
b.contains(a)    # is B holding A?   — identical truth value

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

# each point tagged with the ward it falls in
gpd.sjoin(points, wards, predicate="within")

# each ward tagged with the points it holds
gpd.sjoin(wards, points, predicate="contains")

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

The catch: contains requires at least one point of A in B's interior. A polygon does not contain its own boundary, so a point sitting exactly on the edge is not "within" the ward. With floating-point coordinates and boundaries snapped from the same source, that happens far more often than intuition suggests.

covers is the fix for boundary cases #

square.contains(on_edge)   # False — no interior intersection
square.covers(on_edge)     # True  — "no point of A is outside B"

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:

gpd.sjoin(wards, points, predicate="covers")

touches and overlaps are mutually exclusive #

a.touches(b)     # they share boundary and NOTHING else
a.overlaps(b)    # their interiors intersect, but neither contains the other

Two adjacent parcels sharing a fence line: touches is True, overlaps is False, intersects is True. Two parcels with a digitising error creating a 2 cm sliver of overlap: touches is now False, overlaps is True. That flip is a useful topology check — see topology in GIS explained.

# find parcels that should be adjacent but actually overlap
neighbours = gpd.sjoin(parcels, parcels, predicate="overlaps")
neighbours = neighbours[neighbours.index != neighbours.index_right]
print(f"{len(neighbours)} overlapping pairs — expected 0")

crosses is for mixed dimensions #

road.crosses(ward)     # a line passing through a polygon boundary
river.crosses(county)

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.

Code examples #

Example 1: choosing a predicate by asking the right question #

QUESTIONS = {
    "which ward is this point in?":            ("points", "wards", "within"),
    "which points are in this ward?":          ("wards", "points", "contains"),
    "which wards does this road pass through?": ("roads", "wards", "intersects"),
    "which parcels share a fence with this one?": ("parcels", "parcels", "touches"),
    "which parcels wrongly overlap?":          ("parcels", "parcels", "overlaps"),
    "which points are on or in the ward?":     ("wards", "points", "covers"),
}

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."""
    joined_within = gpd.sjoin(points, polys, predicate="within")
    joined_covers = gpd.sjoin(points, polys, predicate="covers", how="left")

    on_boundary = len(joined_covers) - len(joined_within)
    return {
        "points": len(points),
        "matched_within": len(joined_within),
        "matched_covers": len(joined_covers),
        "on_boundary": on_boundary,
    }

print(boundary_report(points, wards))
# {'points': 4118, 'matched_within': 4096, 'matched_covers': 4118, 'on_boundary': 22}

Twenty-two points sit exactly on a ward line. Whether they belong 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 covers-matched by both neighbouring polygons, producing two rows. That is correct behaviour and usually not what you want.

joined = gpd.sjoin(points, wards, predicate="covers", how="left")

dupes = joined.index.duplicated(keep=False)
print(f"{dupes.sum()} rows from boundary points")

# keep one deterministically — alphabetical, so re-runs agree
joined = (
    joined.sort_values(["index_right"])
          .loc[~joined.sort_values(["index_right"]).index.duplicated(keep="first")]
)

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

Explanation #

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: intersects is cheap, the others are not necessarily. GeoPandas uses the spatial index to find bounding-box candidates for every predicate, then runs the exact test. The exact test for intersects can short-circuit on the first shared point; within and contains must examine the full relationship. On large joins that difference shows up, which is another reason to use intersects unless a stricter predicate answers a question you actually asked.

Edge cases or notes #

• 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.
• Predicates say nothing about distance. For "within 500 m", buffer first or use sjoin_nearest — no predicate expresses proximity.
• CRS must match. A predicate between frames in different CRS is meaningless and GeoPandas will warn; a mismatch usually returns zero rows. See spatial join returns empty results.
• dwithin exists in GeoPandas 1.0+ for distance joins, mirroring PostGIS's ST_DWithin.
• Empty geometries return False for everything including intersects, and never raise. Filter them before joining.

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
• Topology in GIS explained: slivers, gaps and shared boundaries — why touches flips to overlaps
• 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? #

Usually slightly, because the exact test can stop at the first shared point. All predicates use the spatial index for candidate selection first, so the difference is small next to the index lookup.

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.