How to Do a Spatial Join in PostGIS with SQL #

Problem statement #

You know the GeoPandas version:

joined = gpd.sjoin(parcels, wards, predicate="within", how="left")

The SQL version is not obvious, because SQL has no sjoin. There is no spatial join operator — a spatial join is an ordinary JOIN whose ON clause happens to be a spatial predicate. Once that clicks, everything transfers.

What does not transfer automatically is the behaviour around the edges, and it is where the results go quietly wrong:

SELECT p.*, w.ward_name
FROM parcels p JOIN wards w ON ST_Intersects(p.geom, w.geom);

4,318,552 rows

from 4,012,884 parcels. Three hundred thousand extra rows, because parcels on a ward boundary intersect two wards and the join duplicates them. And parcels matching no ward have vanished entirely, because an inner join drops them.

Both are correct SQL. Neither is what was wanted.

Quick answer #

-- one row per parcel, keeping parcels with no match
SELECT p.id, p.geom, w.ward_name
FROM parcels p
LEFT JOIN LATERAL (
    SELECT w.ward_name
    FROM wards w
    WHERE ST_Intersects(p.geom, w.geom)
    ORDER BY ST_Area(ST_Intersection(p.geom, w.geom)) DESC
    LIMIT 1
) w ON TRUE;

From Python:

import geopandas as gpd
from sqlalchemy import create_engine

engine = create_engine("postgresql+psycopg://user@localhost/gis")
gdf = gpd.read_postgis(SQL, engine, geom_col="geom")

Step-by-step solution #

1. Make sure both tables are indexed and in the same SRID #

CREATE INDEX IF NOT EXISTS parcels_geom_idx ON parcels USING GIST (geom);
CREATE INDEX IF NOT EXISTS wards_geom_idx   ON wards   USING GIST (geom);
ANALYZE parcels; ANALYZE wards;

SELECT 'parcels' AS t, ST_SRID(geom) AS srid, COUNT(*) FROM parcels GROUP BY 2
UNION ALL
SELECT 'wards',        ST_SRID(geom),         COUNT(*) FROM wards   GROUP BY 2;

    t    | srid  |  count
---------+-------+---------
 parcels | 27700 | 4012884
 wards   | 27700 |     215

A mixed-SRID join does not return zero rows quietly the way GeoPandas does — PostGIS raises:

ERROR:  ST_Intersects: Operation on mixed SRID geometries (MultiPolygon, 27700) != (Polygon, 4326)

That is a feature. The error names both SRIDs and both geometry types, and it appears immediately rather than after you have published a map of nothing. Fix it by transforming, in the right direction:

-- transform the small table, not the big one
SELECT p.id, w.ward_name
FROM parcels p
JOIN (SELECT ward_name, ST_Transform(geom, 27700) AS geom FROM wards) w
  ON ST_Intersects(p.geom, w.geom);

Wrapping the indexed column in ST_Transform disables its index; wrapping the small table's column costs 215 transformations. See PostGIS spatial indexes explained.

2. Choose the predicate deliberately #

ST_Intersects(a, b)   -- share any point at all, boundaries included
ST_Within(a, b)       -- a is entirely inside b
ST_Contains(a, b)     -- b is entirely inside a
ST_Covers(a, b)       -- like Contains, but boundary-touching counts
ST_Overlaps(a, b)     -- partial overlap, same dimension, neither contains the other
ST_Touches(a, b)      -- share a boundary but no interior
ST_DWithin(a, b, d)   -- within d units — the index-aware proximity test

For point-in-polygon, ST_Intersects and ST_Within differ only for points exactly on a boundary — rare with real coordinates, and decisive when your points were snapped to a grid. For polygon-in-polygon, the difference is large: ST_Intersects matches a parcel touching a ward at one corner, ST_Within matches only parcels wholly inside.

The full comparison is in spatial predicates explained.

3. Decide what "one row per input" means #

A plain join returns every matching pair. If a parcel intersects three wards, it appears three times.

-- 4,318,552 rows from 4,012,884 parcels
SELECT p.id, w.ward_name
FROM parcels p JOIN wards w ON ST_Intersects(p.geom, w.geom);

Three ways to get one row per parcel, each answering a different question:

(a) The largest overlap — the usual intent for polygons:

SELECT p.id, p.geom, w.ward_name
FROM parcels p
LEFT JOIN LATERAL (
    SELECT w.ward_name
    FROM wards w
    WHERE ST_Intersects(p.geom, w.geom)
    ORDER BY ST_Area(ST_Intersection(p.geom, w.geom)) DESC
    LIMIT 1
) w ON TRUE;

