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. Here are HM Land Registry's 293,729 INSPIRE parcels for Leeds joined to Leeds's 33 wards:

def q(sql, **params):
    """Run a query and return a DataFrame."""
    return pd.read_sql(text(sql), engine, params=params)

print(q("""
    SELECT (SELECT COUNT(*) FROM parcels) AS parcels,
           COUNT(*)                        AS joined_rows,
           COUNT(DISTINCT p.id)            AS parcels_in_result
    FROM parcels p JOIN wards w ON ST_Intersects(p.geom, w.geom)
""").to_string(index=False))

 parcels  joined_rows  parcels_in_result
  293729       295766             293711

Over two thousand more rows than parcels, because parcels on a ward boundary intersect two wards and the join repeats them. And 18 fewer distinct parcels than went in: the ones matching no ward have vanished, 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, any of them goes straight into a GeoDataFrame:

gdf = gpd.read_postgis(SQL, engine, geom_col="geom")

tablename        indexname
  parcels idx_parcels_geom
    wards   idx_wards_geom
    wards   wards_geom_idx
    wards   idx_wards_geom
      t  srid  count
parcels 27700 293729
  wards 27700     33

Indexed tables in one SRID, a deliberate predicate

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

print(q("SELECT tablename, indexname FROM pg_indexes WHERE tablename IN ('parcels', 'wards') ORDER BY 1").to_string(index=False))
with engine.begin() as con:
    con.execute(text("ANALYZE parcels"))
    con.execute(text("ANALYZE wards"))

print(q("""
    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
""").to_string(index=False))

tablename        indexname
  parcels idx_parcels_geom
    wards   idx_wards_geom
    wards   wards_geom_idx
    wards   idx_wards_geom
      t  srid  count
parcels 27700 293729
  wards 27700     33

Both tables already have a GiST index: to_postgis builds one, idx_<table>_geom, when it creates a table, so there is nothing to add but ANALYZE. Tables loaded some other way may have none — this is the check.

Now suppose the wards had arrived in EPSG:4326. The common advice is that PostGIS, unlike GeoPandas, raises on mixed SRIDs. It depends on the plan:

with engine.begin() as con:
    con.execute(text("CREATE TABLE wards_4326 AS SELECT ward_name, ST_Transform(geom, 4326) AS geom FROM wards"))

print("indexed join:", q("SELECT COUNT(*) FROM parcels p JOIN wards_4326 w ON ST_Intersects(p.geom, w.geom)").iloc[0, 0], "rows")
try:
    q("SELECT ST_Intersects(p.geom, w.geom) FROM (SELECT geom FROM parcels LIMIT 1) p, wards_4326 w")
except Exception as exc:
    message = str(exc)
    print("direct call:", message[message.index("ST_Intersects:"):].splitlines()[0])

indexed join: 0 rows
direct call: ST_Intersects: Operation on mixed SRID geometries (Polygon, 27700) != (Polygon, 4326)

The join uses the GiST index, and the index's bounding-box test (&&) runs first. Degrees and metres never overlap, so no candidate survives to reach ST_Intersects, and the query returns zero rows without an error — exactly the GeoPandas failure. Only when ST_Intersects actually evaluates a mixed pair does PostGIS raise. Check the SRIDs yourself, and fix a mismatch 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_4326) 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 33 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, which step 1's problem statement showed. Two ways to get one row per parcel, each answering a different question, timed on all 293,729:

import time

