How to Read a Large PostGIS Table into Python Without Running Out of Memory

Problem statement

The query is fine. The read is not. Here are the 293,729 HM Land Registry parcels for Leeds — a 131 MB table — read into GeoPandas three ways, each in a fresh Python process so its peak memory can be measured:

import subprocess
import sys

PROBE = '''
import sys
import geopandas as gpd
from sqlalchemy import create_engine, text


def memory(field):                      # this process only: a fresh exec starts from zero
    with open("/proc/self/status") as f:
        return next(int(line.split()[1]) for line in f if line.startswith(field)) // 1024


engine = create_engine("postgresql+psycopg://user@localhost/gis")
mode, chunk, sql = sys.argv[1], int(sys.argv[2]), sys.argv[3]
rows, wkb, keep = 0, 0, []
if mode == "whole table":
    gdf = gpd.read_postgis(sql, engine, geom_col="geom")
    rows, wkb, keep = len(gdf), int(gdf.geometry.to_wkb().map(len).sum()), [gdf]
elif mode == "chunksize only":
    for part in gpd.read_postgis(sql, engine, geom_col="geom", chunksize=chunk):
        rows, wkb = rows + len(part), wkb + int(part.geometry.to_wkb().map(len).sum())
elif mode == "server-side cursor":
    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", chunksize=chunk):
            rows, wkb = rows + len(part), wkb + int(part.geometry.to_wkb().map(len).sum())
print(rows, wkb // 1_000_000, memory("VmRSS:"), memory("VmHWM:"))
'''
with open("probe.py", "w") as f:
    f.write(PROBE)


def peak_mb(mode, chunk=50_000, sql="SELECT id, geom FROM parcels"):
    rows, wkb_mb, rss, peak = map(int, subprocess.run(
        [sys.executable, "probe.py", mode, str(chunk), sql],
        capture_output=True, text=True, check=True).stdout.split())
    return rows, peak, wkb_mb, rss


baseline = peak_mb("imports only")[1]
for mode in ("whole table", "chunksize only"):
    rows, peak, _, _ = peak_mb(mode)
    print(f"{mode:15} {rows:,} rows   peak {peak:5,} MB  ({peak - baseline:,} MB above a bare interpreter)")

whole table     293,729 rows   peak 1,314 MB  (1,149 MB above a bare interpreter)
chunksize only  293,729 rows   peak   722 MB  (557 MB above a bare interpreter)

Reading the whole table costs over a gigabyte for a 131 MB table. And chunksize does not fix it on its own — it halves the peak here, but no more:

for chunk in gpd.read_postgis(SQL, engine, geom_col="geom", chunksize=50_000):
    process(chunk)          # still buffers the whole result first

The 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 arrived — not how much of it arrives at once. On a table ten or fifty times this size, that is the MemoryError.

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(text("""
    SELECT id, ward_code, geom FROM parcels
    WHERE ward_code = :w
"""), engine, geom_col="geom", params={"w": "E05011397"})       # Headingley & Hyde Park
print(f"{len(gdf):,} parcels")

7,337 parcels

2. Stream with a server-side cursor when you genuinely need every row:

with engine.connect().execution_options(stream_results=True, yield_per=5_000) as con:
    for chunk in gpd.read_postgis(text(SQL), con, geom_col="geom", chunksize=5_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.

Measure the table, reduce in SQL, stream properly

1. Find out how big it actually is

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': 293729, 'on_disk': '131 MB', 'geom_bytes': '113 MB', 'avg_vertices': Decimal('23'), 'max_vertices': 5526}

Two numbers decide the strategy. geom_bytes is the wire payload — and in Python it grows, because each geometry becomes a Shapely object with per-object overhead (Example 3 measures it). And max_vertices says whether one feature can blow a chunk budget: a chunk that happens to contain the largest features 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:

VARIANTS = {
    "full geometry":     ("SELECT id, geom FROM parcels", True),
    "area, no geometry": ("SELECT id, ST_Area(geom) AS area_m2 FROM parcels", False),
    "simplified to 5 m": ("SELECT id, ST_SimplifyPreserveTopology(geom, 5) AS geom FROM parcels", True),
    "one ward":          ("SELECT id, geom FROM parcels WHERE ward_code = 'E05011397'", True),
}
for name, (sql, has_geom) in VARIANTS.items():
    vertices = "SUM(ST_NPoints(v.geom))" if has_geom else "0"
    size = q(f"SELECT COUNT(*) AS n, SUM(pg_column_size(v.*)) AS bytes, {vertices} AS vertices FROM ({sql}) v")
    vertex_count = round(int(size.vertices[0]), -3)    # simplifying gives a few more or fewer each run
    print(f"{name:18} {size.n[0]:>8,} rows  {size.bytes[0] / 1e6:6.1f} MB  {vertex_count:>10,} vertices")

full geometry       293,729 rows   116.0 MB   6,638,000 vertices
area, no geometry   293,729 rows    11.7 MB           0 vertices
simplified to 5 m   293,729 rows    46.9 MB   1,713,000 vertices
one ward              7,337 rows     2.1 MB     104,000 vertices

Vertex counts are printed to the thousand because the simplified one is not repeatable: three runs on the same table gave 1,712,689, 1,712,691 and 1,712,692. The area-only query is the biggest lever and the most often missed: if the analysis needs areas and centroids rather than shapes, computing them in SQL leaves the geometry in the database. Simplification is next — lossy, and so right for display and screening, wrong for anything where the exact boundary matters.

3. Stream properly with a server-side cursor

def stream_postgis(sql, engine, *, geom_col="geom", chunk=5_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:,} rows, {total_area / 1e6:,.1f} km² of parcels")

for chunk in (5_000, 20_000, 50_000):
    streamed = peak_mb("server-side cursor", chunk)[1]
    sliced = peak_mb("chunksize only", chunk)[1]
    print(f"chunks of {chunk:>6,}: server-side cursor {streamed - baseline:4,} MB,  "
          f"chunksize only {sliced - baseline:4,} MB above baseline")

293,729 rows, 587.8 km² of parcels
chunks of  5,000: server-side cursor   77 MB,  chunksize only  315 MB above baseline
chunks of 20,000: server-side cursor  177 MB,  chunksize only  405 MB above baseline
chunks of 50,000: server-side cursor  346 MB,  chunksize only  558 MB above baseline

With a server-side cursor, peak memory follows the chunk size. Without one it cannot fall below the whole result, however small the chunks. Three details make it work:

• stream_results=True switches psycopg to a named (server-side) cursor. Without it, the driver buffers everything before yielding anything.
• Pass the connection that carries the option. stream_results is set on the connection's execution options; handing read_postgis a plain engine reads with an ordinary cursor, which is the "chunksize only" row.
• text(sql) plus params keeps 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, reading a buffer around each tile and keeping only the features whose representative point falls in the tile's core:

from shapely.geometry import box


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,
    })


