Preventing Memory Leaks When Processing Large GeoTIFF Rasters

TL;DR: Read GeoTIFF data in tile-aligned windows, assign None and del to every GDAL dataset after use, remove layers from QgsProject.instance() before deletion, and call gdal.FlushCache() plus gc.collect() at batch boundaries — or native C++ buffers will outlive their Python proxies and grow unbounded.

This page is part of the Memory Management and Garbage Collection for GIS Objects guide, itself a section of the PyQGIS Core Architecture & Data Handling reference.

Complete Runnable Template #

The function below processes any GeoTIFF in memory-safe chunks, applies a threshold operation, and guarantees deterministic cleanup whether or not an exception occurs. Drop it into a QGIS 3.x plugin or a standalone script — it has no dependency on a running QgsApplication.

import gc
from osgeo import gdal

def process_large_geotiff_chunked(
    input_path: str,
    output_path: str,
    threshold: float = 50.0,
    chunk_size: int = 0,  # 0 = auto-align to native tile size
) -> None:
    """
    Process a large GeoTIFF with memory-safe windowed I/O.

    Reads data in tile-aligned blocks, applies a threshold operation,
    and writes results to a new Cloud-Optimised GeoTIFF. Guarantees
    GDAL handle cleanup via explicit None assignment in a finally block.

    Args:
        input_path:  Absolute path to the source GeoTIFF.
        output_path: Absolute path for the destination file.
        threshold:   Pixel values above this become 1; others become 0.
        chunk_size:  Block dimension in pixels (0 = read from raster metadata).
    """
    gdal.SetCacheMax(256 * 1024 * 1024)  # cap at 256 MB for this process

    ds = gdal.Open(input_path, gdal.GA_ReadOnly)
    if ds is None:
        raise RuntimeError(f"GDAL could not open: {input_path}")

    band = ds.GetRasterBand(1)
    width, height = ds.RasterXSize, ds.RasterYSize
    geotransform = ds.GetGeoTransform()
    projection = ds.GetProjection()

    # Align to native tile boundaries to avoid cache thrashing
    if chunk_size == 0:
        tile_w, tile_h = band.GetBlockSize()
        chunk_size = max(tile_w, tile_h, 512)

    driver = gdal.GetDriverByName("GTiff")
    creation_opts = ["TILED=YES", f"BLOCKXSIZE={chunk_size}", f"BLOCKYSIZE={chunk_size}", "COMPRESS=DEFLATE"]
    out_ds = driver.Create(output_path, width, height, 1, band.DataType, creation_opts)
    out_ds.SetGeoTransform(geotransform)
    out_ds.SetProjection(projection)
    out_band = out_ds.GetRasterBand(1)

    try:
        for y in range(0, height, chunk_size):
            for x in range(0, width, chunk_size):
                w = min(chunk_size, width - x)
                h = min(chunk_size, height - y)

                chunk = band.ReadAsArray(x, y, w, h)
                if chunk is None:
                    continue

                processed = (chunk > threshold).astype(chunk.dtype)
                out_band.WriteArray(processed, x, y)

                # Release NumPy arrays immediately; do not accumulate in outer scope
                del chunk, processed
                gc.collect()  # force cyclic-reference sweep inside tight loops

        out_band.FlushCache()
    finally:
        # Assign None before del — SIP bindings may not honour del alone
        out_ds = None
        ds = None
        gc.collect()

Processing Flow #

The diagram below maps the decision points from opening the raster to final handle release. Every path through the loop ends with explicit dereferencing before the next iteration begins.

Architecture Breakdown #

Why Python Scope Exit Is Not Enough #

PyQGIS wraps GDAL through SIP/PyBind11 bindings. Python’s reference counter tracks the proxy object, but the underlying C++ GDALDataset, tile cache, and memory-mapped file handle live in native heap memory that the Python GC cannot see. When a QgsRasterLayer or gdal.Open() result falls out of scope, the proxy’s reference count reaches zero — but only if no other Python object holds a reference. QGIS’s internal working with QgsProject and layer registry keeps strong references to every layer added via QgsProject.instance().addMapLayer(). The vector and raster data access patterns guide documents a parallel effect for QgsVectorLayer providers — the same ownership rules apply here.

band.GetBlockSize() — Tile Alignment Is Not Optional #

GDAL reads raster data in native blocks (strips or tiles) encoded in the TIFF’s TileWidth / TileHeight tags. When your read window does not align with these boundaries, GDAL must decompress and cache overlapping blocks to satisfy even a single ReadAsArray call. This causes cache thrashing: tiles are loaded, partially used, evicted, then reloaded for the next misaligned window. On a 10 GB file with 256 × 256 native tiles, a 1000 × 1000 read window forces GDAL to cache up to 20 extra tiles per call.

Always query the native block dimensions first:

band = ds.GetRasterBand(1)
tile_w, tile_h = band.GetBlockSize()
# Use the larger dimension; fall back to 512 px if metadata is absent
chunk_size = max(tile_w, tile_h) or 512

