⚙️ How-To PostGIS & databases Joins & relationships Spatial SQL & PostGIS · step 9 of 15 Updated

How to Run Spatial SQL Queries from Python with PostGIS

Problem statement

The data is in PostGIS. The code still reads all of it.

parcels = gpd.read_postgis("SELECT * FROM parcels", engine, geom_col="geom")
wards   = gpd.read_postgis("SELECT * FROM wards", engine, geom_col="geom")

leith = wards[wards.name == "Leith"]
joined = gpd.sjoin(parcels, leith, predicate="within")
result = joined.groupby("ward")["area_m2"].sum()

Four million parcels crossed the network so that eight hundred could be kept. The spatial index the database maintains was never consulted. The machine running the script needs 6 GB of RAM to answer a question the database could have answered in 40 ms.

This is the most common way a PostGIS migration fails to pay off: the database becomes an expensive file server. The fix is not to abandon GeoPandas — it is to move the selection, the join and the aggregation into SQL and let GeoPandas receive the answer.

Quick answer

Push the filter, the join and the aggregation down; keep the modelling in Python:

import geopandas as gpd
from sqlalchemy import create_engine

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

sql = """
    SELECT p.parcel_id,
           w.name AS ward,
           ST_Area(p.geom) AS area_m2,
           p.geom
    FROM parcels p
    JOIN wards w ON ST_Within(p.geom, w.geom)
    WHERE w.name = %(ward)s
      AND p.surveyed >= %(since)s
"""

gdf = gpd.read_postgis(
    sql, engine, geom_col="geom",
    params={"ward": "Leith", "since": "2024-01-01"},
)

The rules of thumb:

Do in SQL Do in Python
WHERE, LIMIT, bbox filters anything with scikit-learn in it
spatial and attribute joins iterative or recursive logic
GROUP BY and aggregation plotting and reporting
ST_Intersection, ST_Buffer on many rows one-off geometry work on a handful
anything that reduces the row count anything that needs a Python library

And one rule that is not a rule of thumb: never build SQL with f-strings. Pass params and let the driver handle quoting.

The two shapes of the same analysis

Two panels comparing read-everything-then-filter in Python against filter-and-join-in-SQL.
Same answer. One moves 4 million rows, the other moves 812.

Step-by-step solution

Vertical steps from writing the query with placeholders through EXPLAIN, reading into GeoPandas and verifying the CRS.
EXPLAIN before you optimise. The index is usually the whole story.

1. Parameterise, always 

# never — a ward called O'Brien's Field breaks this, and worse things break it deliberately
sql = f"SELECT * FROM parcels WHERE ward = '{ward}'"

# psycopg / SQLAlchemy with a raw connection: %(name)s placeholders
gdf = gpd.read_postgis(
    "SELECT * FROM parcels WHERE ward = %(ward)s AND surveyed > %(since)s",
    engine, geom_col="geom",
    params={"ward": ward, "since": since},
)

Placeholders are for values, not identifiers. A table or column name cannot be parameterised — if it must be dynamic, validate it against a known list rather than interpolating user input.

2. Use the spatial predicates that use the index

ST_Intersects, ST_Within, ST_Contains and ST_DWithin are index-accelerated: PostGIS filters by bounding box using the GIST index, then runs the exact test on the survivors. ST_Distance in a WHERE clause is not — it computes a distance for every row.

-- scans every row, computes 4 million distances
WHERE ST_Distance(p.geom, :point) < 500

-- uses the index, then tests exactly
WHERE ST_DWithin(p.geom, :point, 500)

ST_DWithin is the single most valuable substitution in spatial SQL. Same answer, different order of magnitude.

3. Confirm with EXPLAIN, not with a stopwatch 

with engine.connect() as conn:
    plan = conn.exec_driver_sql(
        "EXPLAIN (ANALYZE, BUFFERS) " + sql, {"ward": "Leith", "since": "2024-01-01"}
    ).fetchall()
    print("\n".join(r[0] for r in plan))

You are reading for two words:

Seq Scan on parcels          ← no index used; 4 million rows examined
Index Scan using parcels_geom_idx   ← the index did its job

A sequential scan on a large table almost always means a missing GIST index, a predicate that cannot use one, or an SRID mismatch that forced a cast.

4. Aggregate in the database when you do not need the geometry

If the answer is a table of numbers, do not fetch geometry at all — use pandas.read_sql, not read_postgis.

import pandas as pd

