Fixing memory leaks in GDAL Python bindings

This page answers one question: why does a long-running process using GDAL’s Python bindings grow to gigabytes of resident memory, and how do you enforce the dataset-close discipline that stops it — with a CI leak gate that fails the build before the leak ships? It sits inside the Memory Management in Geospatial Extensions section of the Geospatial C-Extension Fundamentals & ABI Architecture reference, and gives you the ownership rules, the correct teardown order, and a valgrind/tracemalloc gate that catches regressions.

Why a dropped GDAL dataset reference leaks the native raster block cache Two columns compare teardown. On the left, a Python variable holding a GDAL Dataset goes out of scope but a lingering band reference keeps the C++ GDALDataset alive, so its block cache of raster tiles is never freed and native memory grows. On the right, explicit ds.Close after dropping the band reference destroys the GDALDataset, frees the block cache, and returns memory to the allocator. leak: lingering reference fixed: ordered teardown band = ds.GetRasterBand(1) del ds # band still holds it GDALDataset stays alive block cache of tiles never freed RSS climbs each iteration band = None # drop first ds.Close() # explicit GDALDataset destroyed block cache freed to allocator RSS flat across iterations

Context & Root Cause

GDAL’s Python bindings (the SWIG-generated osgeo.gdal module) wrap C++ objects whose lifetime is not governed by Python’s garbage collector in the way maintainers assume. A Dataset owns a block cache — decoded raster tiles held in native memory — that can reach hundreds of megabytes per open file. When you write the common pattern of opening a dataset in a loop and relying on del or scope exit to clean up, two things go wrong. First, any surviving reference to a child object (a Band, which internally holds its parent Dataset) keeps the whole C++ object graph alive, so the native cache is never released. Second, GDAL flushes writes only on explicit close, so a dropped write-mode dataset can both leak and corrupt output.

The result is a process whose Python heap looks healthy under tracemalloc while its resident set climbs relentlessly, because the leak lives in the native allocator arena, not the interpreter’s. This is the arena-boundary problem that Memory Management in Geospatial Extensions frames in general; here the concrete rule is that GDAL objects require ordered, explicit teardown, not reference-count roulette.

Solution / Fix

This targets GDAL 3.6–3.9 with the osgeo Python bindings. The rule set is small but strict.

1. Close datasets explicitly, in child-before-parent order 

from osgeo import gdal
gdal.UseExceptions()

def process(path):
    ds = gdal.Open(path)
    band = ds.GetRasterBand(1)
    stats = band.ComputeStatistics(False)
    band = None          # drop the child reference FIRST
    ds.Close()           # GDAL 3.7+: deterministic close + flush
    return stats

ds.Close() (added in GDAL 3.7) is the deterministic teardown entry point; on older versions assign ds = None only after every child (band, layer, feature) is already None. The ordering is non-negotiable because the child holds the parent.

2. Cap the block cache so a leak is bounded, not unbounded 

# Bound GDAL's global raster block cache (bytes). A hard cap turns a slow
# leak into a visible plateau you can alarm on.
gdal.SetCacheMax(256 * 1024 * 1024)   # 256 MB

3. Use context managers to make ordering automatic 

from contextlib import contextmanager

@contextmanager
def open_raster(path):
    ds = gdal.Open(path)
    try:
        yield ds
    finally:
        ds.Close()        # runs even on exception — no leaked handle

with open_raster("scene.tif") as ds:
    arr = ds.ReadAsArray()   # arr is a NumPy copy; safe after close

GDAL 3.8+ makes Dataset itself a context manager, so with gdal.Open(path) as ds: works directly. Either form guarantees the close runs on the exception path, which manual del does not.

Verification 

# 1. RSS must be flat across many iterations, not climbing
python - <<'PY'
import os, resource
from mymod import process
for i in range(500):
    process("scene.tif")
    if i % 100 == 0:
        rss = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss // 1024
        print(f"iter {i}: {rss} MB")
PY
# expected: RSS plateaus after warm-up; a rising staircase means a leak remains

# 2. Valgrind gate for native leaks (run on a native arch, not emulation)
PYTHONMALLOC=malloc valgrind --leak-check=full --error-exitcode=1 \
  python -c "from mymod import process; process('scene.tif')"
# expected: "definitely lost: 0 bytes" and exit code 0

PYTHONMALLOC=malloc disables pymalloc’s pooled allocator so Valgrind sees real malloc/free pairs — without it, false negatives hide the leak. Wire this command into CI as a nightly gate, as the checklist in Memory Management in Geospatial Extensions recommends.

Pitfalls & Alternatives

Relying on del ds alone. del only decrements the reference count; if any band, layer, or feature still references the dataset, the native object survives and the cache leaks. Drop children first, then close.

Holding a Band past its Dataset. A Band is a view into its parent’s memory. Keeping the band and closing the dataset is both a leak (the dataset stays alive) and a latent use-after-free once you do close it. Copy the data you need into a NumPy array and drop the band.

Assuming write flushes happen automatically. A write-mode dataset that is garbage-collected instead of closed may lose buffered blocks, producing a truncated GeoTIFF. Always Close() writers explicitly. For where these buffers sit relative to the interpreter heap and the PROJ transform buffers that share the arena, see the parent guide.

Two Heaps, One Process

The reason this leak is hard to see is that the process contains two independent memory regions with separate accounting, and the tools most Python developers reach for can only see one of them. tracemalloc, sys.getsizeof, gc.get_objects and every memory profiler built on them instrument CPython’s allocator. GDAL’s block cache, PROJ’s grid cache, GEOS’s geometry arenas and every buffer returned by a driver live in the C allocator, which those tools cannot observe at all.

