How to Run Spatial SQL Queries from Python with PostGIS

Problem statement

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

t0 = time.perf_counter()
parcels = gpd.read_postgis("SELECT * FROM parcels", engine, geom_col="geom")
wards = gpd.read_postgis("SELECT * FROM wards", engine, geom_col="geom")
ward = wards[wards["name"] == "Oakham North East"]
joined = gpd.sjoin(parcels, ward, predicate="within")
python_s = time.perf_counter() - t0
print(f"{len(parcels):,} parcels fetched to keep {len(joined):,}; {python_s:.2f} s")

21,253 parcels fetched to keep 1,570; 0.52 s

Every parcel in the county crossed the connection so that one ward's worth could be kept, and the spatial index the database maintains was never consulted. Rutland is small; the same pattern on a national table needs gigabytes of memory to answer a question the database could answer in milliseconds.

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:

sql = """
    SELECT p."INSPIREID" AS 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 ST_Area(p.geom) >= %(min_area)s
"""
t0 = time.perf_counter()
gdf = gpd.read_postgis(sql, engine, geom_col="geom", params={"ward": "Oakham North East", "min_area": 0})
print(f"{len(gdf):,} rows fetched in {time.perf_counter() - t0:.2f} s | CRS {gdf.crs.to_epsg()}")

1,570 rows fetched in 0.03 s | CRS 27700

The same rows, a fraction of the traffic. The rules of thumb:

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

Parameterise, index-friendly predicates, EXPLAIN

1. Parameterise, always

A ward name with an apostrophe is enough to break an f-string query — and a deliberately crafted one is enough to do worse:

ward = "Ryhall & Casterton' OR '1'='1"
unsafe = f"SELECT count(*) FROM wards WHERE name = '{ward}'"
print("f-string:", pd.read_sql(unsafe, engine).iloc[0, 0], "wards matched")
safe = pd.read_sql("SELECT count(*) FROM wards WHERE name = %(ward)s", engine, params={"ward": ward})
print("params:  ", safe.iloc[0, 0], "wards matched")

f-string: 15 wards matched
params:   0 wards matched

The f-string turned a name into a condition that matched every ward. 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:

school = "(SELECT geom FROM schools WHERE name = 'Uppingham School')"
for label, where in (("ST_Distance < 500", f"ST_Distance(p.geom, {school}) < 500"),
                     ("ST_DWithin 500", f"ST_DWithin(p.geom, {school}, 500)")):
    plan = [row[0] for row in pd.read_sql(f"EXPLAIN ANALYZE SELECT count(*) FROM parcels p WHERE {where}", engine).values]
    scan = next(line.strip().lstrip("-> ") for line in plan if "on parcels" in line)
    print(f"{label:18} {plan[-1]:28} {scan[:60]}")