kept_total, read_total, tiles = 0, 0, 0
for bounds in tile_bounds(engine, "parcels", 10_000):
    gdf = read_tile(engine, "parcels", bounds, buffer=500)
    if gdf.empty:
        continue
    tiles += 1
    read_total += len(gdf)
    x0, y0, x1, y1 = bounds
    points = gdf.geometry.representative_point()
    core = (points.x >= x0) & (points.x < x1) & (points.y >= y0) & (points.y < y1)
    kept_total += int(core.sum())
print(f"{tiles} tiles, {read_total:,} features read, {kept_total:,} kept once each")

12 tiles, 357,013 features read, 293,729 kept once each

The 500 m buffer reads features beyond each tile so that neighbour-dependent operations see what they need; the core filter then keeps each feature exactly once, which is why the kept total equals the table's row count while the read total is larger. 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 every tile costs about the same.

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 os


def export(dsn, sql, out_path, fmt="GPKG"):
    cmd = ["ogr2ogr", "-f", fmt, str(out_path), f"PG:{dsn}", "-sql", sql, "-nln", "parcels", "-overwrite"]
    subprocess.run(cmd, check=True)


export("host=localhost dbname=gis user=user", "SELECT id, ward_code, geom FROM parcels", "parcels.gpkg")
import pyogrio
print(f"parcels.gpkg: {pyogrio.read_info('parcels.gpkg')['features']:,} rows, "
      f"{os.path.getsize('parcels.gpkg') / 1e6:.1f} MB")

parcels.gpkg: 293,729 rows, 160.5 MB

ogr2ogr reads through a server-side cursor in batches and writes as it goes, so its memory does not grow with the table. With a GDAL built with the Arrow/Parquet driver, -f Parquet writes GeoParquet the same way, and the result can be read back lazily — see GeoParquet and columnar storage explained.

Streaming with bounded memory

The data used on this page

HM Land Registry's INSPIRE Index Polygons for Leeds, loaded into a local PostGIS 3.4 (the postgis/postgis:16-3.4 Docker image; connection details from PGPORT and PGPASSWORD). Each parcel is tagged in SQL with the ONS December 2024 ward its representative point falls in, so there is an attribute to filter and group by.

import io
import os
import time

import geopandas as gpd
import pandas as pd
import requests
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, **params):
    return pd.read_sql(text(sql), engine, params=params)


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", "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"}).rename_geometry("geom")

with engine.begin() as con:
    con.execute(text("DROP TABLE IF EXISTS parcels, wards_tag CASCADE"))
