Optimizing Buffer Operations for Large Datasets in Python

Buffering millions of geometries in one GeoDataFrame.buffer() call routinely blows up with MemoryError or silently returns distances in degrees, and this guide fixes both. It is for anyone scaling proximity zones past the point where a single-threaded, single-allocation buffer fits in RAM. It sits under Proximity & Buffer Analysis in Spatial Analysis & Advanced Query Techniques, and it is the scaling companion to that guide's coverage of buffer cap and join parameters.

Why This Approach / What Goes Wrong

Two independent failures collide at scale, and treating them together is what makes the job tractable.

The first is coordinate units. A raw GeoPandas layer in geographic WGS84 (EPSG:4326) stores coordinates in degrees, so buffer(500) grows the geometry by 500 degrees of longitude — nonsense near the equator, catastrophically distorted near the poles. A buffer distance in metres is only meaningful in a projected, metric CRS. The naive tempation is to reproject to Web Mercator (EPSG:3857) because it is "already in metres," but Mercator's scale factor stretches with latitude, so a 500 m buffer in Oslo is measured in badly inflated units. Pick a local UTM zone via PyProj-backed CRS estimation instead.

The second is memory. buffer() is vectorized over Shapely 2.0's C geometry array, but it still allocates a brand-new geometry for every input row and holds the source column at the same time. At a few million polygons the peak resident set doubles, the allocator thrashes, and the process is OOM-killed. Splitting the frame into chunks caps peak memory to one chunk's worth of new geometries, and because each chunk is independent, the same split doubles as a unit of parallelism across CPU cores. Invalid input geometries — self-intersecting bowties, collapsed rings — raise TopologicalError mid-run and abort the whole job unless you clean them with Shapely's make_valid first.

There is a third failure that only appears at scale, and it is the one that turns a fifteen-minute job into an overnight one: the step after the buffer. Buffering is embarrassingly parallel and its cost grows linearly with row count, but dissolving or unioning the result is neither. union_all() over three million overlapping rings is a single-threaded GEOS operation whose cost grows super-linearly, and it cannot be chunked naively because rings that straddle a chunk boundary must be merged with their neighbours. If the pipeline ends in a union, plan for it explicitly: union within each spatial partition first, then union the (far fewer) partition results, and pass grid_size to snap coordinates so GEOS is not fighting floating-point noise on every shared edge.

import shapely   # buffered_chunks: list[GeoDataFrame], already in the metric CRS

# Chunk-local union first, then one cheap union of the partial results
partials = [chunk.geometry.union_all(grid_size=0.01) for chunk in buffered_chunks]
merged = shapely.union_all(partials, grid_size=0.01)

That only produces the right answer when the chunks are spatially coherent — a random row split scatters neighbouring rings across every chunk and defeats the optimisation entirely. Splitting by a grid cell, a region code or a geohash prefix is what makes the partial-union trick work, and the same partitioning logic is what Partitioning Strategies for Dask-GeoPandas formalises.

Prerequisites

• geopandas>=1.0 — vectorized buffer over the Shapely 2.0 geometry engine
• shapely>=2.0 — provides make_valid and C-backed buffer
• pyproj>=3.6 — drives estimate_utm_crs() for a correct metric projection
• numpy>=1.26 — array_split for chunking

pip install "geopandas>=1.0" "shapely>=2.0" "pyproj>=3.6" "numpy>=1.26"

Step-by-Step Implementation

1. Choose one metric CRS for the whole dataset

Estimate the UTM zone once, on the full frame, and reuse it for every chunk. Estimating per-chunk would pick different zones for chunks that happen to fall in different longitudes, making the outputs inconsistent.

import geopandas as gpd

parcels = gpd.read_file("parcels_nationwide.gpkg")
assert parcels.crs is not None, "Source layer has no CRS — set one before buffering"

# One projected, metric CRS for the entire job (e.g. EPSG:32633 for UTM 33N)
metric_crs = parcels.estimate_utm_crs()
print(f"[setup] buffering in {metric_crs.to_epsg()} ({metric_crs.name})")