ST_Distance < 500  Execution Time: 76.769 ms    Parallel Seq Scan on parcels p  (cost=0.00..158028.33 rows=4
ST_DWithin 500     Execution Time: 9.426 ms     Index Scan using idx_parcels_geom on parcels p  (cost=0.40..

Same answer, different plan: the distance filter read every parcel, ST_DWithin asked the index. ST_DWithin is the single most valuable substitution in spatial SQL.

3. Confirm with EXPLAIN, not with a stopwatch

with engine.connect() as conn:
    plan = conn.exec_driver_sql("EXPLAIN ANALYZE " + sql, {"ward": "Oakham North East", "min_area": 0}).fetchall()
print("\n".join(row[0] for row in plan if "Scan" in row[0] or "Execution" in row[0]))

  ->  Seq Scan on wards w  (cost=0.00..2.19 rows=1 width=27877) (actual time=0.004..0.005 rows=1 loops=1)
  ->  Index Scan using idx_parcels_geom on parcels p  (cost=0.28..37.56 rows=1 width=684) (actual time=0.201..4.599 rows=1570 loops=1)
Execution Time: 4.934 ms

You are reading for two words: Seq Scan means every row was examined, Index Scan means the GIST index did its job. A sequential scan on a large table almost always means a missing GIST index or a predicate that cannot use one.

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:

totals = pd.read_sql("""
    SELECT w.name AS ward,
           count(*)                        AS parcels,
           round(sum(ST_Area(p.geom)))     AS total_m2,
           round(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)
print(totals.head(4).to_string(index=False))
print(len(totals), "rows over the wire")

                     ward  parcels   total_m2  mean_m2
Braunston & Martinsthorpe     1821 43132099.0  23686.0
               Lyddington     1190 32657858.0  27444.0
                    Exton     1037 29945053.0  28877.0
                 Greetham     1079 25926059.0  24028.0
15 rows over the wire

Fifteen rows over the wire instead of twenty thousand parcels. 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. Most functions keep their input's SRID; geometries built in SQL start with none:

for label, query in (("ST_Buffer(geom, 50)", "SELECT ST_Buffer(geom, 50) AS geom FROM schools LIMIT 1"),
                     ("ST_MakePoint(x, y)", "SELECT ST_MakePoint(492000, 300000) AS geom"),
                     ("ST_SetSRID(ST_MakePoint(...), 27700)", "SELECT ST_SetSRID(ST_MakePoint(492000, 300000), 27700) AS geom")):
    print(f"{label:38} -> CRS {gpd.read_postgis(query, engine, geom_col='geom').crs}")

ST_Buffer(geom, 50)                    -> CRS EPSG:27700
ST_MakePoint(x, y)                     -> CRS None
ST_SetSRID(ST_MakePoint(...), 27700)   -> CRS EPSG:27700

Then assert on arrival: assert gdf.crs is not None catches the second case before it reaches a join. Comparing geometries in two SRIDs is an error — when PostGIS gets as far as comparing them:

for label, x, y in (("Oakham in degrees", -0.72, 52.67), ("metre-sized numbers", 492000, 300000)):
    query = f"SELECT count(*) FROM parcels WHERE ST_Intersects(geom, ST_SetSRID(ST_MakePoint({x}, {y}), 4326))"
    try:
        print(f"{label:20} -> {pd.read_sql(query, engine).iloc[0, 0]} rows, no error")
    except Exception as exc:
        message = str(exc)
        print(f"{label:20} -> {message[message.index('Operation on'):].splitlines()[0]}")

Oakham in degrees    -> 0 rows, no error
metre-sized numbers  -> Operation on mixed SRID geometries (Polygon, 27700) != (Point, 4326)

The first query is the dangerous one. Its point, in degrees, lies nowhere near any parcel's bounding box in metres, so the index filter found no candidates, the exact test never ran, and the mismatch was never noticed: zero rows, no error. Only when the numbers happen to overlap does PostGIS compare the geometries and refuse. Transform the literal to the table's SRID in the query — ST_Transform(ST_SetSRID(ST_MakePoint(...), 4326), 27700) — and never trust a zero from a query you have not checked.

6. Stream results that do not fit in memory

read_postgis accepts a chunksize and then returns an iterator of GeoDataFrames:

chunks = gpd.read_postgis('SELECT "INSPIREID", geom FROM parcels', engine, geom_col="geom", chunksize=5_000)
sizes = [len(chunk) for chunk in chunks]
print(len(sizes), "chunks of", sizes[:2], "... last", sizes[-1])

for chunk in pd.read_sql('SELECT "INSPIREID", ST_AsBinary(geom) AS wkb FROM parcels', engine, chunksize=5_000):
    first = gpd.GeoDataFrame(chunk.drop(columns="wkb"), geometry=shapely.from_wkb(chunk["wkb"]), crs=27700)
    break
print("manual WKB route:", len(first), "rows, CRS set by hand:", first.crs.to_epsg())

5 chunks of [5000, 5000] ... last 1253
manual WKB route: 5000 rows, CRS set by hand: 27700

ST_AsBinary plus shapely.from_wkb is the manual version of what read_postgis does; set crs explicitly, because ST_AsBinary drops the SRID — ST_AsEWKB keeps it.

A spatial join that stays in the database

The data used on this page

HM Land Registry's INSPIRE parcels for Rutland (© Crown copyright and database rights; OGL), Rutland's wards from the ONS December 2024 boundaries (OGL), and school sites from Geofabrik's OpenStreetMap extract of 1 January 2026 (© OpenStreetMap contributors, ODbL), loaded into a throwaway PostGIS in Docker, as other guides on this site do. to_postgis creates a GIST index on each geometry column.

import os
import re
import time
from pathlib import Path

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

Path("data").mkdir(exist_ok=True)
downloads = {"data/Rutland.zip": "https://use-land-property-data.service.gov.uk/datasets/inspire/download/Rutland.zip",
             "data/rutland.shp.zip": "https://download.geofabrik.de/europe/united-kingdom/england/rutland-260101-free.shp.zip"}
for target, url in downloads.items():
    if not Path(target).exists():
        Path(target).write_bytes(requests.get(url, timeout=300).content)
if not Path("data/wards.geojson").exists():
    Path("data/wards.geojson").write_bytes(requests.get(
        "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
        "Wards_December_2024_Boundaries_UK_BFC/FeatureServer/0/query",
        params={"where": "LAD24CD='E06000017'", "outFields": "WD24NM", "outSR": 27700, "f": "geojson"}, timeout=120).content)

os.environ.update(PGPORT="55433", PGPASSWORD="guide")
if os.system("docker inspect guide-postgis >/dev/null 2>&1") != 0:
    os.system("docker run -d --name guide-postgis -e POSTGRES_USER=user -e POSTGRES_PASSWORD=guide "
              "-e POSTGRES_DB=gis -p 127.0.0.1:55433:5432 postgis/postgis:16-3.4 >/dev/null")
    time.sleep(10)
engine = create_engine("postgresql+psycopg://user:guide@localhost:55433/gis")

layers = {
    "parcels": gpd.read_file("/vsizip/data/Rutland.zip/Land_Registry_Cadastral_Parcels.gml", columns=["INSPIREID"], DOWNLOAD_SCHEMA="NO"),
    "wards": gpd.read_file("data/wards.geojson").rename(columns={"WD24NM": "name"}),
    "schools": gpd.read_file("/vsizip/data/rutland.shp.zip/gis_osm_pois_a_free_1.shp", where="fclass = 'school'",
                             columns=["name", "fclass"]).to_crs(27700)[["name", "geometry"]],
}
for name, frame in layers.items():
    frame.rename_geometry("geom").to_postgis(name, engine, if_exists="replace", index=False)
with engine.begin() as conn:
    conn.exec_driver_sql("ANALYZE parcels; ANALYZE wards; ANALYZE schools")
    version = conn.exec_driver_sql("SELECT postgis_lib_version()").scalar()
print("PostGIS", version, {name: len(frame) for name, frame in layers.items()})

PostGIS 3.4.3 {'parcels': 21253, 'wards': 15, 'schools': 26}

Example 1: a spatial join that stays in the database

nearest_school = gpd.read_postgis("""
    SELECT p."INSPIREID" AS parcel_id,
           w.name AS ward,
           s.name AS nearest_school,
           round(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
""", engine, geom_col="geom", params={"ward": "Uppingham"})
print(nearest_school.groupby("nearest_school")["school_m"].agg(["count", "median"]).to_string())

                                            count  median
nearest_school                                           
Leighfield Primary School                      13     5.0
Uppingham Church of England Primary School    690   207.5
Uppingham Community College                    17   444.0
Uppingham School                             1298   160.0

The <-> operator is a distance the GIST index supports directly, so ORDER BY ... LIMIT 1 finds the nearest school without measuring all of them. sjoin_nearest is the GeoPandas equivalent, 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:

from dataclasses import dataclass

PARCELS_IN_WARD = """
    SELECT p."INSPIREID" AS 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::float IS NULL OR ST_Area(p.geom) >= %(min_area)s::float)
"""

@dataclass
class Queries:
    engine: object

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

q = Queries(engine)
print(len(q.parcels_in_ward("Ketton")), "parcels in Ketton;", len(q.parcels_in_ward("Ketton", min_area=10_000)), "of at least a hectare")
empty = q.parcels_in_ward("Nowhere")
print("no rows:", len(empty), "| CRS", empty.crs)

1673 parcels in Ketton; 121 of at least a hectare
no rows: 0 | CRS None

The (%(min_area)s IS NULL OR ...) pattern gives one query an optional filter without string concatenation. The ::float cast is needed: with a bare NULL, PostgreSQL cannot work out the parameter's type and refuses the query.

Example 3: writing results back without a round trip

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

with engine.begin() as conn:
    conn.execute(text("""
        DROP TABLE IF EXISTS school_buffers;
        CREATE TABLE school_buffers AS
        SELECT name, ST_Buffer(geom, 400)::geometry(Polygon, 27700) AS geom FROM schools;
    """))
    count, srid = conn.execute(text("SELECT count(*), max(ST_SRID(geom)) FROM school_buffers")).one()
print(count, "buffers written in the database, SRID", srid)

26 buffers written in the database, SRID 27700

Nothing crossed the network. For bulk geometry work on data that stays in the database, this beats reading, buffering in Python and writing back.

How the bounding-box filter works

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.

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 defeat the index or the query:

• An SRID mismatch. It raises Operation on mixed SRID geometries when the exact test runs — and when the index filter finds nothing to test, it quietly returns no rows, as step 5 showed. 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.

geom_col and empty results

• 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.

Internal links

• How to connect GeoPandas to PostGIS — engines, drivers, connection strings
• How to write a GeoDataFrame to PostGIS — getting the data in with types and indexes intact
• PostGIS explained: when a spatial database beats a folder of files — why push work down at all
• PostGIS write fails on SRID or geometry type — the SRID mismatch that also breaks queries
• Spatial indexes explained: R-trees and why spatial joins are fast — what GIST is doing underneath
• How to perform a spatial join in Python (GeoPandas) — the in-memory equivalent, and when it is the right tool
• How to find the nearest point in GeoPandas — the Python version of the <-> example
• Fixing memory errors in GeoPandas when working with large files — the symptom that sends people here

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.

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.