Reprojecting Large Datasets Without Memory Errors

Calling .to_crs() on a multi-gigabyte layer loads the whole thing into RAM and frequently dies with a MemoryError. This guide reprojects datasets larger than memory by streaming them in Arrow record batches, never holding more than one batch at a time. It is for anyone reprojecting national or continental vector layers on a normal machine, and it sits under Coordinate Reference System Transformations in Geospatial Data Ingestion & Processing Workflows. It builds directly on the GeoDataFrame reprojection API and the axis-order rules covered in Coordinate Systems with PyProj.

Why This Approach / What Goes Wrong

The naive one-liner, gpd.read_file(...).to_crs(...), is a three-way memory spike. The full source is materialised into a single GeoDataFrame, .to_crs() allocates a second reprojected copy of every geometry, and the writer buffers its output on top — three copies of a huge dataset resident at once. On a 40-million-feature parcel layer that is tens of gigabytes of live objects, so the process is killed long before it finishes.

The fix is to stream: read a bounded batch of features, reproject just that batch, append it to the output, release it, and repeat. Peak memory is then governed by the batch size, not the dataset size, so a 60 GB layer reprojects in the same footprint as a 600 MB one. Modern GDAL exposes an Arrow-native batch reader through pyogrio.open_arrow, and OGR writers such as FlatGeobuf and GeoPackage accept incremental appends, which makes the loop clean and fast.

It is worth doing the arithmetic once, because the number that matters is not the file size on disk. A Shapely 2 geometry column holds a NumPy array of pointers to GEOS objects; each polygon costs 16 bytes per vertex for the coordinate sequence plus roughly 150–250 bytes of GEOS and Python object overhead. A cadastral parcel averaging 40 vertices is therefore around 800–900 bytes live, against maybe 250 bytes in a compressed FlatGeobuf. Forty-two million of them is close to 37 GB before a single attribute column is counted — and object-dtype string columns in pandas add another 50–60 bytes per value on top. That is the first copy. .to_crs() builds a complete second geometry array before releasing the first, and the OGR writer buffers a third. A "12 GB file" routinely needs 90 GB of RAM to reproject in one shot, which is why the failure feels so disproportionate to the input.

Two correctness traps recur when people first stream a reprojection:

• Never assume the source has a CRS. Many large public datasets — national cadastres, bulk building-footprint dumps — ship without an embedded CRS. Reprojecting a layer whose .crs is None raises, or worse, silently mislabels coordinates. Read the source CRS once from the file metadata and set_crs it onto any batch that arrives without one.
• Let GeoPandas own the axis order. A hand-built pyproj.Transformer without always_xy=True swaps latitude and longitude and flips your output across the globe — the single most common transformation bug, diagnosed in Fixing PyProj CRS Transformation Errors. GeoDataFrame.to_crs() already handles axis order correctly, so call it per batch rather than reimplementing the transform.

Prerequisites

• geopandas>=1.0 — GeoDataFrame.from_arrow and native GeoArrow I/O land in 1.0
• pyogrio>=0.8 — the Arrow batch reader (open_arrow) and incremental write_dataframe
• pyarrow>=15 — the Arrow tables the batches flow through

conda install -c conda-forge "geopandas=1.0.*" "pyogrio=0.8.*" "pyarrow=15.*"

Install from conda-forge, not pip, so GDAL, PROJ, and the Python bindings stay ABI-compatible — the same environment discipline that avoids the PROJ-path failures seen across Coordinate Systems with PyProj.

Step-by-Step Implementation

1. Inspect the source CRS and feature count without loading any geometry. read_info reads only the header, so it is instant even on a 60 GB file.

import pyogrio

info = pyogrio.read_info("national_parcels.fgb")
print("CRS:", info["crs"], "| features:", info["features"])
# CRS: EPSG:4326 | features: 41872330

source_crs = info["crs"]
if source_crs is None:
    source_crs = "EPSG:4326"   # assert the known source CRS; do NOT guess silently

2. Pick a metric target CRS, not Web Mercator. Reprojection for any distance or area work needs an equal-area or transverse-Mercator CRS in metres. For a single UTM zone that is ETRS89 / UTM 32N (EPSG:25832); if the layer spans several zones, choose the UTM zone programmatically or use a national equal-area CRS. Avoid EPSG:3857 (Web Mercator) — its scale distortion makes every downstream measurement wrong.