totals = pd.read_sql("""
    SELECT w.name AS ward,
           count(*)                AS parcels,
           sum(ST_Area(p.geom))    AS total_m2,
           avg(ST_Area(p.geom))    AS mean_m2
    FROM parcels p
    JOIN wards w ON ST_Within(p.geom, w.geom)
    GROUP BY w.name
    ORDER BY total_m2 DESC
""", engine)

Twelve rows over the wire instead of four million. This is the biggest single win available and the one most often missed, because read_postgis is the function everyone reaches for.

5. Keep the CRS honest on the way back

read_postgis reads the SRID from the geometry and sets the GeoDataFrame's CRS from it — but only if the SRID is set. Geometries built inline in SQL often have SRID 0.

-- loses the SRID: the result has no CRS in Python
SELECT ST_Buffer(geom, 50) AS geom FROM parcels

-- keeps it
SELECT ST_SetSRID(ST_Buffer(geom, 50), ST_SRID(geom)) AS geom FROM parcels

Then assert on arrival:

gdf = gpd.read_postgis(sql, engine, geom_col="geom")
assert gdf.crs is not None, "query returned geometry with no SRID"
assert gdf.crs.to_epsg() == 27700

6. Stream results that do not fit in memory 

# read_postgis has no chunksize; go through pandas and rebuild
from shapely import from_wkb

sql = "SELECT parcel_id, ward, ST_AsBinary(geom) AS wkb FROM parcels"

for chunk in pd.read_sql(sql, engine, chunksize=50_000):
    gdf = gpd.GeoDataFrame(
        chunk.drop(columns="wkb"),
        geometry=from_wkb(chunk["wkb"]),
        crs=27700,
    )
    process(gdf)

ST_AsBinary plus shapely.from_wkb is the manual version of what read_postgis does, and it is the way to get chunking. Set crs explicitly, because ST_AsBinary drops the SRID — ST_AsEWKB keeps it if you would rather not hard-code.

Code examples

Example 1: a spatial join that stays in the database 

-- parcels tagged with their ward, plus distance to the nearest school
SELECT p.parcel_id,
       w.name AS ward,
       s.name AS nearest_school,
       ST_Distance(p.geom, s.geom) AS school_m,
       p.geom
FROM parcels p
JOIN wards w
  ON ST_Within(p.geom, w.geom)
CROSS JOIN LATERAL (
    SELECT name, geom
    FROM schools
    ORDER BY schools.geom <-> p.geom     -- index-backed nearest neighbour
    LIMIT 1
) s
WHERE w.name = %(ward)s;

The <-> operator is a KNN distance that the GIST index supports directly, so ORDER BY ... LIMIT 1 finds the nearest school without measuring all of them. There is no clean GeoPandas equivalent that is this efficient — sjoin_nearest is close, but it needs both layers in memory first.

Example 2: a reusable query layer

Inline SQL scattered through a codebase becomes unmaintainable. Give queries names and signatures.

# src/queries.py
from dataclasses import dataclass
import geopandas as gpd
import pandas as pd

PARCELS_IN_WARD = """
    SELECT p.parcel_id, w.name AS ward, ST_Area(p.geom) AS area_m2, p.geom
    FROM parcels p
    JOIN wards w ON ST_Within(p.geom, w.geom)
    WHERE w.name = %(ward)s
      AND (%(min_area)s IS NULL OR ST_Area(p.geom) >= %(min_area)s)
"""

WARD_TOTALS = """
    SELECT w.name AS ward, count(*) AS parcels, sum(ST_Area(p.geom)) AS total_m2
    FROM parcels p JOIN wards w ON ST_Within(p.geom, w.geom)
    GROUP BY w.name
"""

@dataclass
class Queries:
    engine: object

    def parcels_in_ward(self, ward: str, min_area: float | None = None) -> gpd.GeoDataFrame:
        gdf = gpd.read_postgis(
            PARCELS_IN_WARD, self.engine, geom_col="geom",
            params={"ward": ward, "min_area": min_area},
        )
        if gdf.crs is None:
            raise ValueError("query returned geometry without an SRID")
        return gdf

    def ward_totals(self) -> pd.DataFrame:
        return pd.read_sql(WARD_TOTALS, self.engine)

The (%(min_area)s IS NULL OR ...) pattern gives one query an optional filter without string concatenation. It is the small trick that stops a query module growing an f-string.

Example 3: writing results back without a round trip

When the output is derived from data already in the database, keep it there.

from sqlalchemy import text

