Spatial Join Points to Polygons in PyQGIS

"Which ward is each incident in?" is a spatial join, and it is probably the most-used spatial operation there is. The algorithm is straightforward; the decisions around it are not. A point on a boundary can match two polygons. A point in a gap matches none. And the difference between joining attributes onto points and summarising points onto polygons is the difference between two completely different output layers — both of which people describe as "a spatial join".

This recipe belongs to Vector Data Manipulation in PyQGIS. It covers the join algorithm and its predicates, the two join types, summarising the reverse direction, and what to do about unmatched and doubly-matched points.

The two cases that need a decisionPoints are scattered over three adjacent polygons. Most fall clearly inside one polygon and take its attributes. One point sits exactly on a shared boundary and satisfies the intersects predicate for both neighbours. One point lies outside every polygon and matches nothing, so the join type decides whether it survives.Most points are easy — plan for the two that are notW1W2W3on the boundary — matches W1 and W2two output rows unless you take the firstoutside every polygon — no matchdropped, or kept with null attributesCount the input and output features — if they differ, one of these two cases is why

Prerequisites

  • QGIS 3.34 LTR (bundled Python 3.12) or newer.
  • A point layer and a polygon layer. They may be in different CRSs — the algorithm reprojects — but matching them beforehand avoids surprises.
  • Valid polygon geometry; an invalid polygon can fail a predicate that visually should succeed.

Join polygon attributes onto points

import processing

processing.run("native:joinattributesbylocation", {
    "INPUT": "/data/incidents.gpkg|layername=incidents",
    "JOIN": "/data/wards.gpkg|layername=wards",
    "PREDICATE": [0],                     # intersects
    "JOIN_FIELDS": ["ward_code", "ward_name"],
    "METHOD": 1,                          # first matching feature only
    "DISCARD_NONMATCHING": False,
    "PREFIX": "",
    "OUTPUT": "/data/output/incidents_ward.gpkg",
})

Breakdown: INPUT is the layer that keeps its geometry and gains columns — points, here. PREDICATE is a list of integer codes; 0 is intersects, and for point-in-polygon work it is the correct and most permissive choice. JOIN_FIELDS limits the columns copied across; leaving it empty copies everything, which is how a point layer ends up with forty inherited attributes nobody wanted. METHOD: 1 takes only the first matching polygon, which is what makes a boundary point produce one row instead of two. DISCARD_NONMATCHING: False keeps unmatched points with null values — usually right, because a vanished point is much harder to notice than a null.

Choose the predicate deliberately

The predicate codes are positional and easy to get wrong:

CodePredicatePoint-in-polygon behaviour
0intersectsMatches inside and exactly on the boundary. The usual choice.
1containsJoin feature contains the input; equivalent for points strictly inside.
2disjointMatches points outside — useful for finding the strays.
3equalsOnly identical geometry; not useful for point-in-polygon.
4touchesBoundary only — matches a point exactly on the edge and nothing inside.
5overlapsRequires partial overlap of same-dimension geometry; never matches a point in a polygon.
6withinInput is within the join feature — the strict "inside" test.
7crossesRequires the geometries to cross; not applicable to points in polygons.

Breakdown: intersects and within differ only on the boundary: a point exactly on the edge intersects but is not within. That distinction accounts for a small number of features in most datasets and for an enormous share of the confusion. Pick intersects unless you have a reason, and use disjoint as a quick audit to list the points that fall outside your coverage.

Summarise in the other direction

The reverse question — "how many incidents in each ward?" — keeps the polygons and adds counts. That is a different algorithm.

processing.run("native:joinbylocationsummary", {
    "INPUT": "/data/wards.gpkg|layername=wards",
    "JOIN": "/data/incidents.gpkg|layername=incidents",
    "PREDICATE": [1],                      # contains
    "JOIN_FIELDS": ["severity"],
    "SUMMARIES": [0, 6],                   # count, mean
    "DISCARD_NONMATCHING": False,
    "OUTPUT": "/data/output/wards_incidents.gpkg",
})

Breakdown: Here the polygons are the INPUT, so the output has one row per ward with its geometry intact. SUMMARIES selects which statistics to compute for each field in JOIN_FIELDS — 0 is count and 6 is mean, and the output fields are named severity_count and severity_mean. Wards with no incidents keep a count of zero rather than disappearing, which is what DISCARD_NONMATCHING: False guarantees and what any subsequent choropleth map needs in order to render every ward.

Which layer keeps its geometry decides the algorithmJoining attributes by location keeps the input layer's geometry, so points stay points and gain ward columns. Joining by location with a summary also keeps the input's geometry, but with polygons as the input the result is one row per polygon carrying counts and means of the points inside it.The INPUT layer is the one that survivesjoinattributesbylocationINPUT: pointsJOIN: wardsone row per incident, geometry unchangednew columns: ward_code, ward_namejoinbylocationsummaryINPUT: wardsJOIN: pointsone row per ward, geometry unchangednew columns: severity_count, severity_mean

Audit the matches

A spatial join that silently loses or duplicates rows is a data-quality problem wearing an algorithm's clothes. Check the counts.

from qgis.core import QgsVectorLayer

points = QgsVectorLayer("/data/incidents.gpkg|layername=incidents", "in", "ogr")
joined = QgsVectorLayer("/data/output/incidents_ward.gpkg|layername=incidents_ward", "out", "ogr")

print("input", points.featureCount(), "output", joined.featureCount())

