Spatial SQL or GeoPandas? Choosing Where the Work Happens

Problem statement

The same spatial join can be written two ways, and they look equivalent. Here are HM Land Registry's Leeds parcels joined to the Environment Agency's flood zones, for the stretch of the Aire valley from Kirkstall to the city centre:

AIRE = box(425_000, 432_000, 431_000, 436_000)


def in_python():
    parcels = gpd.read_postgis("SELECT * FROM parcels", engine, geom_col="geom")
    floods = gpd.read_postgis("SELECT * FROM flood_zones", engine, geom_col="geom")
    parcels = parcels[parcels.intersects(AIRE)]
    return gpd.sjoin(parcels, floods, predicate="intersects")


def in_the_database():
    return gpd.read_postgis(text("""
        SELECT p.*, f.flood_zone
        FROM parcels p
        JOIN flood_zones f ON ST_Intersects(p.geom, f.geom)
        WHERE ST_Intersects(p.geom, ST_GeomFromText(:wkt, 27700))
    """), engine, geom_col="geom", params={"wkt": AIRE.wkt})


for name, fn in [("in Python", in_python), ("in the database", in_the_database)]:
    start = time.perf_counter()
    rows = len(fn())
    print(f"{name:16} {rows:>6,} rows  {time.perf_counter() - start:6.2f} s")

in Python         3,745 rows    5.32 s
in the database   3,745 rows    0.44 s

The same rows, and very different costs. The Python version reads all 293,729 parcels across the wire and parses every one into a Shapely object, then throws almost all of them away. On a national table the first version exhausts memory before it reaches the join.

The reverse is also true. Fitting a model, plotting a map, calling a Python library on each geometry, or doing anything the SQL standard has no word for is far easier and often faster in Python. Neither tool wins in general; the decision has a structure, and it is worth knowing it rather than defaulting.

Quick answer

The governing principle is one sentence: push the reduction into the database, and pull only what you will actually use.

# not this
gdf = gpd.read_postgis("SELECT * FROM parcels", engine, geom_col="geom")
subset = gdf[gdf.intersects(area)]

# this
gdf = gpd.read_postgis(text("""
    SELECT * FROM parcels
    WHERE ST_Intersects(geom, ST_GeomFromText(:wkt, 27700))
"""), engine, geom_col="geom", params={"wkt": area.wkt})

Where the data shrinks, what each side does well, transfer cost

1. Find where the data volume drops

Every workflow has a point where the row count or byte count collapses. Everything before that point belongs in the database; everything after belongs in Python.

print(q(f"""
    SELECT (SELECT COUNT(*) FROM parcels)                                       AS all_of_leeds,
           (SELECT COUNT(*) FROM parcels
             WHERE ST_Intersects(geom, ST_GeomFromText('{AIRE.wkt}', 27700)))    AS in_study_area,
           (SELECT COUNT(DISTINCT p.id) FROM parcels p
             JOIN flood_zones f ON ST_Intersects(p.geom, f.geom)
             WHERE ST_Intersects(p.geom, ST_GeomFromText('{AIRE.wkt}', 27700)))  AS at_risk
""").to_string(index=False))

 all_of_leeds  in_study_area  at_risk
       293729          32682     1622

That is the funnel: filter to the study area in SQL, join to the flood zones in SQL, and only then hand the parcels at risk to Python for a risk index, a plot, a model. Moving the boundary one step earlier transfers the whole study area instead of the parcels at risk — survivable. Moving it to the top transfers the whole of Leeds, which the problem statement timed.

2. Know what each side is actually good at

PostGIS is faster when:

• The operation is indexed. A bounding-box-first spatial join on an indexed table is a different algorithm from a linear scan, not a faster version of one — see PostGIS spatial indexes explained.
• The data does not fit in memory. Postgres streams from disk; GeoPandas does not.
• The result is much smaller than the input. Aggregations, filters and joins that reduce.
• Several processes need the same answer. The database caches; each Python process does not.

GeoPandas is faster when:

• The data is already in memory and small enough. No round trip, no serialisation.
• The operation is vectorised NumPy over an array of geometries.
• You are iterating on the analysis. Re-running a cell beats re-running a query and re-transferring.
• The operation has no SQL expression, so doing it in SQL means a plpython function or a loop.