TARGET_EPSG = 25832   # ETRS89 / UTM 32N — metric, suited to area and distance work
BATCH = 250_000       # features per batch; the single knob that bounds peak RAM

3. Stream Arrow batches, reproject each, and append to a single GeoPackage. open_arrow yields fixed-size record batches; each is converted to a GeoDataFrame, reprojected, and written incrementally. Only one batch is ever resident.

import geopandas as gpd
import pyogrio
import pyarrow as pa

SOURCE = "national_parcels.fgb"
OUTPUT = "parcels_utm.gpkg"

with pyogrio.open_arrow(SOURCE, use_pyarrow=True, batch_size=BATCH) as (meta, reader):
    first = True
    for record_batch in reader:
        batch = gpd.GeoDataFrame.from_arrow(pa.Table.from_batches([record_batch]))
        if batch.crs is None:                    # source shipped without a CRS
            batch = batch.set_crs(source_crs)    # set_crs asserts; it does NOT reproject
        reprojected = batch.to_crs(epsg=TARGET_EPSG)
        pyogrio.write_dataframe(
            reprojected, OUTPUT, layer="parcels",
            append=not first,                    # create on the first batch, append after
        )
        first = False
        del batch, reprojected                   # release before pulling the next batch

4. Prefer a partitioned GeoParquet folder for columnar consumers. If the output feeds DuckDB or Dask, write one GeoParquet file per batch and treat the folder as a single dataset — no format needs to support in-place append, and the writes parallelise trivially.

import os
import geopandas as gpd
import pyogrio
import pyarrow as pa

os.makedirs("parcels_utm_parts", exist_ok=True)
with pyogrio.open_arrow(SOURCE, use_pyarrow=True, batch_size=BATCH) as (meta, reader):
    for i, record_batch in enumerate(reader):
        batch = gpd.GeoDataFrame.from_arrow(pa.Table.from_batches([record_batch]))
        if batch.crs is None:
            batch = batch.set_crs(source_crs)
        batch.to_crs(epsg=TARGET_EPSG).to_parquet(
            f"parcels_utm_parts/part_{i:05d}.parquet"
        )

5. Make the run resumable. A reprojection that takes four hours will eventually be interrupted — a full disk, a preempted spot instance, an OOM killer taking the wrong process. With the part-per-batch layout, resumption is a set lookup: record each completed batch in an append-only manifest and skip the ones already on disk.

import json
from pathlib import Path

MANIFEST = Path("parcels_utm_parts/_manifest.jsonl")
done = set()
if MANIFEST.exists():
    done = {json.loads(line)["batch"] for line in MANIFEST.open()}

with pyogrio.open_arrow(SOURCE, use_pyarrow=True, batch_size=BATCH) as (meta, reader):
    for i, record_batch in enumerate(reader):
        if i in done:
            continue                                    # already written
        batch = gpd.GeoDataFrame.from_arrow(pa.Table.from_batches([record_batch]))
        if batch.crs is None:
            batch = batch.set_crs(source_crs)
        batch.to_crs(epsg=TARGET_EPSG).to_parquet(
            f"parcels_utm_parts/part_{i:05d}.parquet"
        )
        with MANIFEST.open("a") as fh:                   # durable after each batch
            fh.write(json.dumps({"batch": i, "rows": len(batch)}) + "\n")

A resumed run still reads the skipped batches — the reader is a forward-only scan and cannot seek to batch 900 — but it does not transform or write them. On a 60 GB source the re-read costs minutes against hours of transform and write, so the trade is worth taking. Do not attempt the same trick against an appended GeoPackage: without a manifest there is no way to tell how many features from the interrupted batch actually committed, and an unconditional retry duplicates them.

6. Parallelise by row range, not by bounding box. When one core is the bottleneck and the source format supports random access, hand each worker a disjoint slice of features. skip_features/max_features partitions exactly; a spatial bbox filter does not, because every feature straddling a tile edge is returned to both neighbours and lands in the output twice.

from concurrent.futures import ProcessPoolExecutor
import geopandas as gpd
import pyogrio

