Spatial SQL or GeoPandas? Choosing Where the Work Happens #

Problem statement #

The same spatial join can be written two ways, and they look equivalent:

# in Python
parcels = gpd.read_postgis("SELECT * FROM parcels", con, geom_col="geom")
floods  = gpd.read_postgis("SELECT * FROM flood_zones", con, geom_col="geom")
result  = gpd.sjoin(parcels, floods, predicate="intersects")

# in the database
result = gpd.read_postgis("""
    SELECT p.*, f.zone_class
    FROM parcels p
    JOIN flood_zones f ON ST_Intersects(p.geom, f.geom)
""", con, geom_col="geom")

On 5,000 parcels both finish in under a second and the choice does not matter. On 4 million parcels the first one transfers 6 GB over the network, exhausts memory, and takes twenty minutes if it survives at all. The second returns 180,000 rows in eleven seconds.

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", con, geom_col="geom")   # 6 GB
subset = gdf[gdf.intersects(area)]                                       # 40 MB kept

# this
gdf = gpd.read_postgis("""
    SELECT * FROM parcels
    WHERE ST_Intersects(geom, ST_GeomFromText(%(wkt)s, 27700))
""", con, geom_col="geom", params={"wkt": area.wkt})                     # 40 MB transferred

Step-by-step solution #

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.

4,000,000 parcels
      ↓  filter to the study area          ← SQL
  180,000 parcels
      ↓  join to flood zones               ← SQL
   12,400 parcels at risk
      ↓  compute a risk index              ← Python
      ↓  plot, export, model               ← Python

Moving the boundary one step earlier transfers 180,000 rows instead of 12,400 — survivable. Moving it to the top transfers 4 million and 6 GB of geometry, which is not.

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 4 million 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:

import time, geopandas as gpd
from sqlalchemy import create_engine, text

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

def timed(sql, label, **kw):
    t0 = time.perf_counter()
    gdf = gpd.read_postgis(sql, engine, geom_col="geom", **kw)
    dt = time.perf_counter() - t0
    mb = gdf.memory_usage(deep=True).sum() / 1e6
    print(f"{label:<22} {len(gdf):>9,} rows  {mb:>7.1f} MB  {dt:>6.2f} s")
    return gdf

timed("SELECT * FROM parcels", "everything")
timed("""SELECT * FROM parcels
         WHERE ST_Intersects(geom, ST_MakeEnvelope(380000,395000,400000,410000,27700))""",
      "filtered in SQL")

everything             4,012,884   6183.4 MB   412.70 s
filtered in SQL          178,204    271.9 MB    11.30 s

Thirty-seven times faster, and the difference is almost entirely transfer and parsing. This 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 40,000 vertices costs the same to transfer whether you need its detail or not.

# select only the columns you use
"SELECT id, class, 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, class, ST_Area(geom) AS area_m2 FROM parcels WHERE …"

That last line 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 can cut the payload by three orders of magnitude.

5. Write results back rather than round-tripping #

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

# creates a table from a query without the data ever entering Python
with engine.begin() as con:
    con.execute(text("""
        CREATE TABLE parcels_at_risk AS
        SELECT p.id, p.geom, f.zone_class
        FROM parcels p
        JOIN flood_zones f ON ST_Intersects(p.geom, f.geom)
    """))
    con.execute(text("CREATE INDEX ON parcels_at_risk USING GIST (geom)"))

The index on the new table matters as much as the table — a derived table without one turns the next query into a sequential scan. See how to write a GeoDataFrame to PostGIS for the reverse direction.

Code examples #

Example 1: the same analysis, both ways, measured #

import time, geopandas as gpd
from shapely.geometry import box
from sqlalchemy import create_engine

engine = create_engine("postgresql+psycopg://user@localhost/gis")
area = box(380_000, 395_000, 400_000, 410_000)

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 zone_class, geom FROM flood_zones",
                               engine, geom_col="geom")
    parcels = parcels[parcels.intersects(area)]
    joined = gpd.sjoin(parcels, floods, predicate="intersects")
    return len(joined), time.perf_counter() - t0

def sql_side():
    t0 = time.perf_counter()
    joined = gpd.read_postgis("""
        SELECT p.id, p.geom, f.zone_class
        FROM parcels p
        JOIN flood_zones f ON ST_Intersects(p.geom, f.geom)
        WHERE ST_Intersects(p.geom, ST_GeomFromText(%(wkt)s, 27700))
    """, engine, geom_col="geom", params={"wkt": area.wkt})
    return len(joined), time.perf_counter() - t0

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

all in Python      12,417 rows   487.10 s
reduce in SQL      12,417 rows    11.94 s

Identical results, forty times the 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 4 million geometries it then discards.

Example 2: the hybrid that is usually right #

Reduce in SQL, then do the things SQL is bad at in Python:

import geopandas as gpd, numpy as np
from sklearn.ensemble import RandomForestClassifier
from sqlalchemy import create_engine

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

# ── SQL: the join, the filter, and the aggregate — all reductions ──────────
features = gpd.read_postgis("""
    SELECT
        p.id,
        p.geom,
        p.building_type,
        ST_Area(p.geom)                              AS area_m2,
        ST_Perimeter(p.geom)                         AS perimeter_m,
        COUNT(t.id)                                  AS trees_within_50m,
        MIN(ST_Distance(p.geom, r.geom))             AS dist_to_road_m
    FROM parcels p
    LEFT JOIN trees t ON ST_DWithin(p.geom, t.geom, 50)
    LEFT JOIN roads r ON ST_DWithin(p.geom, r.geom, 500)
    WHERE p.ward_code = %(ward)s
    GROUP BY p.id, p.geom, p.building_type
""", engine, geom_col="geom", params={"ward": "E05011368"})

print(f"{len(features):,} parcels, {features.memory_usage(deep=True).sum()/1e6:.1f} MB")

# ── Python: the model, which SQL has no word for ──────────────────────────
X = features[["area_m2", "perimeter_m", "trees_within_50m", "dist_to_road_m"]].fillna(9_999)
X["compactness"] = 4 * np.pi * X.area_m2 / X.perimeter_m ** 2
y = features["building_type"]

model = RandomForestClassifier(n_estimators=200, random_state=0).fit(X, y)
features["predicted"] = model.predict(X)

ax = features.plot(column="predicted", legend=True, figsize=(10, 10))

2,847 parcels, 18.3 MB

The COUNT and MIN are doing the heavy lifting: two joins that could produce millions of intermediate rows collapse to one row per parcel inside the database. Pulling the trees and roads into Python to count them there would transfer hundreds of megabytes to produce two columns of integers.

ST_DWithin rather than ST_Distance(...) < 50 is deliberate — only the former uses the spatial index.

Example 3: chunked reads when the result is genuinely large #

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

import geopandas as gpd
from sqlalchemy import create_engine, text

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

def stream_postgis(sql, engine, geom_col="geom", chunk=50_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:
        for df in gpd.read_postgis(text(sql), con, geom_col=geom_col,
                                   params=params, chunksize=chunk):
            yield df

total_area = 0.0
rows = 0
for part in stream_postgis("SELECT id, geom FROM parcels WHERE ward_code = :w",
                           engine, params={"w": "E05011368"}):
    total_area += part.area.sum()
    rows += len(part)
    print(f"  {rows:,} rows, running area {total_area/1e6:,.1f} 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.

Explanation #

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.

Edge cases or notes #

• gpd.read_postgis needs geom_col to match the geometry column's name, or the result is a plain DataFrame with WKB strings.
• 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.