How to Process a Very Large GeoPackage in Chunks with Python

On this page

Problem statement

One file, three quarters of a million buildings, and a machine — or a container, or a CI runner — with less memory than reading it takes. The obvious line does not work:

%%writefile read_all.py
import geopandas as gpd

gdf = gpd.read_file("data/buildings.gpkg")
gdf["area_m2"] = gdf.area
print(len(gdf), "buildings")

import subprocess, sys

# a 500 MB memory ceiling, the way a container or a small cloud instance imposes one
run = subprocess.run(["systemd-run", "--user", "--scope", "-q", "-p", "MemoryMax=500M", "-p", "MemorySwapMax=0",
                      sys.executable, "read_all.py"], capture_output=True, text=True)
print("exit status", run.returncode, "(-9 is SIGKILL: the kernel killed it)" if run.returncode == -9 else "")

exit status -9 (-9 is SIGKILL: the kernel killed it)

There is no bad file to skip and no loop to fix — the layer simply does not fit in the memory available. 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

The file on this page is every building in Geofabrik's OpenStreetMap extract of the West Midlands dated 1 January 2026 — 745,849 polygons, 183 MB as a GeoPackage in British National Grid.

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="w" for the first write, "a" after)
delete the chunk and report progress
for neighbour-dependent work, chunk by space with bbox= instead of by row

%%writefile chunked.py
from pathlib import Path
import geopandas as gpd
import pyogrio

SRC = "data/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, skip_features=start, max_features=CHUNK, columns=["osm_id", "type"])
    gdf["area_m2"] = gdf.area                      # British National Grid is already in metres
    gdf.to_file(OUT, layer="buildings", driver="GPKG", mode="w" if start == 0 else "a")
    print(f"[{start + len(gdf):>8,}/{total:,}] {(start + len(gdf)) / total:5.1%}", flush=True)
    del gdf

run = subprocess.run(["systemd-run", "--user", "--scope", "-q", "-p", "MemoryMax=500M", "-p", "MemorySwapMax=0",
                      sys.executable, "chunked.py"], capture_output=True, text=True)
print("exit status", run.returncode)
print(run.stdout)

exit status 0
745,849 features, CRS EPSG:27700
[ 100,000/745,849] 13.4%
[ 200,000/745,849] 26.8%
[ 300,000/745,849] 40.2%
[ 400,000/745,849] 53.6%
[ 500,000/745,849] 67.0%
[ 600,000/745,849] 80.4%
[ 700,000/745,849] 93.9%
[ 745,849/745,849] 100.0%

The same job under the same ceiling finishes. read_info() is a header read, so the feature count costs nothing. columns= narrows the read at the driver, so unwanted columns never become Python objects — though with only two short attribute columns here, geometry is most of the memory.

Row chunks or spatial chunks

Four ways to cut a dataset, and which operations each one supports.

Measure, then chunk by rows or by space

Six steps, and the last one is what keeps memory flat.

Inspect before you read

import time
import geopandas as gpd
import pyogrio

SRC = "data/buildings.gpkg"
start = time.perf_counter()
info = pyogrio.read_info(SRC)
print(f"read_info   : {(time.perf_counter() - start) * 1000:.1f} ms")
print(f"layers      : {pyogrio.list_layers('data/buildings.gpkg')[:, 0].tolist()}")
print(f"features    : {info['features']:,}")
print(f"geometry    : {info['geometry_type']}")
print(f"crs         : {info['crs']}")
print(f"fields      : {list(info['fields'])}")
print(f"bounds      : {[round(v) for v in info['total_bounds']]}")

read_info   : 2.2 ms
layers      : ['buildings']
features    : 745,849
geometry    : MultiPolygon
crs         : EPSG:27700
fields      : ['osm_id', 'type']
bounds      : [386081, 271985, 439207, 307281]

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

Measure a sample — but measure the process, not the DataFrame, and in a fresh process so memory freed by earlier work does not hide the growth. memory_usage(deep=True) counts the geometry column as one pointer per row and misses the geometries themselves:

%%writefile measure_sample.py
import psutil
import geopandas as gpd