(b) The containing polygon by representative point — faster, and exact when the coverage is clean:

SELECT p.id, p.geom, w.ward_name
FROM parcels p
LEFT JOIN wards w
  ON ST_Intersects(ST_PointOnSurface(p.geom), w.geom);

(c) All matches, collapsed into an array — when losing the others is not acceptable. Note that an unmatched left row produces {NULL}, an array containing a null, rather than an empty array:

SELECT p.id, p.geom,
       array_agg(w.ward_name ORDER BY w.ward_name) AS wards
FROM parcels p
LEFT JOIN wards w ON ST_Intersects(p.geom, w.geom)
GROUP BY p.id, p.geom;

Option (b) is dramatically cheaper — one point-in-polygon test instead of an intersection area per candidate — and it is correct whenever the polygons form a proper coverage.

It is not, however, a guarantee of one row per parcel. A representative point landing exactly on a shared boundary intersects both neighbours, which on PostgreSQL 16 / PostGIS 3.4 produces two rows:

 id | ward_name
----+-----------
  1 | west
  1 | east          ← the point sits at x = 10, the shared edge

Real coordinates rarely land exactly on an edge, but snapped or gridded data does it routinely. Use ST_Within instead of ST_Intersects to exclude boundary hits, or keep the lateral form in option (a) when one row per parcel must be certain. Details in spatial join cardinality explained.

4. Keep the rows that match nothing #

JOIN silently drops them. LEFT JOIN keeps them with NULL on the right, which is the only way to distinguish "no match" from "excluded":

SELECT COUNT(*) FILTER (WHERE w.ward_name IS NULL) AS unmatched,
       COUNT(*) FILTER (WHERE w.ward_name IS NOT NULL) AS matched
FROM parcels p
LEFT JOIN wards w ON ST_Intersects(ST_PointOnSurface(p.geom), w.geom);

 unmatched | matched
-----------+---------
       412 | 4012472

412 parcels outside every ward. That number is a data-quality finding — offshore features, a ward layer with gaps, or parcels in a neighbouring authority — and an inner join would have hidden it completely.

5. Verify the row count before using the result #

WITH j AS (
    SELECT p.id, w.ward_name
    FROM parcels p
    LEFT JOIN wards w ON ST_Intersects(p.geom, w.geom)
)
SELECT
    (SELECT COUNT(*) FROM parcels)                    AS input_rows,
    COUNT(*)                                          AS output_rows,
    COUNT(DISTINCT id)                                AS distinct_inputs,
    COUNT(*) FILTER (WHERE ward_name IS NULL)         AS unmatched
FROM j;

 input_rows | output_rows | distinct_inputs | unmatched
------------+-------------+-----------------+-----------
    4012884 |     4318552 |         4012884 |       412

output_rows > input_rows proves duplication. distinct_inputs < input_rows would prove loss. Both checks are two lines and both catch errors that survive review otherwise.

Code examples #

Example 1: count points in polygons, correctly #

The most common spatial join in practice, and the version most often written wrong:

-- ✅ every ward appears, including those with no incidents
SELECT
    w.ward_code,
    w.ward_name,
    w.geom,
    COUNT(i.id)                                        AS incident_count,
    COUNT(i.id) / NULLIF(ST_Area(w.geom) / 1e6, 0)     AS per_km2
FROM wards w
LEFT JOIN incidents i
       ON ST_Intersects(w.geom, i.geom)
GROUP BY w.ward_code, w.ward_name, w.geom
ORDER BY incident_count DESC;

Three details do the work. LEFT JOIN keeps wards with zero incidents — an inner join drops them, and a map missing its zero regions is a map that lies. COUNT(i.id) counts non-null right-hand rows, so an unmatched ward gets 0; COUNT(*) would count the row itself and give 1. And NULLIF(..., 0) prevents a division-by-zero error on a degenerate polygon.

From Python:

import geopandas as gpd
from sqlalchemy import create_engine

engine = create_engine("postgresql+psycopg://user@localhost/gis")
wards = gpd.read_postgis(SQL, engine, geom_col="geom")
ax = wards.plot(column="per_km2", scheme="quantiles", k=5, legend=True)

The aggregation happens in the database, so what crosses the boundary is 215 rows rather than every incident — the reduction argued for in spatial SQL or GeoPandas.

Example 2: nearest-neighbour join with the KNN operator #

"The nearest road to each parcel" is a join no ordinary predicate expresses. PostGIS has a dedicated operator:

SELECT
    p.id,
    p.geom,
    r.road_name,
    ROUND(ST_Distance(p.geom, r.geom)::numeric, 1) AS distance_m
FROM parcels p
CROSS JOIN LATERAL (
    SELECT r.road_name, r.geom
    FROM roads r
    ORDER BY p.geom <-> r.geom      -- KNN: index-assisted distance ordering
    LIMIT 1
) r
WHERE p.ward_code = 'E05011368';

<-> is the distance operator, and a GiST index makes ORDER BY … <-> … an index scan rather than a sort over every row. That is what makes this practical: without it, one nearest-neighbour lookup means computing 400,000 distances.

CROSS JOIN LATERAL runs the subquery once per parcel with p in scope. Use LEFT JOIN LATERAL … ON TRUE instead if a parcel with no road within any distance should still appear.

Cap the search when the data is sparse — an unbounded KNN in an empty region walks a long way:

CROSS JOIN LATERAL (
    SELECT r.road_name, r.geom
    FROM roads r
    WHERE ST_DWithin(p.geom, r.geom, 2000)     -- index-assisted cut-off
    ORDER BY p.geom <-> r.geom
    LIMIT 1
) r

Full treatment of the ties and radius questions in nearest-neighbour joins explained.

Example 3: a reusable spatial join from Python #

import geopandas as gpd
from sqlalchemy import create_engine, text

engine = create_engine("postgresql+psycopg://user@localhost/gis")

PREDICATES = {"intersects": "ST_Intersects", "within": "ST_Within",
              "contains": "ST_Contains", "covers": "ST_Covers"}

def spatial_join(engine, left, right, right_cols, *,
                 predicate="intersects", how="left", one_to_one=True,
                 left_geom="geom", right_geom="geom", where=None):
    """Run a spatial join in the database and return a GeoDataFrame."""
    if predicate not in PREDICATES:
        raise ValueError(f"predicate must be one of {sorted(PREDICATES)}")
    fn = PREDICATES[predicate]
    cols = ", ".join(f"r.{c}" for c in right_cols)
    filt = f"WHERE {where}" if where else ""

    if one_to_one:
        sql = f"""
            SELECT l.*, {cols}
            FROM {left} l
            LEFT JOIN LATERAL (
                SELECT {', '.join(right_cols)}, {right_geom}
                FROM {right} r
                WHERE {fn}(l.{left_geom}, r.{right_geom})
                ORDER BY ST_Area(ST_Intersection(l.{left_geom}, r.{right_geom})) DESC
                LIMIT 1
            ) r ON TRUE
            {filt}
        """
    else:
        join = "LEFT JOIN" if how == "left" else "JOIN"
        sql = f"""
            SELECT l.*, {cols}
            FROM {left} l
            {join} {right} r ON {fn}(l.{left_geom}, r.{right_geom})
            {filt}
        """

    gdf = gpd.read_postgis(sql, engine, geom_col=left_geom)
    with engine.begin() as con:
        n_left = con.execute(text(f"SELECT COUNT(*) FROM {left} {filt}")).scalar()
    print(f"{n_left:,} left rows → {len(gdf):,} joined "
          f"({len(gdf) - n_left:+,}), {gdf[right_cols[0]].isna().sum():,} unmatched")
    return gdf

parcels = spatial_join(engine, "parcels", "wards", ["ward_code", "ward_name"],
                       predicate="intersects", one_to_one=True,
                       where="l.ward_code IS NULL")

412 left rows → 412 joined (+0), 0 unmatched

The printed reconciliation is the reason to wrap this rather than writing the SQL inline each time. +0 proves no duplication; a positive number means the one-to-one guarantee failed. Making that check automatic is what stops a join quietly inflating a total three steps downstream.

Note that left and right are interpolated into the SQL string, so this function must only ever receive table names from your own code — never from user input. Parameters are bound with params= precisely because identifiers cannot be.

Explanation #

A spatial join in SQL is not a special construct. It is a JOIN whose ON clause is a function returning a boolean, and every rule about joins applies unchanged: inner joins drop non-matching rows, left joins keep them, and a join between two tables returns one row per matching pair. All the surprises come from applying those familiar rules to a predicate that matches many-to-many far more often than an equality does.

That is the essential difference from an attribute join. ON a.ward_code = b.ward_code against a unique key matches at most once. ON ST_Intersects(a.geom, b.geom) has no such guarantee: a boundary parcel touches two wards, a road crosses twenty, a buffered point contains hundreds. Duplication is not an error condition here — it is the normal output of a correct query, and the row count is the only signal that it happened.

