How to Process a Very Large GeoPackage in Chunks with Python

Problem statement

One file, 34 million buildings, 41 GB. The obvious line does not work:

>>> gdf = gpd.read_file("data/raw/buildings.gpkg")
Killed

There is no bad file to skip and no loop to fix — the dataset simply does not fit in memory. And the usual batch pattern does not apply either, because batching assumes many files; here there is one.

Chunking is the answer: read a slice, process it, write it, release it, repeat. Peak memory then depends on the chunk size rather than the file size, and the job becomes resumable and parallelisable almost for free.

The complications are real, though:

some operations are row-independent (buffer, reproject, attribute maths) and chunk trivially
some need neighbours (dissolve, spatial join, nearest) and need spatial chunking instead
some need the whole dataset (global statistics, deduplication) and need two passes
writing from many chunks into one output has its own rules

Quick answer

Read by row window with pyogrio, append to the output, and never hold two chunks at once:

get the feature count from the header — do not read the data to count it
loop skip_features / max_features in windows of 50k–250k rows
process the chunk and append it to the output (mode="a" after the first write)
delete the chunk and report progress
for neighbour-dependent work, chunk by space with bbox= instead of by row

from pathlib import Path
import geopandas as gpd
import pyogrio

SRC = "data/raw/buildings.gpkg"
OUT = Path("data/out/buildings_area.gpkg")
CHUNK = 100_000

info = pyogrio.read_info(SRC)
total = info["features"]
print(f"{total:,} features, CRS {info['crs']}")

OUT.parent.mkdir(parents=True, exist_ok=True)
OUT.unlink(missing_ok=True)

for start in range(0, total, CHUNK):
    gdf = gpd.read_file(SRC, engine="pyogrio", skip_features=start, max_features=CHUNK,
                        columns=["id", "type", "geometry"])
    gdf["area_m2"] = gdf.to_crs(gdf.estimate_utm_crs()).area
    gdf.to_file(OUT, layer="buildings", driver="GPKG", mode="w" if start == 0 else "a")
    print(f"[{min(start + CHUNK, total):>10,}/{total:,}] {(start + len(gdf))/total:5.1%}", flush=True)
    del gdf

read_info() is a header read, so the feature count costs nothing. columns= is the other half of the win — narrowing the read at the driver means the unwanted columns never become Python objects.

Row chunks or spatial chunks

Grid comparing row-window, bbox-tile, attribute-partition and SQL-filter chunking. — Four ways to cut a dataset, and which operations each one supports.

Step-by-step solution

Vertical steps: read info, plan chunks, read window, transform, append, release. — Six steps, and the last one is what keeps memory flat.

Inspect before you read

import pyogrio

info = pyogrio.read_info("data/raw/buildings.gpkg")
print(f"layers      : {[name for name, _ in pyogrio.list_layers('data/raw/buildings.gpkg')]}")
print(f"features    : {info['features']:,}")
print(f"geometry    : {info['geometry_type']}")
print(f"crs         : {info['crs']}")
print(f"fields      : {len(info['fields'])} — {list(info['fields'])[:8]}")
print(f"bounds      : {info['total_bounds']}")

Feature count, schema and extent all come from the header. That is enough to plan the whole job without reading a single geometry.

Choose a chunk size from memory, not from habit

import pyogrio, geopandas as gpd

sample = gpd.read_file(SRC, engine="pyogrio", max_features=5_000,
                       columns=["id", "type", "geometry"])
bytes_per_feature = sample.memory_usage(deep=True).sum() / len(sample)
print(f"~{bytes_per_feature:,.0f} bytes per feature")

BUDGET_GB = 2
chunk = int(BUDGET_GB * 1e9 / (bytes_per_feature * 3))     # ×3 headroom for intermediates
print(f"chunk size ≈ {chunk:,} features")

Measuring a small sample beats guessing. The factor of three is for the temporary copies that operations like to_crs and buffer create.

Chunk by row window for row-independent work

def iter_row_chunks(path, layer=None, chunk=100_000, columns=None, where=None):
    """Yield successive row windows of a layer as GeoDataFrames."""
    total = pyogrio.read_info(path, layer=layer)["features"]
    for start in range(0, total, chunk):
        gdf = gpd.read_file(path, layer=layer, engine="pyogrio",
                            skip_features=start, max_features=chunk,
                            columns=columns, where=where)
        if gdf.empty:
            break
        yield start, total, gdf

for start, total, gdf in iter_row_chunks(SRC, chunk=100_000, columns=["id", "type", "geometry"]):
    ...

Reprojection, attribute maths, buffering, validity repair, format conversion and filtering are all row-independent — each feature's result depends only on itself, so the window boundary is irrelevant.

Chunk by space when neighbours matter

A dissolve or a spatial join spanning a window boundary gives the wrong answer. Tile the extent instead, and overlap the tiles.