process = psutil.Process()
gpd.read_file("data/buildings.gpkg", max_features=10)      # load the libraries first
before = process.memory_info().rss
sample = gpd.read_file("data/buildings.gpkg", max_features=50_000, columns=["osm_id", "type"])
after = process.memory_info().rss
print(sample.memory_usage(deep=True).sum() / len(sample), (after - before) / len(sample))

naive, measured = map(float, subprocess.run([sys.executable, "measure_sample.py"],
                                             capture_output=True, text=True).stdout.split())
print(f"memory_usage(deep=True): {naive:,.0f} bytes per feature")
print(f"process memory growth:   {measured:,.0f} bytes per feature")

BUDGET_GB = 0.25
chunk = int(BUDGET_GB * 1e9 / (measured * 3))     # x3 headroom for intermediates
print(f"chunk size for a {BUDGET_GB} GB budget: about {chunk:,} features")

memory_usage(deep=True): 39 bytes per feature
process memory growth:   1,070 bytes per feature
chunk size for a 0.25 GB budget: about 77,908 features

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, skip_features=start, max_features=chunk,
                            columns=columns, where=where)
        if gdf.empty:
            break
        yield start, total, gdf

for start in (0, 300_000, 700_000):
    t = time.perf_counter()
    gpd.read_file(SRC, skip_features=start, max_features=50_000)
    print(f"window at row {start:>7,}: {time.perf_counter() - t:.2f} s")

window at row       0: 0.10 s
window at row 300,000: 0.09 s
window at row 700,000: 0.11 s

In a GeoPackage a window near the end costs the same as one at the start, so a row-window loop is linear. 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 spatial join spanning a window boundary still works row by row, but anything that looks at neighbours — dissolve, nearest, a join against features read in the same tile — gives the wrong answer at the edges. Tile the extent instead, and overlap the tiles. Here every building is joined to the ONS ward it touches, whole and in 36 overlapping tiles:

import pandas as pd
from shapely.geometry import box

wards = gpd.read_file("data/wards.gpkg")[["WD24CD", "geometry"]]
whole = gpd.sjoin(gpd.read_file(SRC), wards, predicate="intersects")

minx, miny, maxx, maxy = pyogrio.read_info(SRC)["total_bounds"]
NX = NY = 6
OVERLAP = 500          # CRS units — at least the largest feature you expect
dx, dy = (maxx - minx) / NX, (maxy - miny) / NY

def tiled(clip_to_core):
    parts = []
    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)
            tile = gpd.read_file(SRC, bbox=core.buffer(OVERLAP).bounds)
            if tile.empty:
                continue
            result = gpd.sjoin(tile, wards, predicate="intersects")
            if clip_to_core:
                # keep only results whose representative point is in the core tile
                result = result[result.representative_point().within(core)]
            parts.append(result)
    return pd.concat(parts)

pairs = lambda frame: frame["osm_id"].astype(str) + "|" + frame["WD24CD"]
for clip in (False, True):
    result = tiled(clip)
    print(f"clip to core {clip!s:5}: {len(result):,} rows, {pairs(result).duplicated().sum():,} duplicates; "
          f"whole-file join {len(whole):,} rows, same pairs: {set(pairs(result)) == set(pairs(whole))}")

clip to core False: 977,633 rows, 229,918 duplicates; whole-file join 747,715 rows, same pairs: True
clip to core True : 747,715 rows, 0 duplicates; whole-file join 747,715 rows, same pairs: True

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. Without the clip, the overlapping reads emit nearly a quarter of a million joined rows twice.

Measured above. The clipped tiles produce exactly the same building-ward pairs as the whole-file join.

Push filters down to the driver

t = time.perf_counter()
houses = gpd.read_file(SRC, where="type = 'house'")
print(f"where=: {len(houses):,} houses in {time.perf_counter() - t:.2f} s")

t = time.perf_counter()
centre = gpd.read_file(SRC, bbox=(406000, 286000, 408000, 288000))       # 2 km square on Birmingham
print(f"bbox=:  {len(centre):,} buildings in {time.perf_counter() - t:.2f} s")

t = time.perf_counter()
large = gpd.read_file(SRC, sql="SELECT osm_id, type, geom FROM buildings WHERE ST_Area(geom) > 5000",
                      sql_dialect="SQLITE")
print(f"sql=:   {len(large):,} buildings over 5,000 m² in {time.perf_counter() - t:.2f} s")

