How to Read a Large PostGIS Table into Python Without Running Out of Memory
Problem statement
The query is fine. The read is not.
gdf = gpd.read_postgis("SELECT * FROM parcels", engine, geom_col="geom")MemoryError
Or, on a machine with more RAM, it works — and takes eleven minutes, uses 24 GB, and produces a GeoDataFrame you then immediately filter down to 40,000 rows.
The frustrating part is that chunksize does not fix it on its own:
for chunk in gpd.read_postgis(SQL, engine, geom_col="geom", chunksize=50_000):
process(chunk) # still runs out of memoryThe loop looks like streaming. It is not. By default the database driver fetches the entire result set into client memory before handing over the first chunk, so chunksize controls how the data is sliced after it has already arrived — not how much of it arrives at once.
There are four fixes, and they are worth knowing in order, because the first one usually makes the rest unnecessary.
Quick answer
1. Do not read it all. Reduce in SQL first — this solves most cases outright:
gdf = gpd.read_postgis("""
SELECT id, class, geom FROM parcels
WHERE ward_code = %(w)s
""", engine, geom_col="geom", params={"w": "E05011368"})2. Stream with a server-side cursor when you genuinely need every row:
with engine.connect().execution_options(stream_results=True) as con:
for chunk in gpd.read_postgis(text(SQL), con, geom_col="geom", chunksize=50_000):
process(chunk)stream_results=True is the argument that makes chunksize mean anything.
3. Read by spatial tile when the work is per-area rather than per-row.
4. Write to disk instead of memory — ogr2ogr streams PostGIS straight to GeoParquet without Python ever holding the data.
| Approach | Peak memory | Use when |
|---|---|---|
| reduce in SQL | result size | almost always — try this first |
| server-side cursor | one chunk | you need every row, sequentially |
| spatial tiles | one tile | the work is per-area |
| stream to disk | ~nothing | the destination is a file |
Step-by-step solution
1. Find out how big it actually is
from sqlalchemy import create_engine, text
engine = create_engine("postgresql+psycopg://user@localhost/gis")
with engine.begin() as con:
stats = con.execute(text("""
SELECT
(SELECT COUNT(*) FROM parcels) AS rows,
pg_size_pretty(pg_total_relation_size('parcels')) AS on_disk,
pg_size_pretty(SUM(ST_MemSize(geom))::bigint) AS geom_bytes,
ROUND(AVG(ST_NPoints(geom))) AS avg_vertices,
MAX(ST_NPoints(geom)) AS max_vertices
FROM parcels
""")).mappings().one()
print(dict(stats)){'rows': 4012884, 'on_disk': '6284 MB', 'geom_bytes': '4917 MB',
'avg_vertices': 118, 'max_vertices': 402113}
Two numbers decide the strategy. geom_bytes is the wire payload — and in Python it is worse, because each geometry becomes a Shapely object with per-object overhead, typically two to three times the WKB size. And max_vertices of 402,113 means a single feature is around 6 MB on its own; a chunk that happens to contain a few of those is much bigger than the average suggests.
2. Reduce in SQL before anything else
This is not a fallback. It is the answer for most cases:
# spatial filter — uses the GiST index
SQL = """
SELECT id, class, geom
FROM parcels
WHERE ST_Intersects(geom, ST_MakeEnvelope(380000, 395000, 400000, 410000, 27700))
"""
# columns you use, not SELECT *
SQL = "SELECT id, class, geom FROM parcels WHERE class = 'residential'"
# derived values instead of geometry, when the shape is not needed
SQL = "SELECT id, class, ST_Area(geom) AS area_m2 FROM parcels"
# simplified geometry when full detail will not be used
SQL = "SELECT id, ST_SimplifyPreserveTopology(geom, 5) AS geom FROM parcels"That third one is the biggest lever and the most often missed. If the analysis needs areas and centroids rather than shapes, computing them in SQL turns 4.9 GB of geometry into 60 MB of floats.
The fourth is the next best: on parcel data, simplifying to a 5 m tolerance typically removes 80–90% of vertices while leaving the shapes visually identical.
3. Stream properly with a server-side cursor
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 without buffering the whole result client-side."""
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, chunksize=chunk, params=params or {}
)
total_area, rows = 0.0, 0
for part in stream_postgis("SELECT id, geom FROM parcels", engine):
total_area += part.area.sum()
rows += len(part)
print(f" {rows:>10,} rows {total_area/1e6:>12,.1f} km²") 50,000 rows 1,204.3 km²
100,000 rows 2,418.8 km²
…
4,012,884 rows 94,882.1 km²
Peak memory stays at roughly one chunk regardless of table size. Three details make it work:
stream_results=Trueswitches psycopg to a named (server-side) cursor. Without it, the driver buffers everything before yielding anything.- The connection must stay open for the whole iteration. Using
enginedirectly rather than a connection closes it between chunks and defeats the cursor. text(sql)plusparamskeeps values bound rather than interpolated. String formatting into SQL is an injection risk even against your own database.
4. Read by spatial tile when the work is per-area
Sequential streaming is wrong when the operation needs neighbours — a dissolve, a nearest-neighbour search, anything with spatial reach. Chunk by geography instead:
import geopandas as gpd
from sqlalchemy import text
def tile_bounds(engine, table, size, geom="geom"):
with engine.begin() as con:
minx, miny, maxx, maxy = con.execute(text(f"""
SELECT ST_XMin(e), ST_YMin(e), ST_XMax(e), ST_YMax(e)
FROM (SELECT ST_Extent({geom}) AS e FROM {table}) s
""")).one()
y = miny
while y < maxy:
x = minx
while x < maxx:
yield (x, y, min(x + size, maxx), min(y + size, maxy))
x += size
y += size
def read_tile(engine, table, bounds, srid=27700, buffer=0, geom="geom"):
sql = f"""
SELECT * FROM {table}
WHERE ST_Intersects({geom},
ST_Expand(ST_MakeEnvelope(:x0, :y0, :x1, :y1, :srid), :buf))
"""
return gpd.read_postgis(text(sql), engine, geom_col=geom, params={
"x0": bounds[0], "y0": bounds[1], "x1": bounds[2], "y1": bounds[3],
"srid": srid, "buf": buffer,
})
for i, bounds in enumerate(tile_bounds(engine, "parcels", 10_000)):
gdf = read_tile(engine, "parcels", bounds, buffer=500)
if gdf.empty:
continue
result = expensive_operation(gdf)
# keep only features whose representative point is in the core tile
core = gpd.GeoSeries.from_wkt([f"POLYGON(({bounds[0]} {bounds[1]}, …))"])
print(f" tile {i}: {len(gdf):,} read, {len(result):,} kept")The buffer=500 reads 500 m beyond the tile so neighbour-dependent operations see the features they need, and the core filter then keeps each output feature exactly once. That read-wide, write-narrow pattern is the same one used when splitting a large layer into tiles.
The index does the work: ST_Intersects against a small envelope is an index scan, so reading tile 400 costs the same as reading tile 1.
5. Skip Python entirely when the destination is a file
If the goal is a GeoPackage or a Parquet file, Python is an expensive middleman:
import subprocess
def export(dsn, sql, out_path, fmt="Parquet"):
cmd = ["ogr2ogr", "-f", fmt, str(out_path), f"PG:{dsn}",
"-sql", sql, "-progress"]
if fmt == "Parquet":
cmd += ["-lco", "COMPRESSION=ZSTD", "-lco", "GEOMETRY_ENCODING=WKB"]
subprocess.run(cmd, check=True)
export("host=localhost dbname=gis user=gis",
"SELECT id, class, geom FROM parcels", "parcels.parquet")This streams row by row from the database to the file. Peak memory is a few megabytes regardless of table size, and the result can then be read back lazily — see GeoParquet and columnar storage explained.
Code examples
Example 1: streaming with progress, retries and a bounded memory footprint
import time
import geopandas as gpd
import pandas as pd
from sqlalchemy import create_engine, text
def process_large_table(engine, sql, fn, *, geom_col="geom", chunk=50_000,
params=None, expected=None):
"""Apply `fn` to each chunk, accumulating small results only."""
if expected is None:
with engine.begin() as con:
expected = con.execute(
text(f"SELECT COUNT(*) FROM ({sql}) s"), params or {}
).scalar()
results, seen, t0 = [], 0, time.perf_counter()
with engine.connect().execution_options(
stream_results=True, yield_per=chunk
) as con:
for part in gpd.read_postgis(text(sql), con, geom_col=geom_col,
chunksize=chunk, params=params or {}):
results.append(fn(part)) # must return something SMALL
seen += len(part)
rate = seen / (time.perf_counter() - t0)
eta = (expected - seen) / rate if rate else 0
print(f" {seen:>10,} / {expected:,} "
f"{rate:>8,.0f} rows/s eta {eta/60:>5.1f} min", end="\r")
print()
return pd.concat(results, ignore_index=True)
def summarise(chunk):
"""One row per class — small enough to accumulate."""
return (chunk.assign(area=chunk.area)
.groupby("class", as_index=False)
.agg(n=("area", "size"), area_m2=("area", "sum")))
summary = process_large_table(
engine, "SELECT id, class, geom FROM parcels", summarise
).groupby("class", as_index=False).sum()
print(summary) 4,012,884 / 4,012,884 18,204 rows/s eta 0.0 min
class n area_m2
0 agricultural 412884 8.412993e+09
1 residential 2841002 1.204118e+09
2 commercial 758998 4.028841e+08
The load-bearing constraint is in the docstring: fn must return something small. Appending each chunk to a list defeats the entire exercise — you end up holding the whole table in results instead of in one DataFrame. Aggregating per chunk and combining the aggregates keeps memory flat, and works for any associative statistic. Medians and percentiles cannot be combined this way and need a different approach.
The counting query wraps the user's SQL in a subquery. On a table this size that count is itself a scan, so pass expected= when you already know it.
Example 2: reading only what a spatial filter selects
The most effective optimisation is often just a WHERE clause the index can use:
import geopandas as gpd
from shapely.geometry import box
from sqlalchemy import create_engine, text
engine = create_engine("postgresql+psycopg://user@localhost/gis")
def read_within(engine, table, geom, *, srid=27700, columns="*",
simplify=None, geom_col="geom"):
"""Read only the features intersecting `geom`, optionally simplified."""
geom_expr = (f"ST_SimplifyPreserveTopology({geom_col}, :tol) AS {geom_col}"
if simplify else geom_col)
cols = columns if columns == "*" else ", ".join(columns)
sql = f"""
SELECT {cols if columns == '*' else cols + ','} {geom_expr}
FROM {table}
WHERE ST_Intersects({geom_col}, ST_GeomFromText(:wkt, :srid))
"""
params = {"wkt": geom.wkt, "srid": srid}
if simplify:
params["tol"] = simplify
return gpd.read_postgis(text(sql), engine, geom_col=geom_col, params=params)
area = box(380_000, 395_000, 400_000, 410_000)
full = read_within(engine, "parcels", area, columns=["id", "class"])
simple = read_within(engine, "parcels", area, columns=["id", "class"], simplify=5)
for name, gdf in [("full detail", full), ("simplified 5 m", simple)]:
mb = gdf.memory_usage(deep=True).sum() / 1e6
verts = gdf.geometry.apply(lambda g: len(g.exterior.coords)
if g.geom_type == "Polygon" else 0).sum()
print(f"{name:<16} {len(gdf):>8,} rows {mb:>7.1f} MB {verts:>10,} vertices")full detail 178,204 271.9 MB 21,028,104 vertices
simplified 5 m 178,204 41.2 MB 2,884,117 vertices
Same rows, one-seventh the memory. Simplification is lossy, so it is right for display, screening and area-scale statistics, and wrong when the exact boundary matters — a cadastral check or a legal boundary. ST_SimplifyPreserveTopology at least guarantees the result stays valid, unlike plain ST_Simplify.
Example 3: diagnosing where the memory actually goes
import gc, tracemalloc
import geopandas as gpd
from sqlalchemy import create_engine, text
def measure(sql, engine, geom_col="geom", label=""):
gc.collect()
tracemalloc.start()
gdf = gpd.read_postgis(text(sql), engine, geom_col=geom_col)
current, peak = tracemalloc.get_traced_memory()
tracemalloc.stop()
attrs = gdf.drop(columns=[geom_col]).memory_usage(deep=True).sum()
geoms = gdf.geometry.memory_usage(deep=True)
verts = int(gdf.geometry.apply(lambda g: 0 if g is None else len(g.wkb)).sum())
print(f"{label}")
print(f" rows {len(gdf):,}")
print(f" attributes {attrs/1e6:8.1f} MB")
print(f" geometry {geoms/1e6:8.1f} MB (WKB on the wire: {verts/1e6:.1f} MB)")
print(f" peak alloc {peak/1e6:8.1f} MB")
del gdf
gc.collect()
measure("SELECT * FROM parcels WHERE ward_code = 'E05011368'", engine,
label="everything")
measure("SELECT id, class, geom FROM parcels WHERE ward_code = 'E05011368'", engine,
label="three columns")everything
rows 2,847
attributes 18.4 MB
geometry 41.2 MB (WKB on the wire: 14.8 MB)
peak alloc 187.3 MB
three columns
rows 2,847
attributes 0.4 MB
geometry 41.2 MB (WKB on the wire: 14.8 MB)
peak alloc 104.6 MB
Two things are worth reading carefully. The geometry occupies 41 MB in memory against 15 MB on the wire — the roughly 3× Shapely object overhead, which is why "the table is only 5 GB" underestimates what the read costs. And peak allocation is far above the final size, because the WKB buffer and the parsed objects coexist during parsing. Sizing a read against final DataFrame size will under-provision by a factor of two or more.
Explanation
Reading a large table into Python fails for a reason that has nothing to do with PostGIS and everything to do with where the data has to be materialised.
The database is built to avoid materialising. Postgres streams a result set: it produces rows as the plan yields them, holds a bounded amount in its own buffers, and can serve a query against a table far larger than its memory. A DataFrame is the opposite — it is a materialised, in-memory, random-access structure, and constructing one requires every row to exist simultaneously.
The default cursor is what turns a streaming source into a materialised one. By default psycopg uses a client-side cursor: it asks for the whole result, buffers it locally, and only then lets you iterate. Adding chunksize slices that already-complete buffer. stream_results=True switches to a named cursor, where the server keeps the result set and hands over batches on request — and only then does chunksize describe how much is in memory. This is why the "streaming" loop that still runs out of memory is such a common report: the code reads as if it streams and the plumbing underneath does not.
Shapely object overhead is the second surprise. WKB is a compact binary format; a Shapely geometry is a Python object wrapping a GEOS structure, with per-object overhead that is significant for small geometries and unavoidable for all of them. A table whose geometry is 5 GB on the wire routinely becomes 12–15 GB in memory. Any capacity planning based on pg_total_relation_size will be wrong in the dangerous direction.
Which is why reduction beats streaming. Streaming makes an impossible read possible; not reading beats both. Every one of the four reduction techniques — fewer rows, fewer columns, derived values instead of geometry, simplified geometry — attacks the payload before it is ever serialised. Computing ST_Area(geom) in SQL replaces a 6 MB polygon with an 8-byte float, and no amount of clever client-side memory management competes with that.
Finally, the shape of the chunking should match the shape of the work. Sequential chunks are correct for row-independent operations: reprojecting, computing an attribute, filtering. They are wrong for anything with spatial reach, because a dissolve or a nearest-neighbour search needs features that arbitrary row order will have put in a different chunk. Spatial tiling with a buffer is the fix, and its correctness argument — read wide, write narrow, assign each feature to exactly one tile — is identical to the one for tiling a file-based layer. The database is a different source; the geometry of the problem is the same.
Edge cases or notes
chunksizewithoutstream_results=Truedoes not stream. The whole result is buffered client-side first.- The connection must stay open across the iteration. Passing
enginerather than a connection closes it between chunks. - A server-side cursor holds a transaction open, which blocks
VACUUMon the table. Do not leave one open for hours. - Shapely objects cost roughly 2–3× their WKB size. Size reads against that, not against the table size.
ST_Simplifycan produce invalid geometry;ST_SimplifyPreserveTopologywill not.- A single huge geometry can blow a chunk budget. Check
MAX(ST_NPoints(geom))before choosing a chunk size. ST_AsBinarydrops the SRID;ST_AsEWKBkeeps it.read_postgishandles this, but hand-rolled readers often do not.work_memandstatement_timeoutare server settings that can abort a long streaming read; check them for scheduled jobs.COPY (SELECT …) TO STDOUTis the fastest bulk export path and is whatogr2ogruses withPG_USE_COPY.del gdf; gc.collect()actually matters in a loop — a lingering reference keeps a multi-gigabyte frame alive.
Internal links
- Spatial SQL or GeoPandas? — deciding what to read in the first place
- How to run spatial SQL queries from Python with PostGIS — the basics of
read_postgis - PostGIS spatial indexes explained — what makes a filtered read cheap
- Fixing memory errors in GeoPandas when working with large files — the file-based version
- How to process a very large GeoPackage in chunks — the same technique on files
- How to split a large layer into tiles — spatial chunking, in depth
- GeoParquet and columnar storage explained — a better destination than a DataFrame
- My PostGIS spatial query is slow — when the query, not the read, is the problem
FAQ
Why does chunksize not stop the MemoryError?
Because the driver buffers the entire result set client-side before yielding the first chunk. Add stream_results=True on the connection to switch to a server-side cursor.
How much memory will a table need in Python?
Roughly two to three times its WKB size, because Shapely objects carry per-object overhead. Peak allocation during parsing is higher still, since the buffer and the objects coexist.
What is the single most effective fix?
Reduce in SQL. Fewer columns, a spatial filter, or replacing geometry with derived values such as ST_Area(geom) — the last can cut the payload by three orders of magnitude.
Should I chunk by rows or by area?
By rows for row-independent work. By spatial tile when the operation needs neighbours — dissolves, nearest-neighbour searches — because arbitrary row order scatters neighbours across chunks.
Is it safe to leave a server-side cursor open?
Not for long. It holds a transaction open, which blocks VACUUM and can bloat the table. Iterate promptly and close the connection.
Can I avoid Python entirely?
Yes, when the destination is a file. ogr2ogr -f Parquet out.parquet PG:… -sql "…" streams from the database to disk with a memory footprint of a few megabytes.
Does simplifying geometry lose data?
Yes — it removes vertices. Use it for display, screening and area-scale statistics; not for anything where the exact boundary matters. ST_SimplifyPreserveTopology at least guarantees validity.