QgsProject.instance().removeMapLayer() — Registry Detachment #

If you add a QgsRasterLayer to the project registry for rendering or processing, the registry holds a strong C++ reference that Python’s del statement cannot override. The sequence must be:

from qgis.core import QgsProject

layer_id = layer.id()
QgsProject.instance().removeMapLayer(layer_id)  # breaks registry reference
layer = None                                     # drop Python proxy
gc.collect()                                     # sweep any residual cycles

Omitting removeMapLayer is the single most common cause of persistent memory growth in QGIS plugin batch workflows.

gdal.SetCacheMax() — Bounding GDAL’s Block Cache #

GDAL’s global block cache defaults to 5 % of available RAM — often several gigabytes on a workstation. For batch pipelines processing many files sequentially, cap the cache before opening the first dataset:

from osgeo import gdal

# Set a hard 256 MB cap (value is in bytes)
gdal.SetCacheMax(256 * 1024 * 1024)

After each file, call gdal.FlushCache() to evict cached blocks from the previous dataset. Without the flush, GDAL retains tiles from closed datasets until the cache pressure limit forces eviction, which may be far into the next file’s processing.

Registration and Integration Patterns #

Standalone Script (no QgsApplication) #

Pure GDAL operations require no QgsApplication initialisation. Run the template directly:

# standalone_geotiff_processor.py
if __name__ == "__main__":
    process_large_geotiff_chunked(
        input_path="/data/dem_10m.tif",
        output_path="/data/dem_threshold.tif",
        threshold=200.0,
    )

Execute with python standalone_geotiff_processor.py — no QGIS installation needed beyond GDAL.

QGIS Plugin Integration #

When calling from inside a plugin, wrap in a QgsTask to keep the UI responsive. This is covered in detail in the asynchronous task execution with QgsTask guide. The memory management rules above apply unchanged inside QgsTask.run() — QgsTask does not alter GDAL’s ownership model.

from qgis.core import QgsTask, QgsApplication
import gc
from osgeo import gdal

class GeotiffProcessingTask(QgsTask):
    """Run windowed GeoTIFF processing on a background thread."""

    def __init__(self, input_path: str, output_path: str) -> None:
        super().__init__("Process GeoTIFF", QgsTask.CanCancel)
        self.input_path = input_path
        self.output_path = output_path

    def run(self) -> bool:
        try:
            process_large_geotiff_chunked(self.input_path, self.output_path)
            return True
        except Exception as exc:
            self.exception = exc
            return False
        finally:
            gc.collect()

task = GeotiffProcessingTask("/data/input.tif", "/data/output.tif")
QgsApplication.taskManager().addTask(task)

Subprocess Isolation for Enterprise Pipelines #

For workflows processing hundreds of files, process isolation is the safest strategy. Each worker process inherits a clean GDAL context, processes exactly one file, and terminates — giving the OS a guaranteed opportunity to reclaim all native memory regardless of Python reference cycles:

import subprocess
import sys
from pathlib import Path

def process_in_subprocess(input_path: str, output_path: str) -> None:
    """Spawn an isolated worker process; reclaims memory at OS level on exit."""
    script = Path(__file__).parent / "standalone_geotiff_processor.py"
    subprocess.run(
        [sys.executable, str(script), input_path, output_path],
        check=True,
    )

Production Best Practices #

• Windowed reads only. Never call ReadAsArray() without (xoff, yoff, xsize, ysize) arguments on files larger than available RAM.
• Assign None before del. SIP bindings do not always honour del alone; the None assignment breaks the C++ reference independently of Python’s destructor call order.
• Detach from QgsProject first. removeMapLayer() must precede any del or None assignment on QgsRasterLayer objects added to the registry.
• Cap GDAL_CACHEMAX at process start. A predictable upper bound prevents one large file from polluting the cache for subsequent files in a batch.
• Call gdal.FlushCache() between files. Do not rely on cache pressure to evict stale blocks — flush explicitly after closing each dataset.
• Align chunk_size to native tile dimensions. Use band.GetBlockSize() and choose a multiple of the larger dimension to eliminate cache thrashing.
• Use subprocess isolation for >50-file batches. OS-level reclamation is the only guarantee that eliminates all leak vectors, including Qt event-loop and QGIS singleton state.
• Test with tracemalloc. Wrap your batch loop with tracemalloc.start() / tracemalloc.take_snapshot() to confirm per-iteration memory delta is flat before deploying.

Parent Page and Related Topics #

This page is part of the Memory Management and Garbage Collection for GIS Objects guide.

Related:

• Memory Management and Garbage Collection for GIS Objects — parent: full GC and ownership model for all QGIS object types
• Vector and Raster Data Access Patterns — provider lifecycle rules that mirror the GDAL handle model shown here
• Optimizing Feature Iteration with QgsVectorLayer.getFeatures — windowed iteration for vector data, the raster pattern’s direct analogue
• Running Heavy Geoprocessing in the Background Without Freezing the UI — integrating memory-safe raster pipelines into a QgsTask worker