2. Cut the work before you spend a core on it

The cheapest buffer is the one you never compute. Three reductions are almost always available, and together they routinely halve the run:

# a. Drop rows that cannot contribute to the answer
parcels = parcels[parcels.geometry.notna() & ~parcels.geometry.is_empty]

# b. Collapse exact duplicate geometries — common in merged municipal extracts
parcels["wkb"] = parcels.geometry.to_wkb()
before = len(parcels)
parcels = parcels.drop_duplicates("wkb").drop(columns="wkb")
print(f"[prep] {before:,} -> {len(parcels):,} after de-duplication")

# c. Thin the vertices of the *input*, not the output
#    A 0.5 m tolerance is invisible under a 500 m buffer but can halve GEOS work
parcels["geometry"] = parcels.geometry.simplify(0.5, preserve_topology=True)

Simplifying the input needs judgement: the tolerance must be small relative to the buffer distance, or the ring shape changes measurably. A useful rule is tolerance ≤ distance/100. Do it after reprojecting to the metric CRS so the tolerance is in metres, and never on a layer whose vertex positions are themselves the deliverable — cadastral boundaries, for instance, where the shared-edge behaviour matters and the safer approach is snapping and simplifying without gaps.

The fourth reduction is resolution. Buffer output vertex count scales linearly with it, and everything downstream — the concat, the Parquet write, any overlay — scales with the vertex count rather than the row count. Dropping from the default 8 to 4 on a screening layer of three million features removes tens of millions of coordinates from the pipeline at a cost of a few tenths of a percent of area.

3. Write a self-contained chunk worker

A ProcessPoolExecutor pickles the worker and its arguments to a fresh interpreter, so the function must not close over module-level state. It validates, projects to the metric CRS, buffers, and reprojects back to the caller's CRS.

from shapely.validation import make_valid


def buffer_chunk(
    chunk: gpd.GeoDataFrame,
    distance_m: float,
    metric_crs: str,
) -> gpd.GeoDataFrame:
    """Validate, project to metric CRS, buffer, project back."""
    if chunk.empty:
        return chunk

    source_crs = chunk.crs
    chunk = chunk.copy()
    chunk["geometry"] = chunk.geometry.apply(make_valid)

    projected = chunk.to_crs(metric_crs)
    projected["geometry"] = projected.buffer(
        distance_m, cap_style="round", join_style="round"
    )
    return projected.to_crs(source_crs)

4. Orchestrate the chunked, parallel run

np.array_split partitions the frame into n_chunks roughly equal slices. Each is submitted to a worker; a failed chunk falls back to its untouched input rather than aborting the batch. Keep the orchestration inside a main() function so it can live behind the __main__ guard that Windows requires.

import numpy as np
import pandas as pd
from concurrent.futures import ProcessPoolExecutor