import geopandas as gpd
import pyogrio

info = pyogrio.read_info(SRC)
minx, miny, maxx, maxy = info["total_bounds"]

NX = NY = 8
OVERLAP = 250          # CRS units — at least the largest feature you expect

dx, dy = (maxx - minx) / NX, (maxy - miny) / NY
for i in range(NX):
    for j in range(NY):
        bbox = (minx + i*dx - OVERLAP, miny + j*dy - OVERLAP,
                minx + (i+1)*dx + OVERLAP, miny + (j+1)*dy + OVERLAP)
        tile = gpd.read_file(SRC, engine="pyogrio", bbox=bbox)
        if tile.empty:
            continue

        result = gpd.sjoin(tile, zones, predicate="within")

        # keep only results whose representative point is in the *core* tile,
        # so overlapping reads do not produce duplicates
        core = (minx + i*dx, miny + j*dy, minx + (i+1)*dx, miny + (j+1)*dy)
        pts = result.representative_point()
        keep = (pts.x >= core[0]) & (pts.x < core[2]) & (pts.y >= core[1]) & (pts.y < core[3])
        result[keep].to_file(OUT, layer="joined", driver="GPKG",
                             mode="w" if (i == 0 and j == 0) else "a")

The overlap-then-clip-to-core pattern is what makes tiled processing correct: features near a boundary are seen with their neighbours, but each one is emitted exactly once.

Push filters down to the driver

# attribute filter evaluated by GDAL — excluded rows never reach Python
gdf = gpd.read_file(SRC, engine="pyogrio",
                    where="type = 'residential' AND height > 10")

# bounding-box filter using the layer's spatial index
gdf = gpd.read_file(SRC, engine="pyogrio", bbox=(320000, 670000, 340000, 690000))

# a full SQL query, for a GeoPackage
gdf = gpd.read_file(SRC, engine="pyogrio",
                    sql="SELECT id, type, geom FROM buildings WHERE height > 30",
                    sql_dialect="SQLITE")

A where= clause that removes 90% of the rows is worth more than any amount of chunk tuning, because the work never happens at all.

Write once, append after that

from pathlib import Path

OUT = Path("data/out/buildings_area.gpkg")
OUT.unlink(missing_ok=True)        # start clean; appending to a stale file mixes runs

for i, (start, total, gdf) in enumerate(iter_row_chunks(SRC, chunk=CHUNK)):
    processed = transform(gdf)
    processed.to_file(OUT, layer="buildings", driver="GPKG", mode="w" if i == 0 else "a")

Appending requires a matching schema in every chunk, so normalise columns and dtypes inside transform() rather than per chunk. For pipelines that stay in Python, partitioned Parquet is often better than appending:

processed.to_parquet(f"data/out/buildings/part-{i:05d}.parquet", compression="zstd")
# later: gpd.read_parquet("data/out/buildings/")   reads the whole folder as one dataset

Two passes for whole-dataset answers

Some questions cannot be answered chunk by chunk. Accumulate in the first pass, apply in the second.

import numpy as np

# pass 1: accumulate the statistics you need
n, total_area, max_area = 0, 0.0, 0.0
for _, _, gdf in iter_row_chunks(SRC, columns=["geometry"]):
    areas = gdf.to_crs(3857).area
    n += len(areas); total_area += areas.sum(); max_area = max(max_area, areas.max())
mean_area = total_area / n
print(f"{n:,} features, mean area {mean_area:,.1f} m², max {max_area:,.1f} m²")

# pass 2: use the global value per chunk
for i, (_, _, gdf) in enumerate(iter_row_chunks(SRC)):
    gdf["area_ratio"] = gdf.to_crs(3857).area / mean_area
    gdf.to_file(OUT, layer="buildings", driver="GPKG", mode="w" if i == 0 else "a")

Reading the file twice is almost always cheaper than not being able to read it at all.

Code examples

Example 1: a resumable chunked processor

#!/usr/bin/env python3
from pathlib import Path
import json, time
import geopandas as gpd
import pyogrio

SRC = "data/raw/buildings.gpkg"
OUT = Path("data/out/buildings")            # a folder of Parquet parts
STATE = Path("logs/chunk_state.json")
CHUNK = 100_000