Both are slow when the geometry is transferred repeatedly. Serialising hundreds of thousands of polygons to WKB, sending them over a socket and parsing them into Shapely objects costs more than either engine's computation.

3. Measure the transfer, not just the query

The query time in EXPLAIN ANALYZE excludes the part that usually dominates:

def timed(sql, label, **kw):
    t0 = time.perf_counter()
    gdf = gpd.read_postgis(text(sql), engine, geom_col="geom", **kw)
    dt = time.perf_counter() - t0
    wkb = gdf.geometry.to_wkb().map(len).sum() / 1e6
    print(f"{label:<18} {len(gdf):>9,} rows  {wkb:>7.1f} MB of WKB  {dt:>6.2f} s")
    return gdf


with engine.begin() as con:
    plan = con.execute(text("EXPLAIN ANALYZE SELECT * FROM parcels")).scalars().all()
print("EXPLAIN ANALYZE says:", plan[-1])
everything = timed("SELECT * FROM parcels", "everything")
filtered = timed("SELECT * FROM parcels WHERE ST_Intersects(geom, ST_GeomFromText(:wkt, 27700))",
                 "filtered in SQL", params={"wkt": AIRE.wkt})

EXPLAIN ANALYZE says: Execution Time: 45.683 ms
everything           293,729 rows    110.1 MB of WKB    3.98 s
filtered in SQL       32,682 rows     10.6 MB of WKB    0.40 s

The database finishes its part of SELECT * FROM parcels in a fraction of a second; reading the same query into GeoPandas takes several seconds. The difference is transfer and parsing, which is why "the query is fast" is not evidence that the approach is right.

4. Reduce the geometry as well as the rows

Rows are only half the payload. A polygon with thousands of vertices costs the same to transfer whether you need its detail or not.

-- select only the columns you use
SELECT id, ward_code, geom FROM parcels WHERE …            -- not SELECT *

-- simplify when the detail will not be used
SELECT id, ST_SimplifyPreserveTopology(geom, 5) AS geom FROM parcels WHERE …

-- or drop the geometry entirely when you only need attributes
SELECT id, ward_code, ST_Area(geom) AS area_m2 FROM parcels WHERE …

That last one is the one people forget. If the analysis needs areas rather than shapes, computing the area in SQL and transferring a float instead of a polygon removes most of the payload — the large-table guide measured 116 MB of full geometry against 11.7 MB of ids and areas for these same parcels.

5. Write results back rather than round-tripping

If the output feeds another query, keep it in the database:

with engine.begin() as con:
    con.execute(text("DROP TABLE IF EXISTS parcels_at_risk"))
    con.execute(text("""
        CREATE TABLE parcels_at_risk AS
        SELECT DISTINCT ON (p.id) p.id, p.geom, f.flood_zone
        FROM parcels p
        JOIN flood_zones f ON ST_Intersects(p.geom, f.geom)
        ORDER BY p.id, f.flood_zone DESC
    """))
    con.execute(text("CREATE INDEX ON parcels_at_risk USING GIST (geom)"))
print(q("SELECT flood_zone, COUNT(*) FROM parcels_at_risk GROUP BY 1 ORDER BY 1").to_string(index=False))

flood_zone  count
       FZ2   5319
       FZ3   7790

The data never entered Python. The index on the new table matters as much as the table — CREATE TABLE AS copies rows, not indexes, so a derived table without one turns the next query into a sequential scan. DISTINCT ON keeps one row per parcel, choosing the higher of the two zones where a parcel touches both. See how to write a GeoDataFrame to PostGIS for the reverse direction.

One Leeds analysis both ways, measured

The data used on this page

All open data for Leeds, in a local PostGIS 3.4 (the postgis/postgis:16-3.4 Docker image; connection details from PGPORT and PGPASSWORD): HM Land Registry's INSPIRE parcels, each tagged with its ONS December 2024 ward; the Environment Agency's Flood Map for Planning zones 2 and 3 over Leeds; the NaPTAN stops of West Yorkshire; and police.uk's street-level crimes for June 2026 in a box around Leeds.

import io
import os
import time

import geopandas as gpd
import numpy as np
import pandas as pd
import requests
from shapely.geometry import box
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")


def q(sql):
    return pd.read_sql(text(sql), engine)


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)

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