where=: 24,153 houses in 0.12 s
bbox=:  2,135 buildings in 0.01 s
sql=:   1,250 buildings over 5,000 m² in 0.18 s

A where= clause that removes 97% of the rows is worth more than any amount of chunk tuning, because the work never happens at all. Two things to know: the column is called type, and gdf.type is GeoPandas' geometry-type property, so use gdf["type"]; and bbox= is fast because GDAL wrote an R-tree into the GeoPackage when it created the layer.

Write once, append after that

from pathlib import Path

def transform(gdf):
    gdf = gdf.copy()
    gdf["area_m2"] = gdf.area
    gdf["perimeter_m"] = gdf.length
    return gdf[gdf["area_m2"] > 1]

OUT = Path("data/out/buildings_metrics.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=100_000)):
    transform(gdf).to_file(OUT, layer="buildings", driver="GPKG", mode="w" if i == 0 else "a")
print(pyogrio.read_info(OUT)["features"], "features written")

745741 features written

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:

parts = Path("data/out/buildings_parts")
parts.mkdir(parents=True, exist_ok=True)
for i, (start, total, gdf) in enumerate(iter_row_chunks(SRC, chunk=100_000)):
    transform(gdf).to_parquet(parts / f"part-{i:05d}.parquet", compression="zstd")
print(len(gpd.read_parquet(parts)), "features read back from", len(list(parts.glob("*.parquet"))), "parts")

745741 features read back from 8 parts

Two passes for whole-dataset answers

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

# 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=[]):
    areas = gdf.area                                   # metres: the layer is in EPSG:27700
    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:,.0f} m²")

# the same mean in Web Mercator, which is in metres but not equal-area
mercator = sum(g.to_crs(3857).area.sum() for _, _, g in iter_row_chunks(SRC, columns=[])) / n
print(f"in EPSG:3857 the mean would be {mercator:,.1f} m² ({mercator / mean_area:.2f}x)")

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

745,849 features, mean area 148.1 m², max 210,149 m²
in EPSG:3857 the mean would be 399.2 m² (2.70x)

Reading the file twice is almost always cheaper than not being able to read it at all. Measure areas in an equal-area or local projected CRS: Web Mercator inflates them by about 1/cos²(latitude), which in the West Midlands is 2.7.

Chunking the West Midlands' buildings

The data used on this page

Every building polygon in Geofabrik's OpenStreetMap extract of the West Midlands dated 1 January 2026 (© OpenStreetMap contributors, ODbL), with its osm_id and type, reprojected to British National Grid and written to data/buildings.gpkg; the first 200,000 also as data/buildings.geojson in WGS 84 for Example 4. The 301 ONS wards (December 2024 generalised boundaries) that intersect the buildings' extent are in data/wards.gpkg.

import io
import os
from pathlib import Path

import geopandas as gpd
import requests

Path("data").mkdir(exist_ok=True)
if not Path("data/buildings.gpkg").exists():
    osm = Path("data/west-midlands-260101-free.shp.zip")
    if not osm.exists():
        osm.write_bytes(requests.get(
            "https://download.geofabrik.de/europe/united-kingdom/england/west-midlands-260101-free.shp.zip",
            timeout=600).content)
    buildings = gpd.read_file(f"/vsizip/{osm}/gis_osm_buildings_a_free_1.shp", columns=["osm_id", "type"])
    buildings.to_crs(27700).to_file("data/buildings.gpkg", layer="buildings")
    buildings.iloc[:200_000].to_file("data/buildings.geojson", driver="GeoJSON")
    del buildings
if not Path("data/wards.gpkg").exists():
    import pyogrio
    minx, miny, maxx, maxy = pyogrio.read_info("data/buildings.gpkg")["total_bounds"]
    reply = requests.get(
        "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
        "Wards_December_2024_Boundaries_UK_BGC/FeatureServer/0/query",
        params={"geometry": f"{minx},{miny},{maxx},{maxy}", "geometryType": "esriGeometryEnvelope",
                "inSR": 27700, "spatialRel": "esriSpatialRelIntersects", "outFields": "WD24CD,WD24NM,LAD24NM",
                "outSR": 27700, "f": "geojson"}, timeout=300)
    reply.raise_for_status()
    gpd.read_file(io.BytesIO(reply.content)).to_file("data/wards.gpkg")