def transform(gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    metric = gdf.to_crs(gdf.estimate_utm_crs())
    gdf = gdf.copy()
    gdf["area_m2"] = metric.area.values
    gdf["perimeter_m"] = metric.length.values
    return gdf[gdf["area_m2"] > 1]

def main() -> int:
    OUT.mkdir(parents=True, exist_ok=True)
    STATE.parent.mkdir(parents=True, exist_ok=True)
    done = set(json.loads(STATE.read_text())["done"]) if STATE.exists() else set()

    total = pyogrio.read_info(SRC)["features"]
    starts = list(range(0, total, CHUNK))
    print(f"{total:,} features in {len(starts)} chunks; {len(done)} already done")

    t0 = time.perf_counter()
    for n, start in enumerate(starts, start=1):
        if start in done:
            continue
        gdf = gpd.read_file(SRC, engine="pyogrio", skip_features=start, max_features=CHUNK,
                            columns=["id", "type", "geometry"])
        if gdf.empty:
            break
        out = transform(gdf)
        tmp = OUT / f"part-{start:09d}.parquet.tmp"
        out.to_parquet(tmp, compression="zstd")
        tmp.replace(tmp.with_suffix(""))          # atomic: only complete parts are visible

        done.add(start)
        STATE.write_text(json.dumps({"done": sorted(done), "total": total}))
        elapsed = time.perf_counter() - t0
        rate = n / elapsed
        print(f"[{n}/{len(starts)}] {start + len(gdf):,}/{total:,} "
              f"eta {(len(starts) - n)/rate/60:.1f} min", flush=True)
        del gdf, out

    print("all chunks complete")
    return 0

if __name__ == "__main__":
    raise SystemExit(main())

Recording completed chunk offsets makes the job restartable after a crash, an eviction, or a coffee-related power cut.

Example 2: parallel chunks, one output file per worker

from concurrent.futures import ProcessPoolExecutor, as_completed
from pathlib import Path
import geopandas as gpd
import pyogrio

SRC, OUT, CHUNK = "data/raw/buildings.gpkg", Path("data/out/parts"), 100_000

def do_chunk(start: int) -> str:
    gdf = gpd.read_file(SRC, engine="pyogrio", skip_features=start, max_features=CHUNK)
    dest = OUT / f"part-{start:09d}.parquet"
    transform(gdf).to_parquet(dest, compression="zstd")
    return str(dest)

if __name__ == "__main__":
    OUT.mkdir(parents=True, exist_ok=True)
    total = pyogrio.read_info(SRC)["features"]
    starts = list(range(0, total, CHUNK))

    with ProcessPoolExecutor(max_workers=4) as pool:
        futures = {pool.submit(do_chunk, s): s for s in starts}
        for i, fut in enumerate(as_completed(futures), start=1):
            print(f"[{i}/{len(starts)}] {fut.result()}", flush=True)

    merged = gpd.read_parquet(OUT)
    print(f"{len(merged):,} features across {len(starts)} parts")

Every worker reads its own window and writes its own file — no shared handle, no lock, no corruption. Reading is safe to parallelise; writing to one GeoPackage from several processes is not.

Example 3: chunk a spatial join by tile, correctly

import geopandas as gpd
import pyogrio
from shapely.geometry import box

def tiled_sjoin(src, right, out, nx=6, ny=6, overlap=500):
    minx, miny, maxx, maxy = pyogrio.read_info(src)["total_bounds"]
    dx, dy = (maxx - minx) / nx, (maxy - miny) / ny
    first = True

    for i in range(nx):
        for j in range(ny):
            core = box(minx + i*dx, miny + j*dy, minx + (i+1)*dx, miny + (j+1)*dy)
            read_bbox = core.buffer(overlap).bounds
            left = gpd.read_file(src, engine="pyogrio", bbox=read_bbox)
            if left.empty:
                continue

            joined = gpd.sjoin(left, right, how="inner", predicate="intersects")
            keep = joined.representative_point().within(core)
            part = joined[keep]
            if part.empty:
                continue

            part.to_file(out, layer="joined", driver="GPKG", mode="w" if first else "a")
            first = False
            print(f"tile {i},{j}: {len(part):,} matches", flush=True)

Example 4: read a chunk stream with Fiona, no pyogrio

import fiona
import geopandas as gpd
from itertools import islice

def fiona_chunks(path, layer=None, chunk=50_000):
    with fiona.open(path, layer=layer) as src:
        crs = src.crs
        it = iter(src)
        while True:
            batch = list(islice(it, chunk))
            if not batch:
                break
            yield gpd.GeoDataFrame.from_features(batch, crs=crs)

for i, gdf in enumerate(fiona_chunks("data/raw/buildings.gpkg")):
    transform(gdf).to_parquet(f"data/out/part-{i:05d}.parquet")

Fiona streams features, so this works on formats without efficient row-offset reads — GeoJSON, GML, CSV — where skip_features would have to scan from the start each time.

Explanation

read_file() builds one Python object holding every feature, so peak memory is proportional to the dataset. Chunking replaces "size of dataset" with "size of chunk" in that formula, which is the entire trick. Everything else is about correctness: making sure the answer computed piecewise equals the answer you would have got whole.

A large layer split into row windows and into spatial tiles, with a boundary feature highlighted. — Row windows are cheap and arbitrary; spatial tiles respect the neighbourhood a feature needs.

That is why the choice of chunking axis follows the operation. Row-independent work — reprojection, attribute calculation, buffering, validity repair, format conversion — gives identical results whatever the window boundaries, because each feature's output depends only on itself. Neighbour-dependent work does not: a dissolve across a window boundary produces two half-merged polygons, and a nearest-neighbour search inside a window misses the nearer feature just outside it.

Spatial tiling fixes that, with one subtlety. Reading a tile with an overlap gives boundary features their neighbours, but it also means the same feature is read by several tiles. Emitting only the features whose representative point falls in the core rectangle guarantees each one appears exactly once in the output. Choose the overlap to be at least as large as the biggest feature or the largest distance the operation looks across, or you reintroduce the edge error you were avoiding.

Reading is also where the real savings are. columns=, bbox=, where= and sql= are evaluated inside GDAL, so filtered rows are never materialised as Python objects at all. Halving the columns halves the memory and the parse time before any chunking logic runs. In practice a where= clause that drops most of the data is worth more than any amount of chunk-size tuning.

On the write side, appending into one GeoPackage works and requires a stable schema across chunks; writing a folder of Parquet parts is usually faster, compresses better, parallelises safely and can be read back as a single dataset. And because each part is written atomically, a crashed run leaves a set of complete parts plus a record of which offsets are done — which is what makes a chunked job resumable rather than restartable.

Edge cases or notes

skip_features is not free on every format: GeoPackage supports efficient row offsets; GeoJSON and CSV must scan, so a windowed loop over them is quadratic. Stream with Fiona instead.
The order of features is driver-dependent: Do not assume a stable order between runs unless you sort explicitly with sql=.
Appending needs a matching schema: Normalise columns and dtypes inside the transform, or mode="a" fails midway through the run.
estimate_utm_crs() per chunk can differ: For a dataset spanning several UTM zones, compute the target CRS once from the full extent and pass it in.
Don't chunk a dissolve by rows: Group by a spatial key or dissolve per tile and merge the results; a row window will split groups arbitrarily.
Spatial index rebuilds: Each bbox= read uses the source's index — good. Building an index over a huge GeoPackage once (CREATE VIRTUAL TABLE ... rtree) makes tiled reads dramatically faster.
Progress needs flushing: In a job that runs for hours, print(..., flush=True) or a logger is the difference between visible progress and a silent terminal.

Internal links

FAQ

How big should a chunk be?

Large enough to amortise the per-read overhead, small enough that a chunk plus its intermediates fits comfortably in RAM. Measure bytes per feature on a small sample and divide your memory budget by three times that number; 50k–250k features is typical.

Which operations are safe to chunk by row?

Anything where a feature's result depends only on itself: reprojection, attribute maths, buffering, simplification, validity repair, filtering and format conversion. Dissolve, spatial join, nearest-neighbour and aggregation are not.

How do I avoid duplicates when tiling with an overlap?

Read with an overlap so boundary features see their neighbours, then keep only the results whose representative point falls inside the tile's core rectangle. Each feature is then emitted exactly once.

Should I append to a GeoPackage or write Parquet parts?

Parts, if the consumer is Python: they are faster, compress better, parallelise safely and read back as one dataset. Append to a GeoPackage when the deliverable must be a single portable file.

Why is chunking a GeoJSON so slow?

skip_features cannot seek in a text format, so every window re-scans from the beginning. Stream features with Fiona instead, or convert once to GeoPackage or Parquet and chunk that.

How do I compute a dataset-wide statistic?

Two passes: accumulate counts, sums and extremes in the first, then apply the global value per chunk in the second. Reading twice is far cheaper than not being able to read at all.

Can I process chunks in parallel?

Yes for reading — each worker reads its own window. Give each worker its own output file and merge afterwards; several processes writing one GeoPackage will corrupt it.

How to Process a Very Large GeoPackage in Chunks with Python #

Problem statement #

Quick answer #

Row chunks or spatial chunks #

Step-by-step solution #

Inspect before you read #

Choose a chunk size from memory, not from habit #

Chunk by row window for row-independent work #

Chunk by space when neighbours matter #

Push filters down to the driver #

Write once, append after that #

Two passes for whole-dataset answers #

Code examples #

Example 1: a resumable chunked processor #

Example 2: parallel chunks, one output file per worker #

Example 3: chunk a spatial join by tile, correctly #

Example 4: read a chunk stream with Fiona, no pyogrio #

Explanation #

Edge cases or notes #

Internal links #

FAQ #

How big should a chunk be? #

Which operations are safe to chunk by row? #

How do I avoid duplicates when tiling with an overlap? #

Should I append to a GeoPackage or write Parquet parts? #

Why is chunking a GeoJSON so slow? #

How do I compute a dataset-wide statistic? #

Can I process chunks in parallel? #