zones = requests.get("https://environment.data.gov.uk/spatialdata/flood-map-for-planning-flood-zones/wfs",
                     timeout=300, params={
                         "service": "WFS", "version": "2.0.0", "request": "GetFeature",
                         "typeNames": "dataset-04532375-a198-476e-985e-0579a0a11b47:Flood_Zones_2_3_Rivers_and_Sea",
                         "bbox": "419000,425000,446000,450000,urn:ogc:def:crs:EPSG::27700",
                         "outputFormat": "application/json", "count": 100000})
floods = gpd.read_file(io.BytesIO(zones.content)).set_crs(27700, allow_override=True).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")

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"], float(c["location"]["longitude"]), float(c["location"]["latitude"])) for c in r.json()]
crimes = pd.DataFrame(found, columns=["id", "lon", "lat"]).drop_duplicates("id")
incidents = gpd.GeoDataFrame(crimes[["id"]], crs=4326,
                             geometry=gpd.points_from_xy(crimes["lon"], crimes["lat"])).to_crs(27700)
incidents = incidents.rename_geometry("geom")

with engine.begin() as con:
    con.execute(text("DROP TABLE IF EXISTS parcels, flood_zones, stops, incidents, wards_tag, parcels_at_risk CASCADE"))
parcels.to_postgis("parcels", engine, index=False, chunksize=50_000)
floods.to_postgis("flood_zones", engine, index=False)
stops.to_postgis("stops", engine, index=False)
incidents.to_postgis("incidents", engine, index=False)
wards.to_postgis("wards_tag", engine, index=False)
with engine.begin() as con:
    con.execute(text("ALTER TABLE parcels ADD COLUMN ward_code text"))
    con.execute(text("""UPDATE parcels p SET ward_code = w.ward_code FROM wards_tag w
                        WHERE ST_Intersects(ST_PointOnSurface(p.geom), w.geom)"""))
    for table in ("parcels", "flood_zones", "stops", "incidents"):
        con.execute(text(f"ANALYZE {table}"))
print({"parcels": len(parcels), "flood_zones": len(floods), "stops": len(stops), "incidents": len(incidents)})

{'parcels': 293729, 'flood_zones': 4223, 'stops': 16677, 'incidents': 15942}

Example 1: the same analysis, both ways, measured

def python_side():
    t0 = time.perf_counter()
    parcels = gpd.read_postgis("SELECT id, geom FROM parcels", engine, geom_col="geom")
    floods = gpd.read_postgis("SELECT flood_zone, geom FROM flood_zones", engine, geom_col="geom")
    parcels = parcels[parcels.intersects(AIRE)]
    joined = gpd.sjoin(parcels, floods, predicate="intersects")
    return len(joined), joined["id"].nunique(), time.perf_counter() - t0


def sql_side():
    t0 = time.perf_counter()
    joined = gpd.read_postgis(text("""
        SELECT p.id, p.geom, f.flood_zone
        FROM parcels p
        JOIN flood_zones f ON ST_Intersects(p.geom, f.geom)
        WHERE ST_Intersects(p.geom, ST_GeomFromText(:wkt, 27700))
    """), engine, geom_col="geom", params={"wkt": AIRE.wkt})
    return len(joined), joined["id"].nunique(), time.perf_counter() - t0


for name, fn in [("all in Python", python_side), ("reduce in SQL", sql_side)]:
    rows, parcels_hit, dt = fn()
    print(f"{name:<16} {rows:>6,} rows, {parcels_hit:>6,} parcels  {dt:>7.2f} s")

all in Python     3,745 rows,  1,622 parcels     5.09 s
reduce in SQL     3,745 rows,  1,622 parcels     0.44 s

Identical results, and a large multiple in wall clock. Note what the slow version spends its time on: not the join, which GeoPandas does perfectly well with its own spatial index, but reading a whole city of geometries it then discards. More rows than parcels in both, for the same reason: a parcel touching several flood-zone polygons — zones 2 and 3 overlap, and a long parcel can cross more than one of each — matches once for every polygon it touches.

Example 2: the hybrid that is usually right

Reduce in SQL, then do the things SQL is bad at in Python — here, grouping the parcels of one ward by their surroundings with k-means:

from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler

# ── SQL: the joins and the aggregates — all reductions ────────────────────
features = gpd.read_postgis(text("""
    SELECT p.id, p.geom,
           ST_Area(p.geom)                                               AS area_m2,
           ST_Perimeter(p.geom)                                          AS perimeter_m,
           (SELECT COUNT(*) FROM incidents i
             WHERE ST_DWithin(p.geom, i.geom, 50))                       AS crimes_within_50m,
           (SELECT ST_Distance(p.geom, s.geom) FROM stops s
             ORDER BY p.geom <-> s.geom LIMIT 1)                         AS dist_to_stop_m
    FROM parcels p
    WHERE p.ward_code = :ward
"""), engine, geom_col="geom", params={"ward": "E05011397"})                # Headingley & Hyde Park
print(f"{len(features):,} parcels in, {features.drop(columns='geom').memory_usage(deep=True).sum() / 1e6:.1f} MB "
      f"of attributes")

# ── Python: the clustering, which SQL has no word for ────────────────────
X = features[["area_m2", "perimeter_m", "crimes_within_50m", "dist_to_stop_m"]].copy()
X["compactness"] = 4 * np.pi * X.area_m2 / X.perimeter_m ** 2
features["cluster"] = KMeans(n_clusters=4, n_init=10, random_state=0).fit_predict(StandardScaler().fit_transform(X))
print(features.groupby("cluster")[["area_m2", "crimes_within_50m", "dist_to_stop_m"]].median().round(1)
      .assign(parcels=features["cluster"].value_counts().sort_index()).to_string())
ax = features.plot(column="cluster", categorical=True, legend=True, figsize=(10, 10))

7,337 parcels in, 0.3 MB of attributes
         area_m2  crimes_within_50m  dist_to_stop_m  parcels
cluster                                                     
0          123.2                0.0            93.9     3330
1           83.1                0.0           198.7     2865
2        20287.5                2.0            56.3       30
3           74.9                4.0           105.9     1112

The subqueries are doing the heavy lifting: counting crimes within 50 m and finding the nearest stop collapse to one row per parcel inside the database. Pulling the crimes and the stops into Python to count and measure them there would transfer every point to produce two columns of numbers.

ST_DWithin rather than ST_Distance(...) < 50 is deliberate — only the former uses the spatial index — and ORDER BY p.geom <-> s.geom LIMIT 1 is the index-assisted nearest-neighbour form.

Example 3: chunked reads when the result is genuinely large

Sometimes the reduced result is still too big for memory. Stream it:

def stream_postgis(sql, engine, geom_col="geom", chunk=2_000, params=None):
    """Yield GeoDataFrames of `chunk` rows using a server-side cursor."""
    with engine.connect().execution_options(
        stream_results=True, yield_per=chunk
    ) as con:
        yield from gpd.read_postgis(text(sql), con, geom_col=geom_col,
                                    params=params, chunksize=chunk)


total_area = 0.0
rows = 0
for part in stream_postgis("SELECT id, geom FROM parcels WHERE ward_code = :w",
                           engine, params={"w": "E05011397"}):
    total_area += part.area.sum()
    rows += len(part)
    print(f"  {rows:,} rows, running area {total_area / 1e6:,.2f} km²")

  2,000 rows, running area 0.89 km²
  4,000 rows, running area 1.50 km²
  6,000 rows, running area 1.97 km²
  7,337 rows, running area 2.45 km²

stream_results=True is the part that matters. Without it, psycopg buffers the entire result set client-side before the first row is handed over, so a "streaming" loop still needs memory for everything. With it, the server holds the cursor and rows arrive in batches.

This is the honest middle ground: the aggregation could have been one SUM(ST_Area(geom)) in SQL, but when the per-chunk work is something Python-only, streaming keeps memory flat. More on this in how to read a large PostGIS table without running out of memory.

Where the data is and how much moves

The choice between spatial SQL and GeoPandas is usually framed as a performance question, and it is really a question about where the data is and how much of it moves.

A database and a DataFrame make opposite assumptions. Postgres assumes the data is larger than memory, lives on disk, is shared by many clients, and will be queried in ways it cannot predict — so it builds indexes, keeps statistics, plans each query, and streams results. GeoPandas assumes the data is in memory, is yours alone, and will be manipulated in ways no query planner could anticipate — so it stores geometries in a NumPy-backed array and hands you every Python library ever written.

