Geospatial Orchestration Hub

Offloading GDAL Work to a Process Pool

A ProcessPoolExecutor is the right home for GDAL work in a Python pipeline: each worker gets its own interpreter, its own GDAL configuration and its own memory, so the GIL is irrelevant and one worker’s failure does not take the others down. The two decisions that matter are sizing — by peak memory per worker, never by core count — and what crosses the boundary, which should be file paths and small parameter objects rather than arrays or frames. Get those right and the pool is close to free; get them wrong and it is slower than the sequential version while also running out of memory.

When to Use This Pattern

  • The pipeline is CPU-bound in a C library — rasterio.warp.reproject, shapely set operations, gdal.Translate through the Python bindings.
  • The work is per-partition and independent, so tiles or municipalities can be handed out with no coordination.
  • You are already async for the I/O half and need somewhere to put the compute half, as in async execution for heavy GIS tasks.
  • A subprocess would be awkward because the operation is a sequence of library calls rather than one command-line invocation.

Where the work is one command-line invocation, prefer an actual subprocess: it is equally parallel, it can be cancelled by signalling its process group, and its memory is bounded by the OS rather than by your pool’s arithmetic.

Complete Working Example

The pool below is created once for the flow, configured inside each worker, and fed with paths.

from __future__ import annotations

import os
from concurrent.futures import ProcessPoolExecutor
from dataclasses import dataclass
from pathlib import Path
from typing import Iterable


@dataclass(frozen=True)
class WarpJob:
    """Small, picklable, and it names files rather than carrying them."""

    source: Path
    dest: Path
    dst_crs: str = "EPSG:3857"
    resampling: str = "bilinear"
    nodata: float = -9999.0


def _init_worker() -> None:
    """Runs once per pool process, before any job.

    GDAL configuration set in the PARENT is not reliably inherited, and each
    worker needs its own cache budget — GDAL_CACHEMAX is per process, so four
    workers at the default 5% of RAM each is 20% of RAM in caches alone.
    """
    os.environ.setdefault("GDAL_CACHEMAX", "256")          # MB, per worker
    os.environ.setdefault("GDAL_NUM_THREADS", "1")          # the pool is the parallelism
    os.environ.setdefault("CPL_VSIL_CURL_CACHE_SIZE", "67108864")
    os.environ.setdefault("VSI_CACHE", "TRUE")


def warp_one(job: WarpJob) -> tuple[str, int]:
    """The unit of work. Imports are inside so the parent never loads rasterio."""
    import rasterio
    from rasterio.warp import Resampling, calculate_default_transform, reproject

    with rasterio.open(job.source) as src:
        transform, width, height = calculate_default_transform(
            src.crs, job.dst_crs, src.width, src.height, *src.bounds
        )
        profile = src.profile | {
            "crs": job.dst_crs, "transform": transform,
            "width": width, "height": height,
            "nodata": job.nodata,           # explicit: never inherit a missing nodata
            "compress": "DEFLATE", "tiled": True,
        }
        tmp = job.dest.with_suffix(".part")
        with rasterio.open(tmp, "w", **profile) as dst:
            for band in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, band),
                    destination=rasterio.band(dst, band),
                    resampling=Resampling[job.resampling],
                    num_threads=1,
                )
    os.replace(tmp, job.dest)               # promote only on success
    return str(job.dest), width * height


def run_pool(jobs: Iterable[WarpJob], max_workers: int = 4) -> list[tuple[str, int]]:
    # One pool for the whole run: starting a pool costs hundreds of milliseconds
    # per worker, most of it importing GDAL.
    with ProcessPoolExecutor(max_workers=max_workers, initializer=_init_worker) as pool:
        return list(pool.map(warp_one, jobs, chunksize=1))

Sizing the pool is arithmetic rather than judgement. Measure one job’s peak resident memory, add the per-worker GDAL cache, and divide the machine’s usable memory by the total:

