A conventional join needs a shared key. Two tables both hold a customer identifier, the identifier matches or it does not, and the result is unambiguous. Spatial data usually has no such key. A dataset of incidents and a dataset of districts share nothing but the fact that both describe positions on the same planet.
A spatial join solves that by making geometry the join condition. Rows match because one shape contains, touches, crosses, or lies near another. Almost every question that sounds like location analysis is a spatial join underneath: how many customers per district, which properties are inside the hazard zone, what is the nearest depot, how much road runs through each catchment.
The consequence is that the join condition becomes a modelling decision rather than an equality test. Two records either share a key or they do not. Two shapes relate to each other in a way that has to be chosen, and different choices give different, equally valid, mutually contradictory answers.
Joining without a key
The mechanics are the same as a relational join. For each row on the left, find the rows on the right that satisfy a condition, and emit the combination. Only the condition changes: instead of matching identifiers, it asks whether the left geometry falls inside the right one.
That substitution has three consequences worth stating plainly. The condition is not symmetric the way key equality is: point within polygon and polygon contains point describe the same relationship, but which side you join from determines the shape of the result and which side produces nulls.
The condition can also match many rows on both sides. Key joins are frequently many to one by construction. Spatial joins are many to many by default, and the cardinality depends on the data rather than the schema.
And it has edge cases that equality does not. Geometries can touch without overlapping, share a boundary without sharing an interior, or be separated by less than the precision of the coordinates. Each has to resolve to true or false, and the resolution is a convention.
The predicates, and what each one asserts
The standard predicates come from a formal model of how two geometries relate through their interiors, boundaries, and exteriors. A handful cover almost all work.
- Intersects. The two geometries share at least one point, including a boundary point. The loosest useful predicate and the one most engines optimise best.
- Contains and within. One geometry lies entirely inside the other, the same relationship viewed from opposite sides. The strict definition normally excludes cases where the inner geometry touches the boundary of the outer one, which surprises people the first time it matters.
- Covers and covered by. Like contains and within, but permitting boundary contact. Where many points sit exactly on shared edges, the difference between contains and covers changes the counts.
- Touches. The geometries share a boundary but no interior. Used for adjacency: which districts border this one, which parcels abut it.
- Crosses and overlaps. Partial interior intersection, the meaning depending on the dimensions of the two geometries.
- Within a distance. Not strictly a topological predicate. It measures, and it introduces a unit, which is where a great many silent errors enter.
Choosing between these is not stylistic. Counting points in a polygon with intersects rather than within includes boundary points in both neighbouring polygons. Counting with within rather than covers drops points sitting exactly on a boundary from both. The same input produces three different totals, and all three are defensible.
Nearest is a different animal
Nearest-neighbour joins look like the same operation with a different predicate. They are not. Every predicate above is a yes or no test on a pair of geometries. Nearest is a ranking across all candidates, and it always returns something.
That guarantee is the danger. A containment join can return nothing, which is informative: either the point is wrong, the coverage is incomplete, or the answer really is none. A nearest join never returns nothing. It returns the closest record even when that record is fifty kilometres away and has no relationship to the query point.
In dense areas, nearest joins fail worse. When candidates are packed closely, the gap between first and second nearest can be smaller than the positional error of either dataset. The join is then assigning at random, but deterministically, so the same wrong answer appears on every run and looks stable.
Consider matching customer locations to the nearest branch in a city centre where four branches sit within a few hundred metres. The customer positions carry tens of metres of error and some are centroid fallbacks. The nearest computation is precise. Its inputs are not. The resulting attribution of customers to branches is noise presented as a table.
Three mitigations, in order of usefulness. Cap the distance, so a nearest join finding nothing within a sensible radius returns null rather than an absurd match. Return the distance alongside the match and force consumers to look at it. Return the runner-up and its distance, because when the first and second are within the noise floor, the correct output is an ambiguity flag rather than a winner.
Straight-line nearest also ignores the network. The nearest facility by air can be across a motorway with no crossing, or inside a compound whose only gate is on the far side. Where the question is about access rather than proximity, distance has to be measured along the network the movement actually uses.
One to many, and the double count
The most common analytical error in spatial joins is not a wrong predicate. It is forgetting that the join changed the cardinality of the left table.
A point inside a single district joins to one row. A point inside overlapping catchments joins to several. A road crossing five districts joins to five rows. A footprint straddling a boundary joins to both sides. In each case the left table has grown, and every subsequent aggregate over it is weighted by how many right-hand features each left-hand feature happened to hit.
Sum a revenue column after a one to many join and the revenue of every multiply-matched record is counted once per match. The total is inflated, and the inflation concentrates in the records near boundaries or inside overlapping zones, so it looks like a real geographic pattern.
The defence is to decide, before the join, what the correct cardinality is and enforce it explicitly.
- If each left record must attribute to exactly one right record, define the tie-break up front: largest overlap area, smallest containing polygon, highest priority zone. Then assert one row per left key after the join and fail loudly if the assertion breaks.
- If a left record genuinely spans several right records, apportion the value rather than attributing all of it to each. A road crossing five districts contributes its length per district, not its full length five times. A building straddling a boundary contributes its area proportionally, or is assigned by centroid, but the choice is made and recorded.
- If the join is intentionally many to many, aggregate on the right side, and never sum a left-side measure without deduplicating first.
This error persists because the row count after a join is rarely checked. Comparing the count before and after, and knowing which one you expected, catches most of it in one line.
The cost model, and why spatial indexes matter
A join on equality can use a hash table. A join on geometry cannot, because there is no value to hash: the condition depends on shapes, not on any scalar.
Without an index, the engine tests every left geometry against every right geometry. The number of tests is the product of the two table sizes, and each test is an exact predicate over potentially thousands of vertices. This is why an unindexed spatial join on modest tables can run for hours while the equivalent key join finishes instantly.
Spatial indexes make it tractable through a two-phase model that explains most of the performance behaviour people observe.
The filter phase uses approximations. Each geometry is reduced to its bounding box, and the boxes are organised into a tree or a grid so that candidate pairs whose boxes cannot interact are eliminated without any geometric computation. This phase is fast and produces false positives but never false negatives. The refinement phase then runs the exact predicate on the surviving pairs only. That is where the real cost sits, and it scales with vertex count.
Two consequences follow. First, geometries whose bounding boxes approximate their shape poorly defeat the filter phase. A long diagonal river, a coastal administrative area with a huge offshore box, or a multipolygon scattered across a country all produce boxes that overlap almost everything, pushing enormous numbers of pairs into the expensive phase. Splitting such geometries into parts, or indexing the parts, is often the single biggest performance win available.
Second, vertex count is a first-class cost. Joining against boundaries digitised at survey precision when the analysis needs only neighbourhood resolution wastes most of the compute in the refinement phase. Simplifying the right-hand geometries to the tolerance the question requires can change runtime by an order of magnitude, at the price of moving some near-boundary records to the wrong side. That is a trade to make deliberately.
For very large joins, spatial partitioning applies the same principle to the whole job: divide the space into tiles, distribute both sides by tile, join within tiles, and handle straddling geometries explicitly. The straddling cases are where distributed spatial joins produce duplicates, which loops back to the counting problem above.
The point that sits exactly on the line
Coordinates are floating point numbers and boundaries are shared. Together those facts generate a category of bug that is hard to reproduce and easy to dismiss.
When two adjacent polygons were digitised from the same source they usually share vertices exactly, and a point placed on that shared edge is genuinely on the boundary of both. Whether it joins to one, both, or neither depends on the predicate and on the engine's convention, and a point on a vertex can behave differently from a point on an edge.
The situations where this stops being theoretical are common. Addresses interpolated onto a street centreline that is also the administrative boundary put many points exactly on the line. Grid-aligned data joined to grid-aligned boundaries puts entire rows of points on edges. Reprojected and rounded data can land on boundaries it was originally just inside.
The robust approach is not to hope. Pick a predicate with defined boundary behaviour and use it consistently. Where a single definitive assignment is needed, apply a deterministic tie-break so repeated runs agree. Where boundary cases are numerous enough to matter, count them and report them as their own category rather than letting them dissolve into whichever side won.
Two datasets describing the same boundary from different sources will also not agree to the last decimal, and points near that boundary fall inside one version and outside the other. No predicate fixes this, because the disagreement is in the data. The only honest handling is to treat a band around the boundary as uncertain and to know how many records fall in it.
When the point placement is approximate
The classic error is joining points to polygons without asking how the points were placed.
A containment test is exact. It answers, with total precision, whether a coordinate lies inside a shape. It says nothing about whether that coordinate is where the real-world thing is.
If the points came from geocoding, a meaningful share are estimates: interpolated along street segments, taken from parcel centroids, or fallen back to locality centroids. If they came from consumer positioning, they carry tens of metres of error and more in dense urban environments. If they came from an administrative record, they may be the registered address rather than the operating site.
Join those points to fine-grained polygons and the exactness of the operation launders the vagueness of the input. The smaller the zone, the higher the proportion of assignments determined by positional error rather than by real position. Assigning approximately-placed points to census blocks or individual parcels produces results that are arbitrary at the record level while looking authoritative in aggregate.
The rule that follows is often ignored: the granularity of the join must not exceed the accuracy of the inputs. If positions are accurate to a hundred metres, joining to polygons fifty metres across is manufacturing detail. Either join to coarser zones, or carry the uncertainty through and report the assignment as a distribution over candidate zones rather than a single winner.
There is a related trap in units. A distance join written as within 500 against data stored in degrees does not mean five hundred metres. Distance joins belong in a projected coordinate system chosen for the region, or in a geodesic function that knows it is working on a spheroid.
A procedure for a join you can defend
Before running it, answer five questions and write the answers next to the code.
- What is the predicate and why. Name it explicitly, including boundary behaviour, rather than accepting the default.
- What cardinality is expected, and what the row count should be afterwards.
- What happens to non-matches. Inner or outer, and if outer, what a null means: outside all zones, or missing coverage in the zone layer.
- What coordinate system the operation runs in, and if a distance is involved, what unit it is in.
- How the input positions were produced, and what the smallest polygon in the right-hand layer is.
Afterwards, check three things. The row count against the expectation. The count of left records that matched nothing, often the most informative number in the whole operation. And the distribution of matches per left record, because a long tail of records matching twenty polygons usually means the right-hand layer has overlaps nobody knew about.
The honest limit
A spatial join computes a relationship between two geometries exactly. It cannot compute the relationship between two real-world things, because it never sees them, only their recorded shapes.
Every join inherits every error in both inputs: positional error in the points, digitising error in the boundaries, the currency of both, and the modelling decisions that turned real objects into geometry. The operation adds no error of its own, which is precisely why it is so easy to trust. The precision of the output is set by the worst of the inputs, and the output carries no memory of which one that was.