Which memory each tool can see in a GDAL-using process The process is drawn as two regions. The CPython heap holds interpreter objects and small NumPy headers, and is visible to tracemalloc and gc. The native heap holds the GDAL block cache, PROJ grid cache, GEOS arenas and raster buffers, and is visible only to resident-set-size measurements, valgrind and massif. A Python wrapper object of a few hundred bytes can hold hundreds of megabytes in the native heap, which is why a healthy tracemalloc report coexists with a climbing resident set. CPython heap — what your profiler sees Dataset wrapper object · 96 bytes Band wrapper object · 88 bytes ndarray header · 112 bytes tracemalloc: 296 bytes retained — looks fine gc.collect() frees nothing that matters here native heap — where the leak lives GDAL block cache · up to SetCacheMax PROJ grid + context cache · tens of MB decoded raster buffers · hundreds of MB visible only to RSS, valgrind, massif, jemalloc stats freed by Close(), not by the garbage collector owns a 96-byte wrapper can pin 400 MB — which is why the leak is invisible to every tool that counts Python objects and why the only reliable signal is resident set size measured across many iterations of the same work

That asymmetry has a practical consequence for how you look for these leaks. Measuring RSS in a loop is not a crude approximation of a better technique — it is the correct technique, because it is the only measurement that spans both regions. Run the same operation several hundred times, sample RSS every fifty iterations, and look for a line rather than a plateau. A plateau at the cache cap means the cache filled and stopped, which is healthy. A line means something is not being released.

The second consequence is that gc.collect() is not a remedy. Even when the collector does break a reference cycle and destroy a wrapper, GDAL’s C++ destructor only runs if the binding’s __del__ path reaches it, and for datasets opened in write mode a destructor-driven close does not guarantee the flush ordering you want. Explicit Close() is not defensive style; it is the only mechanism with defined semantics.

Reading the RSS Curve

The shape of the resident-set curve, not its height, tells you which of the three common situations you are in. Sampling every fifty iterations for a few hundred iterations is enough to distinguish them at a glance.

Three resident-set curves and what each one means A line chart with iterations on the horizontal axis and resident set size on the vertical axis. A healthy curve rises while the block cache fills and then plateaus flat. A leaking curve rises steadily without ever flattening. A sawtooth curve rises and drops repeatedly, indicating an unbounded cache being trimmed under memory pressure rather than a leak. RSS iterations healthy — cache fills, then flat leaking — never flattens sawtooth — unbounded cache trimmed under pressure, not a leak warm-up window: ignore the first ~50 iterations

Distinguishing the second from the third matters because the remedies differ: a genuine leak needs an ownership fix in your code, while a sawtooth almost always means SetCacheMax was never called and GDAL is sizing its cache as a fraction of available memory.

A Leak Gate That Fails the Build

Leak regressions are easy to reintroduce because the code that causes them looks correct — a helper that returns a Band for convenience, a cache that holds a Dataset “just for its metadata”. A CI gate catches them at review time rather than in production, and it needs no special tooling.

# tests/test_no_leak.py — fails when RSS grows beyond a tolerance
import resource
import pytest
from mymod import process

def rss_mb():
    kb = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss
    return kb / 1024  # Linux reports KiB; macOS reports bytes — normalise per platform

@pytest.mark.slow
def test_rss_is_flat(tmp_raster):
    for _ in range(50):          # warm the caches first
        process(tmp_raster)
    baseline = rss_mb()
    for _ in range(500):
        process(tmp_raster)
    growth = rss_mb() - baseline
    assert growth < 64, f"resident set grew {growth:.0f} MB over 500 iterations"

Three details make the difference between a useful gate and a flaky one. Warm up before taking the baseline, so the first measurement is after the block cache has reached steady state rather than while it is filling. Choose a tolerance in the tens of megabytes rather than zero, because allocator fragmentation and lazily-initialised PROJ state produce real, bounded growth that is not a leak. And cap the cache with gdal.SetCacheMax in the test’s setup, so the plateau is deterministic instead of scaling with whatever memory the runner happens to have.

For a sharper signal, valgrind --leak-check=full --show-leak-kinds=definite run over a short script catches genuinely unreachable native allocations and names the C call site. It is far too slow for every CI run, but a weekly scheduled job is affordable and it finds the class of leak that RSS-watching misses — a small allocation leaked per call, which takes hundreds of thousands of iterations to become visible as RSS but shows up immediately as a definite loss.

Frequently Asked Questions

Does del ds do the same thing as ds.Close()?

No. del removes one reference; the object is destroyed only when the last reference goes, and any live child object holds one. Even then, destruction runs the binding’s finaliser rather than a documented close path, which for a write-mode dataset does not guarantee that pending blocks are flushed before the file handle is released. Close() is explicit, ordered and testable.

Why does memory keep growing even with every dataset closed?

Two usual suspects. The first is the block cache doing its job: it is meant to retain decoded tiles across dataset lifetimes, so growth up to SetCacheMax is expected and stops there. The second is PROJ, which caches grids and coordinate operations per context; a long-running process that builds many transformers accumulates them unless transformers are reused or the context is periodically reset.

Is this different under multiprocessing?

It is more forgiving, because a worker process that exits returns everything to the operating system regardless of what leaked. That is why batch pipelines with short-lived workers rarely notice these bugs, and why the same code in a long-lived service falls over within a day. If a service must use GDAL heavily, recycling workers after a fixed number of tasks is a legitimate and widely used mitigation.

Does NumPy’s ReadAsArray keep the dataset alive?

Not in current GDAL: ReadAsArray copies into a new NumPy array, so the array is independent of the dataset and safe to use after Close(). Zero-copy paths — memory-mapped access or a buffer that wraps native memory — are the exception, and there the array does pin the native allocation, which makes the copy-versus-view distinction worth knowing before you optimise it away.