SOURCE = "national_parcels.fgb"
TARGET_EPSG = 25832
TOTAL = pyogrio.read_info(SOURCE)["features"]
SLICE = 1_000_000

def reproject_slice(job):
    index, start = job
    part = pyogrio.read_dataframe(
        SOURCE, skip_features=start, max_features=SLICE,
    )
    part.to_crs(epsg=TARGET_EPSG).to_parquet(
        f"parcels_utm_parts/slice_{index:04d}.parquet"
    )
    return index, len(part)

jobs = list(enumerate(range(0, TOTAL, SLICE)))
with ProcessPoolExecutor(max_workers=8) as pool:
    for index, n in pool.map(reproject_slice, jobs):
        print(f"slice {index}: {n} features")

Each worker builds its own PROJ context after the fork, which is what you want — a Transformer created in the parent and inherited across a fork is not safe, and the symptom is a segfault or a worker that hangs rather than a clean exception. Size SLICE so one slice still fits in a single worker's share of RAM: eight workers on a 32 GB machine have roughly 4 GB each, not 32.

Throughput past that point is bounded by three different things in turn. PROJ transforms on the order of a million coordinates per second per core for a projection-only conversion, and several times slower when a grid-based datum shift is in the pipeline. The Arrow reader will typically saturate a spinning disk or a network volume long before it saturates the CPU, so on cloud storage the read is the ceiling and adding workers does nothing. And a single appended GeoPackage serialises every write behind one writer — the moment that is the limit, switch to the part-per-batch folder, which is the same reason partitioning strategies for Dask-GeoPandas favour many mid-sized files over one large one.

Verification

Confirm the output carries the target CRS and that no batch was silently dropped — again without reloading everything at once. read_info on the result reads only its header.

import pyogrio

out = pyogrio.read_info("parcels_utm.gpkg", layer="parcels")
print("Output CRS:", out["crs"], "| features:", out["features"])
# Output CRS: EPSG:25832 | features: 41872330

assert out["features"] == 41_872_330, "Feature count changed — a batch was dropped"
assert "25832" in str(out["crs"]), "Output is not in the target CRS"

# Peak RSS should plateau near one batch's worth of geometry, not climb with the
# feature count. Watch it during the run: /usr/bin/time -v python reproject.py
# and confirm "Maximum resident set size" stays flat as batches stream through.

A matching feature count proves nothing about the coordinates, and a streamed job is exactly where a silently wrong transform hides — every batch is wrong in the same way, so nothing looks anomalous. Two cheap checks on the head of the output cover the realistic failures. Coordinates that fell outside the projection come back as inf rather than raising, so a finiteness assertion catches the stray out-of-zone record. And a round-trip back to the source CRS proves the write did not quietly lose precision — a Shapefile hop, or a format that stores single-precision floats, shows up here as a residual of metres instead of nanometres.

import numpy as np
import pyogrio

head = pyogrio.read_dataframe("parcels_utm.gpkg", layer="parcels", max_features=50_000)

# 1. Plausibility: UTM eastings are 100k–900k by construction, and nothing is inf
minx, miny, maxx, maxy = head.total_bounds
assert np.isfinite([minx, miny, maxx, maxy]).all(), "inf coordinates in the output"
assert 100_000 < minx and maxx < 900_000, f"easting outside UTM range: {minx}–{maxx}"

# 2. Round-trip: reproject the sample back and compare against the source features
src_head = pyogrio.read_dataframe("national_parcels.fgb", max_features=50_000)
residual = head.to_crs(src_head.crs).geometry.distance(src_head.geometry, align=False)
assert residual.max() < 1e-9, f"round-trip residual {residual.max()} deg — wrong CRS"
print(f"max round-trip residual: {residual.max():.2e} degrees")

A residual in the tenth decimal place is the pipeline being reversible, as it should be. Note what this check cannot tell you: if step 1 asserted the wrong source CRS, the round trip still closes perfectly, because the same wrong assumption is applied in both directions. Detecting that needs an outside reference — overlay a few hundred output features against a trusted layer in the target CRS, or confirm the output's total_bounds sit inside the known extent of the region the data covers. A national parcel layer whose reprojected bounds are 200 km from the country it describes had its source CRS declared wrong, and no self-consistency check will ever say so.