OPTIONS = {
    "(a) largest overlap": """
        SELECT p.id, 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) point on surface": """
        SELECT p.id, w.ward_name FROM parcels p
        LEFT JOIN wards w ON ST_Intersects(ST_PointOnSurface(p.geom), w.geom)""",
}
for name, sql in OPTIONS.items():
    start = time.perf_counter()
    counts = q(f"SELECT COUNT(*) AS n, COUNT(DISTINCT id) AS ids, "
               f"COUNT(*) FILTER (WHERE ward_name IS NULL) AS unmatched FROM ({sql}) j")
    print(f"{name:22} {counts.n[0]:>7,} rows, {counts.ids[0]:>7,} parcels, "
          f"{counts.unmatched[0]:>4} unmatched, {time.perf_counter() - start:6.1f} s")

(a) largest overlap    293,729 rows, 293,729 parcels,   18 unmatched,   83.2 s
(b) point on surface   293,729 rows, 293,729 parcels,  402 unmatched,    4.0 s

Both return exactly one row per parcel. Option (b) is about twenty times faster — one point-in-polygon test per parcel instead of an intersection area per candidate ward — and it is correct whenever the polygons form a proper coverage. The two disagree on how many are unmatched, 18 against 402: a parcel straddling the Leeds boundary intersects a ward, so (a) assigns it, while its representative point may fall outside every ward, so (b) does not.

Nor is (b) a guarantee of one row per parcel. A representative point landing exactly on a shared boundary intersects both neighbours:

print(q("""
    WITH wards AS (SELECT * FROM (VALUES
            ('west', 'POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))'::geometry),
            ('east', 'POLYGON((10 0, 20 0, 20 10, 10 10, 10 0))'::geometry)) AS t(ward_name, geom)),
         parcels AS (SELECT 1 AS id, 'POLYGON((9 4, 11 4, 11 6, 9 6, 9 4))'::geometry AS geom)
    SELECT p.id, w.ward_name, ST_AsText(ST_PointOnSurface(p.geom)) AS point
    FROM parcels p JOIN wards w ON ST_Intersects(ST_PointOnSurface(p.geom), w.geom)
""").to_string(index=False))

 id ward_name       point
  1      west POINT(10 5)
  1      east POINT(10 5)

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 when one row per parcel must be certain. Details in spatial join cardinality explained.

A third form keeps every match, collapsed into an array — note that an unmatched left row produces {NULL}, an array containing a null, rather than an empty array:

SELECT p.id, 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;

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":

print(q("""
    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)
""").to_string(index=False))

 unmatched  matched
       402   293327

Those parcels lie outside every Leeds ward. Land Registry's "Leeds" file is organised by the local authority's title records, not clipped to its boundary, so it runs over the edge into the neighbouring districts. That number is a data-quality finding, and an inner join would have hidden it completely.

5. Verify the row count before using the result

print(q("""
    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