def optimize_large_buffer(
    gdf: gpd.GeoDataFrame,
    distance_m: float,
    metric_crs: str,
    chunk_size: int = 100_000,
    max_workers: int = 6,
) -> gpd.GeoDataFrame:
    if gdf.empty:
        return gdf.copy()

    n_chunks = max(1, len(gdf) // chunk_size)
    chunks = np.array_split(gdf, n_chunks)
    print(f"[buffer] {len(gdf):,} features -> {len(chunks)} chunks, {max_workers} workers")

    buffered = []
    with ProcessPoolExecutor(max_workers=max_workers) as pool:
        futures = [
            pool.submit(buffer_chunk, chunk, distance_m, metric_crs)
            for chunk in chunks
        ]
        for i, future in enumerate(futures):
            try:
                buffered.append(future.result())
            except Exception as exc:  # noqa: BLE001 — log, keep the batch alive
                print(f"[buffer] chunk {i} failed ({exc}); keeping originals")
                buffered.append(chunks[i])

    return gpd.GeoDataFrame(pd.concat(buffered, ignore_index=True), crs=gdf.crs)

5. Run it and stream the result to disk

Write to GeoParquet so downstream steps read columns and geometry lazily instead of loading the whole result back into RAM.

def main() -> None:
    parcels = gpd.read_file("parcels_nationwide.gpkg")
    metric_crs = parcels.estimate_utm_crs()

    catchments = optimize_large_buffer(
        parcels, distance_m=500, metric_crs=metric_crs,
        chunk_size=100_000, max_workers=6,
    )
    catchments.to_parquet("parcel_catchments_500m.parquet")


if __name__ == "__main__":  # required for the ProcessPool on Windows/macOS spawn
    main()

That version still holds every buffered chunk in buffered before the final pd.concat, which reintroduces a full copy of the result at the end of an otherwise memory-bounded run. When the output is genuinely large, write each chunk as it lands and never assemble the whole frame:

from pathlib import Path


def buffer_to_dataset(gdf, distance_m, metric_crs, out_dir="catchments", **kw):
    """Stream each buffered chunk to its own Parquet file — peak RAM stays at one chunk."""
    out = Path(out_dir)
    out.mkdir(exist_ok=True)
    chunks = np.array_split(gdf, max(1, len(gdf) // kw.get("chunk_size", 100_000)))

    with ProcessPoolExecutor(max_workers=kw.get("max_workers", 6)) as pool:
        futures = {
            pool.submit(buffer_chunk, chunk, distance_m, metric_crs): i
            for i, chunk in enumerate(chunks)
        }
        for future in futures:
            i = futures[future]
            future.result().to_parquet(out / f"part-{i:05d}.parquet", compression="zstd")
    return out


# Downstream readers open the directory as one dataset, or one part at a time
catchments = gpd.read_parquet("catchments/")

A directory of Parquet parts is a first-class dataset for GeoPandas, DuckDB and Dask alike, so nothing downstream has to change — and it makes the job restartable, since a rerun can skip parts that already exist. Keep each part in the range of 50–200 MB; thousands of tiny files cost more in metadata than they save, which is the same sizing logic that governs GeoParquet storage layout.

6. Measure, then tune the two knobs

chunk_size and max_workers are the only levers, and guessing at them wastes more time than measuring. Instrument one representative chunk before launching the full run:

import resource, time

t0 = time.perf_counter()
sample = buffer_chunk(chunks[0], 500, metric_crs)
elapsed = time.perf_counter() - t0
peak_mb = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / 1024   # KB on Linux

print(f"[probe] {len(sample):,} rows in {elapsed:.1f}s, peak RSS {peak_mb:,.0f} MB")
print(f"[probe] projected total: {elapsed * len(chunks) / 6 / 60:.1f} min on 6 workers")
# [probe] 100,000 rows in 11.4s, peak RSS 412 MB
# [probe] projected total: 9.8 min on 6 workers

Multiply the measured peak by max_workers, add a gigabyte of headroom for the parent process, and compare against the machine's RAM. If the product exceeds it, halve chunk_size rather than reducing workers — the throughput loss from smaller chunks is small, while dropping a worker costs a proportional share of the whole run. On ru_maxrss, note the platform difference: Linux reports kilobytes and macOS reports bytes, so the divisor changes.

Verification

Confirm the run preserved every feature, kept the caller's CRS, and produced buffers with the geometrically expected area. A round-cap buffer of a point approaches π·r², so a spot check against that value catches a wrong or missing projection.

import math
import geopandas as gpd
from shapely.geometry import Point

# One point in WGS84; a 500 m buffer must cover ~pi * 500**2 square metres
probe = gpd.GeoDataFrame(geometry=[Point(10.75, 59.91)], crs="EPSG:4326")
metric_crs = probe.estimate_utm_crs()

result = optimize_large_buffer(probe, distance_m=500, metric_crs=metric_crs)

assert len(result) == 1
assert result.crs.to_epsg() == 4326                      # returned in source CRS
area_m2 = result.to_crs(metric_crs).area.iloc[0]
assert math.isclose(area_m2, math.pi * 500**2, rel_tol=0.01)
print(f"[ok] buffer area = {area_m2:,.0f} m2 (expected {math.pi * 500**2:,.0f})")
# [ok] buffer area = 785,145 m2 (expected 785,398)

Edge Cases & Debugging

• Datasets spanning many UTM zones: a single UTM zone distorts features thousands of km away. For continental extents swap estimate_utm_crs() for an equal-area CRS such as EPSG:6933 (World Cylindrical Equal Area), or partition by region and buffer each in its own zone.
• ProcessPoolExecutor still OOMs: each worker holds a full chunk plus its buffered copy. Drop chunk_size to 25_000 and max_workers to 2 — total peak memory is roughly chunk_size × max_workers × 2 geometries.
• make_valid returns a GeometryCollection: cleaning a self-intersecting polygon can yield mixed types; call .buffer(0) afterwards or filter with geometry.geom_type == "Polygon" before writing.
• Pickling overhead dominates on tiny geometries: for point layers the inter-process copy can cost more than the buffer itself — process in-line, or move to Dask-GeoPandas for true out-of-core partitioning instead of manual chunks.
• Windows spawns recursive subprocesses: without the if __name__ == "__main__": guard, every worker re-imports and re-runs the module. Linux/macOS default to fork and tolerate its absence, but keep the guard for portability.
• Every worker pins a core but wall time barely improves: GEOS releases the GIL, so the contention is elsewhere — usually the parent process serialising chunks faster than workers consume them, or the machine's memory bandwidth saturating. Confirm with the probe in step 6; if a single chunk already runs at 100% of one core, more workers cannot help.
• BLAS oversubscription starves the pool: libraries pulled in by numpy may start their own thread pool inside every worker, so six workers become six times N threads fighting for cores. Set OMP_NUM_THREADS=1 in the environment before starting the pool.
• Results differ subtly between runs: array_split is deterministic, but a chunk that fails and falls back to its untouched input leaves unbuffered rows in the output. Assert that every output geometry is a polygon before writing, rather than trusting the row count.
• Disk fills mid-run: buffered polygons are far larger than the input points or lines that produced them, and resolution multiplies that. Estimate output size from the probe chunk (sample.memory_usage(deep=True).sum()) before committing to a full run.

Frequently Asked Questions

When should I switch from manual chunking to Dask-GeoPandas? When the data no longer fits on one machine, when the pipeline has more than one stage that needs partitioning, or when you want the partitioning to persist between steps. Manual chunking is the right tool for a single embarrassingly parallel operation on a frame that fits in RAM — it has no scheduler, no cluster, and nothing to debug. Once a buffer is followed by a join and then a dissolve, the bookkeeping stops being worth it; move to Scaling with Dask-GeoPandas, which keeps spatial partitions and their bounds across the whole graph.

Is threading a cheaper alternative to processes here? Sometimes. Shapely 2.0 releases the GIL for vectorized operations including buffer, so a ThreadPoolExecutor does achieve real parallelism and avoids pickling chunks between processes entirely. That makes threads the better choice when geometries are small and numerous, where serialisation dominates. Processes remain safer when the worker does anything else — file I/O through GDAL, or a library whose thread safety you have not verified.

Does estimate_utm_crs() on a sample give the same zone as on the full frame? Not reliably. It picks the zone from the data's total bounds, so a sample that happens to omit the extremes can land one zone over. Compute it once on the full frame — it only needs the bounding box, which is cheap — and pass the resulting CRS explicitly into every worker, exactly as step 1 does. Deriving the CRS inside the worker is the single most common way to get inconsistent output.

Why write GeoParquet instead of a GeoPackage or Shapefile? Shapefile caps at 2 GB per file and truncates field names, so it is disqualified before performance enters the argument. GeoPackage is a fine single-file format but writes through a single SQLite connection, which becomes the bottleneck when six workers all want to append. Parquet parts are written independently, compress well with zstd, retain CRS metadata, and are read lazily column by column — which is the entire point of not materialising the frame again.