def pool_size_for(peak_mb_per_job: float, gdal_cache_mb: float = 256.0,
                  usable_mb: float = 14_000.0, safety: float = 0.8) -> int:
    """Workers that fit, not workers that would be nice to have."""
    per_worker = peak_mb_per_job + gdal_cache_mb
    return max(1, int((usable_mb * safety) // per_worker))


# A 2.6 GB peak warp on a 16 GB worker: 3 processes, not 8 because it has 8 cores.
print(pool_size_for(peak_mb_per_job=2_600))
Pool size is a memory calculationOn a sixteen-gigabyte machine with a two-point-six-gigabyte peak per job plus a two-hundred-and-fifty-six-megabyte cache, three workers fit inside the safety margin and eight exceed the machine.3 WORKERS — sized by memory2.9 GB2.9 GB2.9 GBheadroom + OS8 WORKERS — sized by core count23.2 GB requested on a 16 GB machine. The OOM killer picks a worker at random, and the batch failssomewhere unrelated to the tile that caused it.
The eight-worker row is not slower — it is unreliable in a way that presents as a random tile failing, which is why it survives so long undiagnosed.

Parameter & Option Reference

Parameter Type Default Spatial notes
max_workers int memory ÷ peak Core count is the wrong input. A 2.6 GB warp on a 16 GB machine is three workers regardless of cores.
initializer callable _init_worker Where GDAL configuration belongs. Parent-process settings are not reliably inherited.
GDAL_CACHEMAX MB 256 Per process. The default is a percentage of total RAM, so N workers each claim it independently.
GDAL_NUM_THREADS int 1 The pool is the parallelism. Letting GDAL thread inside each worker oversubscribes the machine.
chunksize int 1 Tiles differ hugely in cost. Larger chunks batch expensive tiles onto one worker and idle the rest.
payload paths — A WarpJob is a few hundred bytes. An array or frame can be gigabytes, pickled twice.
num_threads in reproject 1 — Same reason as GDAL_NUM_THREADS: one axis of parallelism at a time.

Verification & Testing

Two properties are worth asserting: that the pool actually parallelises, and that memory stays inside the budget the sizing assumed.

import time
import psutil


def test_pool_parallelises(sample_jobs) -> None:
    """Four jobs across four workers should be much faster than four in sequence."""
    start = time.monotonic()
    [warp_one(job) for job in sample_jobs[:4]]
    sequential = time.monotonic() - start

    start = time.monotonic()
    run_pool(sample_jobs[:4], max_workers=4)
    pooled = time.monotonic() - start

    # Not 4x — pickling, start-up and I/O contention eat some — but clearly better.
    assert pooled < sequential * 0.6, f"no parallel speed-up: {pooled:.1f}s vs {sequential:.1f}s"


def test_peak_memory_within_budget(sample_jobs) -> None:
    proc = psutil.Process()
    children_peak = 0
    run = threading.Thread(target=run_pool, args=(sample_jobs, 3))
    run.start()
    while run.is_alive():
        total = sum(c.memory_info().rss for c in proc.children(recursive=True))
        children_peak = max(children_peak, total)
        time.sleep(0.2)
    run.join()
    assert children_peak < 12 * 1024**3, "pool exceeded the memory the sizing assumed"

On a running worker the check is simpler and worth having in the same start-up script that looks for orphaned processes. If the sum of the pool’s resident memory ever approaches the machine’s total, the sizing is wrong for the data currently flowing through — which changes when a new region with denser rasters arrives:

# Total RSS of every child of the flow process, in gigabytes.
ps -o rss= --ppid "$FLOW_PID" | awk '{s+=$1} END {printf "%.1f GB\n", s/1048576}'
The cost of the boundary depends on what crosses itA job description of a few hundred bytes pickles in microseconds. A two-gigabyte array pickles in seconds each way, often exceeding the cost of the computation itself.PATHS CROSS — the boundary is freeWarpJob · 300 Bworker reads the filereturns a path · 8 min of real workARRAYS CROSS — the boundary dominates2 GB arraypickle 14 s+ 2 GB copiedunpickle 14 s on the way backfor 8 min of real workTwenty-eight seconds of pure overhead per job, and three copies of the array alive at the peak.
The lower row also triples peak memory at the moment of transfer — parent copy, pickle buffer and worker copy — which is how a pool sized from a single-job measurement still runs out of memory.

Chunk size deserves a word too, because pool.map’s default is not one. With chunksize above one, jobs are handed out in fixed batches decided before any work starts — which is fine when every job costs the same and pathological when they do not. A batch of ten tiles that happens to contain the three most expensive tiles in the country will keep one worker busy long after the others have finished and gone idle. Spatial workloads are almost never uniform, so chunksize=1 and its slightly higher scheduling overhead is nearly always the right trade.

Why chunksize one, for non-uniform workWith chunksize ten, one worker receives the expensive tiles and runs long after the others finish. With chunksize one, work is redistributed continuously and all four workers finish together.chunksize=10 — batches fixed up frontw1w2w3w4three workers idle for the rest of the runchunksize=1 — dispatched as workers free upw1w2all four finish within a tile of each other
The total work is identical in both rows. The upper one takes almost twice as long because the expensive tiles were assigned before anyone knew they were expensive.

Common Pitfalls

  • Sizing by os.cpu_count(). The most common and most expensive mistake. GDAL work is memory-bound long before it is CPU-bound, and a pool sized from cores overcommits memory on every machine with more cores than it has tens of gigabytes.
  • Leaving GDAL_CACHEMAX at its default. The default is a fraction of total system RAM, applied independently in each worker. Four workers can therefore reserve four times what you expected before doing any work.
  • Configuring GDAL in the parent. gdal.SetConfigOption in the parent process does not reach forked or spawned workers reliably, and on spawn platforms it definitely does not. Use the pool’s initializer.
  • Creating a pool per batch. Each worker start-up imports GDAL and rasterio, which is hundreds of milliseconds. A pool created inside a loop spends more time starting than working.
  • Returning large objects. A worker that returns a numpy array pays the pickle cost on the way back, and the parent then holds every returned array at once. Return paths and counts.

Frequently Asked Questions

Why disable GDAL's own threading?

Because two layers of parallelism multiply. With GDAL_NUM_THREADS=ALL_CPUS inside four pool workers on an eight-core machine, thirty-two threads compete for eight cores and the context switching costs more than it gains. Pick one axis: either a pool of single-threaded workers, or one worker using all cores. The pool is usually better because it isolates failures and bounds memory per unit of work.

Does `fork` or `spawn` matter here?

Yes, more than it looks. On fork, workers inherit the parent’s memory lazily, which is cheap but means an already-large parent makes every worker large. On spawn, each worker is a fresh interpreter that must re-import GDAL, which is slower to start but far more predictable. For GDAL work spawn is usually the better default, and it is the only option on macOS and Windows anyway.

How do I get progress out of pool workers?

Return it, or write it. Workers cannot easily share a progress bar, but pool.map yields results in order as they complete, so counting completions in the parent is straightforward. For per-job detail, have the worker write a small JSON record next to its output; that also survives the worker dying, which an in-memory counter does not.

What happens when a worker dies?

ProcessPoolExecutor raises BrokenProcessPool on the next interaction, and the whole pool becomes unusable — including the jobs that were fine. That is harsh but honest: a worker killed by the OOM killer means the sizing was wrong, and continuing would kill more. Catch it at the flow level, dead-letter the batch, and let the dead-letter queue record which jobs were in flight.

← Async Execution for Heavy GIS Tasks

Related in this section