print(f"buildings.gpkg {os.path.getsize('data/buildings.gpkg') / 2**20:.0f} MB, "
      f"buildings.geojson {os.path.getsize('data/buildings.geojson') / 2**20:.0f} MB, "
      f"{len(gpd.read_file('data/wards.gpkg'))} wards")

buildings.gpkg 183 MB, buildings.geojson 62 MB, 301 wards

Example 1: a resumable chunked processor

%%writefile resumable.py
from pathlib import Path
import json, time
import geopandas as gpd
import pyogrio

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

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

def main(stop_after=None) -> 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")

    for n, start in enumerate(starts, start=1):
        if start in done:
            continue
        if stop_after is not None and len(done) >= stop_after:
            raise RuntimeError(f"simulated crash after {len(done)} chunks")
        gdf = gpd.read_file(SRC, skip_features=start, max_features=CHUNK, columns=["osm_id", "type"])
        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}))
        print(f"[{n}/{len(starts)}] {start + len(gdf):,}/{total:,}", flush=True)
        del gdf, out

    print("all chunks complete")
    return 0

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

import shutil
shutil.rmtree("data/out/resumable", ignore_errors=True)
Path("logs/chunk_state.json").unlink(missing_ok=True)

import resumable
try:
    resumable.main(stop_after=3)             # the first run dies part-way
except RuntimeError as err:
    print("!", err)
resumable.main()                            # the rerun picks up where it stopped
print(len(gpd.read_parquet("data/out/resumable")), "features in", len(list(Path("data/out/resumable").glob("*.parquet"))), "parts")

745,849 features in 8 chunks; 0 already done
[1/8] 100,000/745,849
[2/8] 200,000/745,849
[3/8] 300,000/745,849
! simulated crash after 3 chunks
745,849 features in 8 chunks; 3 already done
[4/8] 400,000/745,849
[5/8] 500,000/745,849
[6/8] 600,000/745,849
[7/8] 700,000/745,849
[8/8] 745,849/745,849
all chunks complete
745741 features in 8 parts

Recording completed chunk offsets makes the job restartable after a crash, an eviction, or a coffee-related power cut: the second run skips the three finished chunks and does the other five.

Example 2: parallel chunks, one output file per worker

A process pool needs its worker function in an importable module:

%%writefile chunk_worker.py
from pathlib import Path
import geopandas as gpd
from resumable import transform

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

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

def append_to_one_file(start: int) -> str:
    gdf = gpd.read_file(SRC, skip_features=start, max_features=CHUNK)
    transform(gdf).to_file("data/out/shared.gpkg", layer="buildings", mode="a")
    return "ok"

import os
from concurrent.futures import ProcessPoolExecutor
import chunk_worker

chunk_worker.OUT.mkdir(parents=True, exist_ok=True)
starts = list(range(0, pyogrio.read_info(SRC)["features"], chunk_worker.CHUNK))

t = time.perf_counter()
for s in starts:
    chunk_worker.do_chunk(s)
print(f"sequential: {time.perf_counter() - t:.1f} s")

t = time.perf_counter()
with ProcessPoolExecutor(max_workers=4) as pool:
    list(pool.map(chunk_worker.do_chunk, starts))
print(f"4 workers:  {time.perf_counter() - t:.1f} s  (load average {os.getloadavg()[0]:.1f} on {os.cpu_count()} CPUs)")
print(f"{len(gpd.read_parquet(chunk_worker.OUT)):,} features across {len(starts)} parts")

sequential: 3.4 s
4 workers:  1.7 s  (load average 3.8 on 24 CPUs)
745,741 features across 8 parts

Every worker reads its own window and writes its own file — no shared handle, no lock. Reading is safe to parallelise. Writing to one GeoPackage from several processes is not, and the way it fails is the dangerous kind:

import sqlite3

Path("data/out/shared.gpkg").unlink(missing_ok=True)
resumable.transform(gpd.read_file(SRC, max_features=1)).iloc[:0].to_file("data/out/shared.gpkg", layer="buildings")
with ProcessPoolExecutor(max_workers=4) as pool:
    results = list(pool.map(chunk_worker.append_to_one_file, starts))