""").to_string(index=False))

 input_rows  output_rows  distinct_inputs  unmatched
     293729       295784           293729         18

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.

Counting points in polygons and nearest neighbours

The data used on this page

All open data for Leeds, loaded into a local PostGIS 3.4 (the postgis/postgis:16-3.4 Docker image; the connection details come from PGPORT and PGPASSWORD): HM Land Registry's INSPIRE Index Polygons for Leeds, the 33 Leeds wards from the ONS December 2024 ward boundaries, every street-level crime police.uk published for June 2026 in a box around Leeds, and the NaPTAN stops of West Yorkshire.

import io
import os
import time

import geopandas as gpd
import numpy as np
import pandas as pd
import requests
from sqlalchemy import create_engine, text

# the server: PostGIS 3.4 in Docker on localhost:55432
if os.system("docker start guide-postgis-55432 >/dev/null 2>&1") != 0:
    os.system("docker run -d --name guide-postgis-55432 -e POSTGRES_USER=user -e POSTGRES_PASSWORD=spatialworkflow "
              "-e POSTGRES_DB=gis -p 127.0.0.1:55432:5432 postgis/postgis:16-3.4 >/dev/null")
for _ in range(60):                                  # the image restarts once while it installs PostGIS
    if os.system("docker exec guide-postgis-55432 pg_isready -h 127.0.0.1 -U user -d gis >/dev/null 2>&1") == 0:
        break
    time.sleep(2)
os.environ.setdefault("PGPORT", "55432")
os.environ.setdefault("PGPASSWORD", "spatialworkflow")
engine = create_engine("postgresql+psycopg://user@localhost/gis")
with engine.begin() as con:
    con.execute(text("DROP TABLE IF EXISTS wards_4326"))

ONS = "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services"
reply = requests.get(f"{ONS}/Wards_December_2024_Boundaries_UK_BFC/FeatureServer/0/query", timeout=300,
                     params={"where": "LAD24NM='Leeds'", "outFields": "WD24CD,WD24NM", "outSR": 27700, "f": "geojson"})
wards = gpd.read_file(io.BytesIO(reply.content)).set_crs(27700, allow_override=True)
wards = wards.rename(columns={"WD24CD": "ward_code", "WD24NM": "ward_name"}).rename_geometry("geom")

reply = requests.Session().get("https://use-land-property-data.service.gov.uk/datasets/inspire/download/Leeds.zip",
                               timeout=600)
reply.raise_for_status()
with open("leeds_inspire.zip", "wb") as f:
    f.write(reply.content)
parcels = gpd.read_file("zip://leeds_inspire.zip!Land_Registry_Cadastral_Parcels.gml", columns=["INSPIREID"])
parcels = parcels.rename(columns={"INSPIREID": "id"}).rename_geometry("geom").set_crs(27700, allow_override=True)

west, south, east, north = wards.to_crs(4326).total_bounds
found = []
for x0, x1 in zip(np.linspace(west, east, 5)[:-1], np.linspace(west, east, 5)[1:]):
    for y0, y1 in zip(np.linspace(south, north, 5)[:-1], np.linspace(south, north, 5)[1:]):
        for attempt in range(5):                        # the API answers bursts with 429
            r = requests.post("https://data.police.uk/api/crimes-street/all-crime", timeout=120,
                              data={"poly": f"{y0},{x0}:{y0},{x1}:{y1},{x1}:{y1},{x0}", "date": "2026-06"})
            if r.status_code != 429:
                break
            time.sleep(2 ** attempt)
        r.raise_for_status()
        found += [(c["id"], c["category"], float(c["location"]["longitude"]), float(c["location"]["latitude"]))
                  for c in r.json()]
crimes = pd.DataFrame(found, columns=["id", "category", "lon", "lat"]).drop_duplicates("id")
incidents = gpd.GeoDataFrame(crimes[["id", "category"]], crs=4326,
                             geometry=gpd.points_from_xy(crimes["lon"], crimes["lat"])).to_crs(27700)
incidents = incidents.rename_geometry("geom")

# NaPTAN area 450 as saved on 9 October 2026 (the register changes weekly); live:
# naptan.api.dft.gov.uk/v1/access-nodes?dataFormat=csv&atcoAreaCodes=450
stops = pd.read_csv(io.StringIO(requests.get("https://data.spatialworkflow.io/naptan-area-450-2026-10-09.csv", timeout=120).text), low_memory=False)
stops = stops.dropna(subset=["Easting", "Northing"])
stops = gpd.GeoDataFrame({"stop_name": stops["CommonName"]}, crs=27700,
                         geometry=gpd.points_from_xy(stops["Easting"], stops["Northing"])).rename_geometry("geom")

for name, gdf in [("wards", wards), ("parcels", parcels), ("incidents", incidents), ("stops", stops)]:
    gdf.to_postgis(name, engine, if_exists="replace", index=False, chunksize=50_000)
print({name: len(gdf) for name, gdf in [("wards", wards), ("parcels", parcels),
                                         ("incidents", incidents), ("stops", stops)]})

{'wards': 33, 'parcels': 293729, 'incidents': 15942, 'stops': 16677}

Example 1: count points in polygons, correctly

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

SQL = """
    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