Neither set of assumptions is better. The mistake is running an operation under assumptions it violates: doing a whole-table scan in Python because the data "should" fit, or trying to express a machine-learning pipeline in SQL because the data "should" stay in the database.

Transfer cost is the term people leave out. A geometry crossing the boundary is serialised to WKB, sent over a socket, and parsed into a Shapely object. For a few thousand features this is invisible. For a few million it dominates everything — often exceeding the sum of both engines' computation. This single fact explains why SELECT * followed by a Python filter is so consistently the wrong shape, and why it is not obviously wrong when you write it: the query is fast, the filter is fast, and the slowness is in a step with no line of code attached to it.

The index is the other asymmetry. PostGIS's GiST index turns a spatial join from O(n×m) into something close to O(n log m), and it exists whether or not this particular query needed it. GeoPandas can build an equivalent R-tree, but it does so per process, per run, over data it has just spent minutes loading. When the same join runs a hundred times a day, the database has already paid a cost the Python version pays every time.

Finally, there is a maintenance argument that outlives any benchmark. Logic in SQL sits next to the data and is available to every consumer — a dashboard, a colleague's script, a scheduled job. Logic in a Python script is available to that script. When a rule ("a parcel is at risk if it intersects a zone-3 polygon") belongs to the data rather than to one analysis, a view or a materialised view expresses it once. When it belongs to one analysis, a script is the honest place for it. Structuring on that basis tends to put the right work in the right place without any benchmarking at all.

geom_col, SELECT * and versions

• gpd.read_postgis needs geom_col to match the geometry column's name. It defaults to geom; a table written by to_postgis calls it geometry, and the mismatch raises ValueError: Query missing geometry column 'geom'.
• SELECT * on a table with several geometry columns returns extra geometry that is silently transferred as WKB.
• PostGIS lacks direct equivalents for some Shapely operations and most of NumPy. Check before committing to a SQL-only design.
• Round-tripping data to Python and back is often worse than either option — it pays the transfer cost twice.
• A materialised view is the middle ground for an expensive query whose inputs change slowly. Refresh it on a schedule.
• ST_DWithin(a, b, d) uses the index; ST_Distance(a, b) < d does not. The same trap exists for ST_Intersects versus ST_Distance = 0.
• EXPLAIN ANALYZE measures the query, not the transfer. Time the whole Python call.
• A derived table has no index until you create one. CREATE TABLE AS copies data, not indexes.
• GeoPandas can read a query, not just a table, so pushing work into SQL never means giving up the DataFrame.
• DuckDB with the spatial extension is a third option for file-based data that is too large for memory but does not warrant a server.

Internal links

• PostGIS explained: when a spatial database beats a folder of files — the prior decision
• PostGIS spatial indexes explained — why the database wins on joins
• How to run spatial SQL queries from Python with PostGIS — the mechanics
• How to read a large PostGIS table into Python — when the result is still too big
• My PostGIS spatial query is slow — diagnosing the SQL side
• Why GeoPandas is slow: the four real bottlenecks — diagnosing the Python side
• How to do a spatial join in PostGIS with SQL — the operation that most rewards moving
• How to perform a spatial join in Python (GeoPandas) — the same join, in memory

FAQ

Is PostGIS always faster than GeoPandas?

No. On data that already fits in memory, GeoPandas is often faster because there is no query planning, no serialisation and no network. PostGIS wins when the data is large, indexed, or shared.

What is the single most important rule?

Push the step that makes the data smaller — the filter, the join, the aggregate — into SQL, and pull only what you will use. Most of the difference comes from that one decision.

Why is SELECT * then filtering in Python so slow?

The slow part is not the filter but transferring and parsing every geometry you then discard. It is invisible in profiling because no single line of your code is doing it.

Can I mix the two?

That is usually the right answer: reduce in SQL, then model, plot and export in Python. Example 2 is the standard shape.

When should I keep results in the database?

Whenever the output feeds another query, or another person. CREATE TABLE AS never brings the data into Python — but remember to create the spatial index afterwards.

Does ST_Distance use the spatial index?

No. Use ST_DWithin(a, b, distance) for proximity tests; it is index-aware and ST_Distance(a, b) < distance is not.

What about very large data with no database?

DuckDB's spatial extension, or GeoParquet with row-group filtering, handles larger-than-memory file data without running a server.