parcels.to_postgis("parcels", engine, index=False, chunksize=50_000)
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)"""))
    con.execute(text("CREATE INDEX ON parcels (ward_code)"))
    con.execute(text("ANALYZE parcels"))
with engine.connect().execution_options(isolation_level="AUTOCOMMIT") as con:
    con.execute(text("VACUUM FULL parcels"))            # the UPDATE left a dead copy of every row
print(q("SELECT COUNT(*) AS parcels, COUNT(ward_code) AS with_ward FROM parcels").to_string(index=False))

 parcels  with_ward
  293729     293327

Example 1: streaming with progress and a bounded memory footprint

import time


def process_large_table(engine, sql, fn, *, geom_col="geom", chunk=5_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)
    print(f"{seen:,} / {expected:,} rows at {rate:,.0f} rows/s")
    return pd.concat(results, ignore_index=True)


def summarise(chunk):
    """One row per ward — small enough to accumulate."""
    return (chunk.assign(area=chunk.area)
                 .groupby("ward_code", dropna=False, as_index=False)
                 .agg(n=("area", "size"), area_m2=("area", "sum")))


summary = (process_large_table(engine, "SELECT id, ward_code, geom FROM parcels", summarise)
           .groupby("ward_code", dropna=False, as_index=False).sum()
           .sort_values("n", ascending=False))
print(summary.head(3).to_string(index=False))
print(f"... {len(summary)} groups, {summary['n'].sum():,} parcels")

293,729 / 293,729 rows at 66,084 rows/s
ward_code     n      area_m2
E05011406 11152 1.244274e+07
E05011395 10787 1.943988e+07
E05011407 10738 1.240526e+07
... 34 groups, 293,729 parcels

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 large table 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:

def read_within(engine, table, geom, *, srid=27700, columns=("id",),
                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)
    sql = f"""
        SELECT {', '.join(columns)}, {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(428_000, 432_000, 432_000, 436_000)             # 4 km around Leeds city centre
full = read_within(engine, "parcels", area, columns=["id", "ward_code"])
simple = read_within(engine, "parcels", area, columns=["id", "ward_code"], simplify=5)

for name, gdf in [("full detail", full), ("simplified 5 m", simple)]:
    mb = gdf.memory_usage(deep=True).sum() / 1e6
    verts = round(int(gdf.geometry.count_coordinates().sum()), -3)
    print(f"{name:<16} {len(gdf):>8,} rows  {mb:>7.1f} MB  {verts:>10,} vertices")

full detail        20,173 rows      0.7 MB     388,000 vertices
simplified 5 m     20,173 rows      0.7 MB     117,000 vertices

Same rows, fewer vertices. 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. Note what does not shrink: memory_usage counts a Shapely geometry as a pointer, so a DataFrame's reported size hides most of what its geometries cost — which is what Example 3 measures.

Example 3: diagnosing where the memory actually goes

The probe from the problem statement reports, for a fresh process, the WKB that crossed the wire, the memory the process still holds after the read, and its peak:

for label, sql in [("one ward", "SELECT id, geom FROM parcels WHERE ward_code = 'E05011397'"),
                   ("all of Leeds", "SELECT id, geom FROM parcels")]:
    rows, peak, wkb_mb, rss = peak_mb("whole table", sql=sql)
    print(f"{label:13} {rows:>8,} rows   WKB {wkb_mb:4,} MB   held after the read {rss - baseline:5,} MB   "
          f"peak {peak - baseline:5,} MB   (above a bare interpreter)")

one ward         7,337 rows   WKB    1 MB   held after the read    38 MB   peak    38 MB   (above a bare interpreter)
all of Leeds   293,729 rows   WKB  110 MB   held after the read 1,076 MB   peak 1,148 MB   (above a bare interpreter)

Two things are worth reading carefully. What the process holds is a multiple of the WKB that crossed the wire — Shapely objects, their GEOS coordinate arrays, the GeoDataFrame and its index all cost more than the compact binary they came from, which is why "the table is only 131 MB" underestimates the read. And the peak is well above what is held, because the driver's copy of the raw result and the parsed objects coexist during parsing. Python's own tracemalloc cannot see most of this — GEOS and the database driver allocate in C — so measure at the process level. Size a read against measured peak memory, not against the table or the final DataFrame.

Where the data gets materialised

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. Leeds's 110 MB of WKB became over a gigabyte held by the process that read it. 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 polygon of up to 5,526 vertices 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.

chunksize without streaming and server cursors

• chunksize without stream_results=True does not stream. The whole result is buffered client-side first.
• Iterate inside the with block. The server-side cursor lives as long as the connection that opened it.
• A server-side cursor holds a transaction open, which blocks VACUUM on the table. Do not leave one open for hours.
• A read costs far more than its WKB. Here 110 MB on the wire was about a gigabyte in the process. Size reads against measured memory, not the table size.
• ST_Simplify can produce invalid geometry; ST_SimplifyPreserveTopology will not.
• A single huge geometry can blow a chunk budget. Check MAX(ST_NPoints(geom)) before choosing a chunk size.
• ST_AsBinary drops the SRID; ST_AsEWKB keeps it. read_postgis handles this, but hand-rolled readers often do not.
• work_mem and statement_timeout are server settings that can abort a long streaming read; check them for scheduled jobs.
• COPY (SELECT …) TO STDOUT is the fastest bulk export path from the server. PG_USE_COPY is the other direction — it makes ogr2ogr load with 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 GPKG out.gpkg PG:… -sql "…" — or -f Parquet with a GDAL that has the Parquet driver — streams from the database to disk in batches, so its memory does not grow with the table.

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.