The index makes this feasible at all, and it does so by changing the algorithm rather than the answer. Without one, a join between 4 million parcels and 215 wards is 860 million geometry comparisons. With GiST indexes, the planner takes each ward's bounding box, asks the index for parcels whose boxes overlap, and runs the exact predicate only on those candidates. Both plans return the same rows; one takes seven minutes and the other takes three seconds. The mechanism is covered in PostGIS spatial indexes explained.

LATERAL is the construct worth learning. An ordinary subquery in a FROM clause cannot see the outer row, so "the best match for this row" is inexpressible. LATERAL lifts that restriction: the subquery runs once per outer row with the outer row's columns in scope, so ORDER BY … LIMIT 1 inside it means "the best one for this row". That is what makes one-to-one spatial joins, nearest-neighbour lookups and top-N-per-group all fall out of the same pattern. LEFT JOIN LATERAL (…) ON TRUE is the form that also keeps rows whose subquery returned nothing.

Where the work happens matters as much as how it is written. A join that produces 200 aggregated rows should produce them in the database; the alternative transfers millions of geometries so Python can count them. A join that produces 4 million rows should probably not be transferred at all — write it to a table with CREATE TABLE AS, index it, and query that.

Finally, the discipline that prevents most spatial-join bugs is arithmetic rather than spatial: compare the output row count to the input row count, every time. More rows means duplication; fewer means loss. Both are usually accidents, both are invisible in a map, and both take one query to detect.

Edge cases or notes #

• Mixed SRIDs raise an error in PostGIS, unlike GeoPandas which returns empty results. Transform the smaller table.
• ST_Transform on the indexed column disables its index. Transform the constant or the small side instead.
• COUNT(*) versus COUNT(right.col) in a LEFT JOIN: the first counts unmatched rows as 1, the second as 0. Almost always you want the second.
• ST_PointOnSurface beats ST_Centroid for one-to-one assignment — a centroid can fall outside a concave polygon.
• ST_DWithin uses the index; ST_Distance(...) < d does not.
• The <-> KNN operator needs an index to be an index scan, and works on geometry and geography alike.
• Invalid geometries cause errors or wrong results. Check ST_IsValid and repair with ST_MakeValid before joining.
• JOIN LATERAL without LEFT drops rows whose subquery returns nothing — the lateral equivalent of an inner join.
• Self-joins need a.id <> b.id, or every feature matches itself.
• array_agg over an unmatched LEFT JOIN row gives {NULL}, not {}. Use array_remove(array_agg(x), NULL) if an empty array is what you want.
• Large results belong in a table: CREATE TABLE AS, then CREATE INDEX, rather than transferring millions of rows.

Internal links #

• PostGIS spatial indexes explained — why an indexed join is a different algorithm
• How to perform a spatial join in Python (GeoPandas) — the same operation in memory
• Spatial join cardinality explained — one-to-one, one-to-many and what to do
• Spatial predicates explained — choosing between intersects, within and covers
• Nearest-neighbour joins explained — distance, ties and search radius
• Spatial SQL or GeoPandas? — deciding where a join should run
• GeoPandas spatial join returns duplicate rows — the same cardinality problem in Python
• How to run spatial SQL queries from Python with PostGIS — getting the result back

FAQ #

Is there an sjoin in SQL? #

No — a spatial join is an ordinary JOIN with a spatial predicate in the ON clause. Everything you know about joins applies unchanged.

Why does my join return more rows than the input table? #

Because features match more than one row on the other side — a parcel on a ward boundary intersects both wards. Use a lateral join with LIMIT 1, or aggregate.

How do I keep rows that match nothing? #

LEFT JOIN. An inner join drops them silently, which hides genuine data-quality findings such as features outside every polygon.

What does LATERAL do? #

It lets a subquery in the FROM clause reference columns from the outer row, which makes "the best match for this row" expressible. LEFT JOIN LATERAL (…) ON TRUE also keeps rows with no match.

What happens if the SRIDs do not match? #

PostGIS raises Operation on mixed SRID geometries and names both SRIDs. Transform the smaller table, and never wrap the indexed column in ST_Transform.

How do I find the nearest feature? #

ORDER BY a.geom <-> b.geom LIMIT 1 inside a lateral join. With a GiST index this is an index scan, not a full distance computation.

Should the join run in the database or in Python? #

In the database when it reduces the data — filters, aggregates, or joins over tables too large for memory. In Python when the data is already small and the next step needs a Python library.