Edge Cases & Debugging

• MemoryError still fires. Lower BATCH — 250k vertex-dense polygons can still be large. Halving it halves peak RAM.
• Source CRS is None. The file has no embedded CRS; set_crs each batch to the known source EPSG before to_crs, never after.
• Axis-flipped output (points land in the ocean). A hand-rolled pyproj.Transformer without always_xy=True. Delete it and call batch.to_crs(...), which orders axes correctly — see Fixing PyProj CRS Transformation Errors.
• Slow throughput. Confirm use_pyarrow=True on open_arrow; the Arrow path is far faster than the legacy per-feature reader.
• Mixed geometry types across batches. GeoPackage tolerates a generic geometry column; FlatGeobuf and Shapefile want one type per layer, so split by geometry type or promote to multi-geometry before writing.
• Need to parallelise the whole reprojection. When one core is the bottleneck, move to Dask-GeoPandas, which partitions and reprojects across workers.
• read_info reports features: None. The source is a streaming format — plain GeoJSON, a compressed CSV of WKT — with no feature count in a header. You lose progress reporting and the row-range parallelism of step 6; convert to FlatGeobuf or GeoPackage first, which pays for itself immediately.
• The job dies at 90% with "no space left on device". GeoPackage grows a -wal journal alongside the file during a long append, so budget roughly double the final size on the output volume. The part-per-batch folder avoids the journal entirely.
• A resumed run duplicates features. The interrupted batch was half-committed to an appending sink. Only the part-per-batch layout is safely idempotent: a partial Parquet file is either replaced whole or absent, never merged.
• Output eastings are inf for a handful of records. Those points sit outside the target zone's valid range — usually a mis-keyed coordinate in an otherwise regional extract. Filter with np.isfinite per batch and write the rejects to a quarantine part rather than dropping them silently.
• Peak memory is flat but the run is far slower than expected. A grid-based datum shift is in the pipeline and PROJ is fetching grid tiles over HTTP per batch. Bake the grids into the image and set PROJ_NETWORK=OFF, or point PROJ_USER_WRITABLE_DIRECTORY at a persistent cache.

Frequently Asked Questions

How do I pick batch_size? Start from the memory you can spare, not from a round number of features. Estimate roughly 16 bytes per vertex plus about 200 bytes of object overhead per geometry, double it because to_crs holds a second copy, and size the batch to fit comfortably in a third of available RAM. For parcel-like polygons that lands near 250,000; for vertex-dense coastlines or river networks, ten thousand can be plenty. If peak RSS climbs across batches rather than plateauing, something is holding a reference — not a batch that is too large.

Why not just use Dask-GeoPandas for this? Because a straight reprojection has no shuffle in it. Every feature is independent, so the work is embarrassingly parallel and a scheduler buys you nothing but a dependency and a diagnostics port. Dask-GeoPandas earns its place when the same pass also joins, dissolves or aggregates across partitions, where the spatial shuffle is the hard part. For reproject-and-write, a batch loop plus a process pool is simpler and usually faster.

Can I let ogr2ogr or DuckDB do the reprojection instead? Yes, and for a plain format-to-format reprojection with no Python logic in the middle, ogr2ogr -t_srs EPSG:25832 is streaming, C-speed, and hard to beat. Reach for the Python loop when each batch needs something else done to it — filtering, attribute repair, validity fixes, quarantining rejects — or when the output must be partitioned by a rule GDAL cannot express. DuckDB's spatial extension can also ST_Transform while reading GeoParquet, covered in Querying GeoParquet with DuckDB Spatial.

Does any of this apply to rasters? The principle does, the code does not. A raster is reprojected by warping windows rather than batches of features, and the memory bound is the window size — the mechanics are in Coordinate Reference System Transformations and, for stacked time series, Reprojecting Raster Cubes with reproject_match.

Do I need to rebuild a spatial index on the output? For GeoPackage, GDAL maintains the R-tree as you append, so the file is queryable when the run ends. A folder of GeoParquet parts has no index at all — consumers rely on per-file bounding-box statistics instead, which only help if the parts are spatially coherent. If the parts will be queried by location rather than scanned whole, sort the source by a spatial key before the run so each part covers a compact area.