"""
wards_counted = gpd.read_postgis(SQL, engine, geom_col="geom")
print(wards_counted[["ward_name", "incident_count", "per_km2"]].head(3).round(1).to_string(index=False))
print(f"... {len(wards_counted)} wards, {wards_counted['incident_count'].sum():,} incidents inside them "
      f"of {q('SELECT COUNT(*) FROM incidents').iloc[0, 0]:,} fetched")
ax = wards_counted.plot(column="per_km2", scheme="quantiles", k=5, legend=True)

                ward_name  incident_count  per_km2
Little London & Woodhouse            1124    213.5
      Hunslet & Riverside             589     66.5
       Gipton & Harehills             572    143.2
... 33 wards, 9,224 incidents inside them of 15,942 fetched

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. The gap between incidents inside the wards and incidents fetched is the box around Leeds catching crimes in neighbouring districts.

The aggregation happens in the database, so what crosses the boundary is 33 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 bus stop to each parcel" is a join no ordinary predicate expresses. PostGIS has a dedicated operator:

with engine.begin() as con:                      # to_postgis indexed stops when it created the table
    con.execute(text("ANALYZE stops"))

nearest = q("""
    SELECT p.id, s.stop_name, ROUND(ST_Distance(p.geom, s.geom)::numeric, 1) AS distance_m
    FROM (SELECT p.* FROM parcels p JOIN wards w ON ST_Within(p.geom, w.geom)
          WHERE w.ward_name = 'Headingley & Hyde Park') p
    CROSS JOIN LATERAL (
        SELECT s.stop_name, s.geom
        FROM stops s
        ORDER BY p.geom <-> s.geom      -- KNN: index-assisted distance ordering
        LIMIT 1
    ) s
""")
print(f"{len(nearest):,} parcels; distance to nearest stop: median {nearest.distance_m.median():.0f} m, "
      f"max {nearest.distance_m.max():.0f} m")
print(nearest.stop_name.value_counts().head(3).to_string())

7,304 parcels; distance to nearest stop: median 132 m, max 462 m
stop_name
Royal Park Road            601
Headingley Stadium         594
Hyde Park Picture House    572

<-> 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, each lookup computes the distance to all 16,679 stops.

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 stop 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 s.stop_name, s.geom
    FROM stops s
    WHERE ST_DWithin(p.geom, s.geom, 2000)     -- index-assisted cut-off
    ORDER BY p.geom <-> s.geom
    LIMIT 1
) s

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

Example 3: a reusable spatial join from Python

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} l {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


shops = "(SELECT * FROM incidents WHERE category = 'shoplifting')"
for one_to_one in (False, True):
    spatial_join(engine, shops, "wards", ["ward_code", "ward_name"], one_to_one=one_to_one)

1,290 left rows → 1,290 joined (+0), 387 unmatched
1,290 left rows → 1,290 joined (+0), 387 unmatched

The printed reconciliation is the reason to wrap this rather than writing the SQL inline each time. For points the two modes agree unless a point sits exactly on a ward boundary; for polygons, the one-to-many mode is where + numbers appear. 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.

What the index and LATERAL change

The index makes this feasible at all, and it does so by changing the algorithm rather than the answer. Without one, joining 293,729 parcels to 33 wards means testing every pair; 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:

def timed_count(index_scans):
    with engine.begin() as con:
        con.execute(text("SET LOCAL max_parallel_workers_per_gather = 0"))
        if not index_scans:
            con.execute(text("SET LOCAL enable_indexscan = off"))
            con.execute(text("SET LOCAL enable_bitmapscan = off"))
        start = time.perf_counter()
        n = con.execute(text("SELECT COUNT(*) FROM parcels p JOIN wards w ON ST_Intersects(p.geom, w.geom)")).scalar()
        return n, time.perf_counter() - start

for index_scans in (True, False):
    n, seconds = timed_count(index_scans)
    print(f"{'with' if index_scans else 'without'} the index: {n:,} rows in {seconds:6.1f} s")

with the index: 295,766 rows in    0.8 s
without the index: 295,766 rows in  196.2 s

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.

Mixed SRIDs, transformed columns and COUNT

• Mixed SRIDs do not reliably raise. An indexed join's bounding-box test filters every candidate out first and returns zero rows; only a direct evaluation raises Operation on mixed SRID geometries. Check ST_SRID on both sides.
• 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 hundreds of thousands 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?

It depends on the plan. An indexed join returns zero rows with no error, because the bounding-box test runs first; a direct evaluation raises Operation on mixed SRID geometries. Check the 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.