rows = sqlite3.connect("data/out/shared.gpkg").execute("SELECT count(*) FROM buildings").fetchone()[0]
print("every append reported:", set(results))
print(f"rows in the table: {rows:,}; feature count GDAL reports: {pyogrio.read_info('data/out/shared.gpkg')['features']:,}")

every append reported: {'ok'}
rows in the table: 745,741; feature count GDAL reports: 545,741

No error, and the rows are all there — but the count GeoPackage keeps in its metadata was updated by four writers at once and is wrong. Anything that plans from that count, including the row-window loop at the top of this page, then silently reads only part of the file.

Example 3: chunk a spatial join by tile, correctly

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, written = True, 0

    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)
            left = gpd.read_file(src, bbox=core.buffer(overlap).bounds)
            if left.empty:
                continue

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

            part.to_file(out, layer="joined", driver="GPKG", mode="w" if first else "a")
            first, written = False, written + len(part)
    return written

Path("data/out/joined.gpkg").unlink(missing_ok=True)
print(tiled_sjoin(SRC, wards, "data/out/joined.gpkg"), "rows written; the whole-file join has", len(whole))

747715 rows written; the whole-file join has 747715

Example 4: stream a format that cannot seek

skip_features cannot seek in a text format. GDAL has to parse a GeoJSON before it can hand over row 150,000, so every window pays for the parse again. Stream it once instead:

import shapely
from itertools import islice
import fiona

GEOJSON = "data/buildings.geojson"
total = pyogrio.read_info(GEOJSON)["features"]

t = time.perf_counter()
for start in range(0, total, 50_000):
    gpd.read_file(GEOJSON, skip_features=start, max_features=50_000)
print(f"row windows:        {time.perf_counter() - t:.1f} s")

t = time.perf_counter()
with pyogrio.open_arrow(GEOJSON, batch_size=50_000, use_pyarrow=True) as (meta, reader):
    geometry_column = meta["geometry_name"] or "wkb_geometry"
    for batch in reader:
        table = batch.to_pandas()
        gdf = gpd.GeoDataFrame(table.drop(columns=geometry_column),
                               geometry=shapely.from_wkb(table[geometry_column]), crs=meta["crs"])
print(f"pyogrio open_arrow: {time.perf_counter() - t:.1f} s")

t = time.perf_counter()
with fiona.open(GEOJSON) as src:
    features = iter(src)
    while batch := list(islice(features, 50_000)):
        gdf = gpd.GeoDataFrame.from_features(batch, crs=src.crs)
print(f"Fiona, per feature: {time.perf_counter() - t:.1f} s")

row windows:        3.6 s
pyogrio open_arrow: 1.5 s
Fiona, per feature: 7.3 s

open_arrow reads the file once and hands it over in batches, so memory stays at one batch and the parse is paid once. Fiona streams too, but it builds a Python dictionary per feature, which here costs more than re-parsing the whole file for every window.

Measured in Example 4, one run. Windows re-parse the file each time; Fiona streams but builds a Python dictionary per feature.

Swapping dataset size for chunk size

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.

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. 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 from one process 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.

Formats that cannot seek, and other chunking costs

skip_features is not free on every format: GeoPackage seeks — a window at row 700,000 took as long as one at row 0 — but GeoJSON and CSV must be parsed up to the window, so a windowed loop over them pays the parse once per window. Stream with pyogrio.open_arrow 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 — or use a national grid, as here.
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.
Check the spatial index exists: GDAL writes an R-tree (rtree_<layer>_<geometry>) into every GeoPackage layer it creates, which is what makes bbox= fast. A GeoPackage from another tool may not have one; ogrinfo -sql "SELECT CreateSpatialIndex('layer', 'geom')" adds it.
One writer per GeoPackage: Concurrent appends reported success here and left a wrong feature count in the file's metadata.
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 the process's memory growth while reading a sample — not memory_usage(deep=True), which misses the geometries — and divide your memory budget by three times the bytes per feature; 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 pays for parsing the file again. Stream it once with pyogrio.open_arrow, or convert once to GeoPackage or Parquet and chunk that. Fiona streams too, but its per-feature dictionaries made it the slowest of the three here.

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. Four processes appending to one GeoPackage all reported success here and left the file's stored feature count wrong, which misleads anything that plans from it.