unmatched = [f["fid"] for f in joined.getFeatures() if f["ward_code"] is None]
print(f"{len(unmatched)} points matched no ward")

Breakdown: Equal counts with some nulls mean unmatched points were kept, which is the expected outcome of the parameters above. A larger output count means METHOD: 0 allowed multiple matches — every boundary point produced one row per polygon it touched. A smaller count means DISCARD_NONMATCHING was True. Printing this pair after every join takes two lines and removes an entire category of quiet errors, in the spirit of the plausibility checks in Handle Errors and Logging in Unattended Scripts.

Do it in the database when the data lives there

If both layers are already in PostGIS, the join belongs there — one indexed query rather than two table reads and an in-memory match.

SELECT i.id, i.geom, w.ward_code
FROM public.incidents i
LEFT JOIN public.wards w ON ST_Intersects(w.geom, i.geom)

Breakdown: LEFT JOIN is the SQL equivalent of keeping non-matching points, and PostGIS uses the GiST index on wards.geom to find candidates rather than testing every pair. Wrapped as a query layer it becomes an ordinary QGIS layer that is always current — see Load a PostGIS Query Layer in PyQGIS.

Rescue the points that fall just outside

Real point data misses. A GPS reading lands two metres into the sea, an address is geocoded to the centre of the road rather than the parcel, a boundary was digitised at a coarser scale than the points. Those features are not errors to discard — they belong somewhere obvious, and a nearest-feature join says where.

import processing

processing.run("native:joinbynearest", {
    "INPUT": "/data/output/unmatched.gpkg|layername=unmatched",
    "INPUT_2": "/data/wards.gpkg|layername=wards",
    "FIELDS_TO_COPY": ["ward_code", "ward_name"],
    "DISCARD_NONMATCHING": False,
    "PREFIX": "",
    "NEIGHBORS": 1,
    "MAX_DISTANCE": 25,
    "OUTPUT": "/data/output/unmatched_nearest.gpkg",
})

Breakdown: MAX_DISTANCE is the honesty control: it caps how far a point may be from a polygon before the algorithm refuses to guess. Twenty-five metres might be reasonable for GPS drift against ward boundaries and absurd for parcel assignment, so the number is a decision about the data rather than a default to accept. NEIGHBORS: 1 returns only the closest match. The output gains a distance field, and keeping it is what separates a defensible assignment from an invented one — a report can then state how many points were matched exactly and how many were snapped, with the largest distance involved.

Run this as a second pass over the unmatched points rather than instead of the containment join. Containment is exact and should win wherever it applies; nearest is a fallback whose results deserve to be labelled as such, so downstream analysis can exclude them if the question demands precision.

Exact first, nearest second, distance recordedThe first pass assigns points by containment. The points that matched nothing go through a second pass that finds the nearest polygon within a capped distance and records how far away it was. Points beyond the cap remain unassigned rather than being guessed at.A guess with a recorded distance is not the same as a match12 480 pointsto assigncontained — 12 431exact, distance = 0unmatched — 49second passnearest within 25 m — 44distance kept · 5 left unassignedKeep the distance column so a later analysis can exclude the snapped points if it needs to

QGIS version compatibility

The examples target QGIS 3.34 LTR (Python 3.12).

QGIS versionPythonNotes
3.22 LTR3.9Both algorithms present; predicate codes identical.
3.28 LTR3.9Behaviour matches this page.
3.34 LTR3.12Baseline for this page.
3.40 / 3.443.12Unchanged; the summary algorithm gained additional statistics.

Troubleshooting

  • Every point matched nothing. The CRSs differ in a way the algorithm could not resolve, or the layers do not actually overlap. Compare extents first.
  • The output has more rows than the input. METHOD: 0 allows multiple matches. Use METHOD: 1 if one row per point is required.
  • Some points vanished. DISCARD_NONMATCHING is True.
  • Points on boundaries are assigned inconsistently. With intersects and one-match-only, the winner is whichever polygon is visited first. Use within to exclude boundaries, or snap the points off the edges deliberately.
  • The join is very slow. The polygon layer has no spatial index, or the polygons are extremely detailed. Build an index, and consider simplifying — see Simplify Geometry in PyQGIS.
  • Joined field names have a prefix you did not ask for. A name collided with an existing field. Set PREFIX explicitly so the collision is visible rather than automatic.

Conclusion

Point-in-polygon work is native:joinattributesbylocation with the intersects predicate, an explicit field list and METHOD: 1 for one row per point; the reverse question is native:joinbylocationsummary with the polygons as input. Keep non-matching features so gaps are visible, compare the input and output counts every time, and push the join into PostGIS when both layers already live there.

Frequently Asked Questions

Which layer should be the INPUT? Whichever one you want the output to look like. The input keeps its geometry and its rows; the join layer only contributes attributes.

How do I find points that fall outside every polygon? Run the join with DISCARD_NONMATCHING: False and filter for nulls, or run a separate join with the disjoint predicate.

Can I join by nearest polygon instead of containment? Yes — native:joinbynearest attaches the closest feature and records the distance, which is the right tool for points that should have been inside but are not.

Does the join work across different CRSs? The algorithm reprojects the join layer to the input's CRS. It is still worth aligning them beforehand so the transformation is a deliberate, documented step.

Why do two points at the same location get different wards? They do not, unless the coordinates differ at a precision you cannot see. Print the coordinates to full precision before concluding the join is inconsistent.