with engine.begin() as conn:
    conn.execute(text("""
        CREATE TABLE IF NOT EXISTS parcel_buffers (
            parcel_id bigint PRIMARY KEY,
            geom geometry(Polygon, 27700) NOT NULL
        );
        TRUNCATE parcel_buffers;
        INSERT INTO parcel_buffers (parcel_id, geom)
        SELECT parcel_id, ST_SetSRID(ST_Buffer(geom, 50), 27700)
        FROM parcels
        WHERE ward = :ward;
    """), {"ward": "Leith"})

Nothing crossed the network. For bulk geometry work on data that stays in the database, this beats reading, buffering in Python and writing back — often by a factor of ten.

Explanation

Flow showing a query narrowing from the full table through the GIST bounding-box filter to exact tests and the returned rows.
The index does the cheap rejection; the exact predicate runs on what survives.

A GIST index on a geometry column stores each geometry's bounding box in an R-tree. When a query asks ST_Intersects(geom, :shape), PostGIS rewrites it internally as "bounding boxes overlap" — an index lookup that discards almost everything — followed by the exact geometric test on the remainder. That two-stage filter is why an intersects query over 40 million rows can return in milliseconds.

The rewrite only happens for predicates PostGIS knows are bounding-box-compatible. ST_Intersects, ST_Within, ST_Contains, ST_Overlaps, ST_Crosses and ST_DWithin qualify. ST_Distance(a, b) < d does not, because there is no bounding box that expresses it — hence the ST_DWithin substitution, which is the same question phrased so the index can answer it.

Two things silently disable the index:

  • An SRID mismatch. Comparing a geometry in 27700 with a literal in 4326 does not fall back to a slow query — it raises Operation on mixed SRID geometries and the query fails. Transform the literal once, in the query, so both sides match.
  • A function wrapping the indexed column. ST_Buffer(geom, 10) in a WHERE clause means the index on geom cannot help, because the thing being tested is not geom.

The wider point is about where the data has to travel. Every row a query returns must be serialised, sent over a socket, deserialised, and turned into Python objects. That per-row cost dwarfs the cost of the geometric test itself. So the fastest spatial query is almost always the one that returns the fewest rows — which is the same as saying: do the reducing work where the data already is.

Edge cases or notes

  • read_postgis needs geom_col to match your column name. The PostGIS convention is geom; GeoPandas defaults to geometry.
  • A query returning zero rows gives a GeoDataFrame with no CRS, because there was no geometry to read an SRID from. Handle empties before asserting on crs.
  • ST_Area on a geographic SRID returns square degrees. Store projected, or use geography and accept the slower spheroid maths.
  • Very long queries belong in .sql files, loaded at import. Multi-hundred-line strings inside Python are unreviewable.
  • Connection pooling matters under a scheduler. Create one engine at module level; do not build one per task.
  • LIMIT without ORDER BY is non-deterministic. For sampling, use TABLESAMPLE SYSTEM (1).
  • EXPLAIN ANALYZE runs the query. On a DELETE or UPDATE, wrap it in a transaction you roll back.
  • Keep transformation out of the join. ST_Within(p.geom, ST_Transform(w.geom, 27700)) transforms every ward for every parcel. Store both tables in one CRS instead.

FAQ

Should I use read_postgis or read_sql?

read_postgis when you need the geometry, read_sql when you need numbers. Fetching geometry you then discard is the most expensive mistake in this whole workflow.

Why is my spatial query slow even with an index?

Run EXPLAIN ANALYZE. The usual causes are a missing GIST index, ST_Distance in the WHERE clause instead of ST_DWithin, an SRID mismatch forcing a per-row transform, or a function wrapped around the indexed column.

How do I pass a geometry as a query parameter?

Send it as WKT or WKB and let PostGIS parse it: ST_GeomFromText(%(wkt)s, 27700) with params={"wkt": shape.wkt}.

Can I use GeoPandas syntax against PostGIS directly?

No — GeoPandas operates on in-memory frames. Tools like ibis or GeoAlchemy2's ORM layer let you build queries in Python, but the work still happens as SQL.

Is sjoin in GeoPandas ever better than a SQL join?

Yes, when both layers are already in memory and small, or when the join needs a predicate PostGIS lacks. Otherwise the database wins, because it never has to load either side.

How do I stop SQL injection when the table name is dynamic?

You cannot parameterise identifiers. Validate against an allow-list of known table names, or use psycopg.sql.Identifier, which quotes them safely.

Why does my buffered result lose its CRS?

ST_Buffer returns geometry with the input's SRID in most versions, but geometry constructed from scratch has SRID 0. Wrap results in ST_SetSRID(..., ST_SRID(geom)) when in doubt, and assert on gdf.crs after reading.