Reprojecting Large Raster Stacks Without Memory Spikes #

The scenario: rioxarray.reproject() on an 8,760-band hourly GHI stack raises a MemoryError on a 64 GB machine, and the obvious fix — a bigger machine — buys one more band before it fails again. The problem is not the machine; it is that a naive warp holds the source array, the destination array and an intermediate simultaneously, so peak memory is roughly three times the larger of the two. This page reprojects the same stack in bounded memory, and it is the scaling detail behind coordinate reference systems for energy projects.

Root-cause analysis #

Three compounding causes turn a routine warp into an out-of-memory failure.

1. Whole-array semantics. reproject on an in-memory array allocates the destination before it writes into it, and the resampling kernel needs the source resident at the same time. For a 61,000 by 58,000 float32 band that is 14.2 GB twice over before any intermediate.
2. A destination grid derived implicitly. When the target transform and shape are not specified, the warp computes them from the source bounds, which can produce a destination substantially larger than the source — a rotation of a few degrees expands the bounding box, and a poorly chosen target resolution can multiply cell count.
3. Band-major iteration on a band-minor file. Reading band by band from a file chunked by tile means every band read touches every tile, so the I/O cost multiplies by the band count even when memory is under control.

Pre-flight validation #

The peak is computable before the run, and computing it is faster than discovering it.

import numpy as np
import rasterio
from rasterio.warp import calculate_default_transform


def estimate_warp_memory(
    src_path: str,
    dst_crs: str,
    *,
    dst_resolution: float | None = None,
    window_px: int = 2048,
) -> dict:
    """Peak resident bytes for a whole-array warp versus a windowed one."""
    with rasterio.open(src_path) as src:
        transform, width, height = calculate_default_transform(
            src.crs, dst_crs, src.width, src.height, *src.bounds, resolution=dst_resolution
        )
        itemsize = np.dtype(src.dtypes[0]).itemsize
        src_band = src.width * src.height * itemsize
        dst_band = width * height * itemsize
        window_bytes = window_px * window_px * itemsize
        return {
            "bands": src.count,
            "src_band_gb": src_band / 1e9,
            "dst_band_gb": dst_band / 1e9,
            "whole_array_peak_gb": (src_band + dst_band) * 1.2 / 1e9,
            "windowed_peak_gb": window_bytes * 4 * 1.2 / 1e9,
            "dst_shape": (height, width),
            "dst_transform": transform,
        }

The 1.2 factor is the intermediate and bookkeeping overhead; the point of the function is not precision but the ratio, which is typically three to four orders of magnitude.

Fix implementation #

The fix has two halves: define the destination grid explicitly, then stream windows into it. Writing directly to a rasterio dataset means the destination never has to be resident either.

import rasterio
from rasterio.warp import Resampling, calculate_default_transform, reproject
from rasterio.windows import Window


def reproject_stack_windowed(
    src_path: str,
    dst_path: str,
    dst_crs: str,
    *,
    dst_resolution: float | None = None,
    resampling: Resampling = Resampling.bilinear,
    window_px: int = 2048,
    compress: str = "LZW",
) -> dict:
    """Warp every band through bounded-memory windows into a tiled GeoTIFF."""
    with rasterio.open(src_path) as src:
        transform, width, height = calculate_default_transform(
            src.crs, dst_crs, src.width, src.height, *src.bounds, resolution=dst_resolution
        )
        profile = src.profile.copy()
        profile.update(
            crs=dst_crs, transform=transform, width=width, height=height,
            tiled=True, blockxsize=512, blockysize=512, compress=compress,
            BIGTIFF="IF_SAFER",
        )

        written = 0
        with rasterio.open(dst_path, "w", **profile) as dst:
            for band in range(1, src.count + 1):
                for row in range(0, height, window_px):
                    for col in range(0, width, window_px):
                        w = Window(col, row,
                                   min(window_px, width - col),
                                   min(window_px, height - row))
                        dst_arr = rasterio.band(dst, band)
                        reproject(
                            source=rasterio.band(src, band),
                            destination=dst_arr,
                            src_transform=src.transform,
                            src_crs=src.crs,
                            dst_transform=dst.window_transform(w),
                            dst_crs=dst_crs,
                            dst_nodata=src.nodata,
                            resampling=resampling,
                            num_threads=2,
                        )
                        written += 1
        return {"bands": src.count, "windows_written": written,
                "dst_shape": (height, width), "dst_crs": dst_crs}

Specifying tiled=True on the destination is not cosmetic: an untiled (striped) GeoTIFF forces every windowed read afterwards to touch whole rows, which undoes the memory discipline at the next stage.

Fallback routing and performance tuning #

• Set GDAL_CACHEMAX explicitly. GDAL’s block cache defaults to a share of RAM and will grow to fill it; a 512 MB cap keeps the peak predictable and costs almost nothing in throughput.
• Prefer VRT for a purely lazy warp. gdal.BuildVRT plus a warped VRT gives a virtual reprojected dataset that materialises nothing until read — ideal when only a subset will ever be consumed.
• Choose the window from the destination, not the source. Windows are written in destination space; sizing them from source tiles produces uneven work per window.
• Reproject once, reuse many times. A warped national product is expensive and static; a pipeline that warps on every run is paying that cost repeatedly for an artefact that could be cached.
• Watch the resampling kernel’s read amplification. Lanczos reads a 6 by 6 source window per destination cell, so its I/O is 36 times nearest’s — which on a network-backed source dominates everything else.

Downstream validation #

import numpy as np
import rasterio


def assert_warp_integrity(src_path: str, dst_path: str, *, sample_px: int = 512) -> None:
    """Cheap post-warp checks that catch the failures a MemoryError would have hidden."""
    with rasterio.open(src_path) as src, rasterio.open(dst_path) as dst:
        assert dst.count == src.count, f"band count changed: {src.count} → {dst.count}"
        assert dst.dtypes[0] == src.dtypes[0], f"dtype changed: {src.dtypes[0]} → {dst.dtypes[0]}"
        assert dst.nodata == src.nodata, "nodata value was not carried through the warp"
        assert dst.crs.to_epsg() is not None, "destination CRS did not resolve to an EPSG code"

        w = rasterio.windows.Window(0, 0, min(sample_px, dst.width), min(sample_px, dst.height))
        sample = dst.read(1, window=w, masked=True)
        assert sample.count() > 0, "top-left window is entirely nodata — check the destination bounds"
        finite = sample.compressed()
        assert np.isfinite(finite).all(), "non-finite values introduced by the warp"

Choosing the destination grid deliberately #

The most consequential decision in a warp is one that is usually left to a default: what the destination grid should be. calculate_default_transform proposes a grid that covers the reprojected footprint at a resolution derived from the source, and its proposal is frequently wrong in two ways.

It is wrong in extent when the analysis only needs a study area. A national source reprojected in full produces a national destination, most of which will be clipped away immediately — and the warp paid for every cell. Passing explicit destination bounds cuts both the write and the read.

It is wrong in resolution when the source and target frames have different units or the reprojection stretches one axis. A 0.0417-degree source becomes something like 4,632 metres in a metric frame, and rounding that to a convenient 5,000 or 4,000 metres is usually preferable — an awkward resolution propagates into every downstream product and makes alignment with other layers harder than it needs to be.

Set both explicitly, record them with the output, and the warp becomes reproducible: two people warping the same source with the same declared grid get byte-identical results, which is not true when both rely on a default that depends on the source extent.

Frequently asked questions #

Should I reproject the raster or the vector? #

Almost always the vector. Moving a few thousand geometries is microseconds; resampling a few billion pixels is minutes and lossy. Reproject the raster only when the analysis is raster-on-raster and the two grids genuinely have to align — and then reproject the smaller one.

Does dask solve this? #

It manages it rather than solving it. rioxarray with a Dask-backed array will chunk the warp and keep peak memory bounded, which is the same win as windowing, with a scheduler attached. For a single machine the windowed loop above is simpler and has fewer failure modes; Dask earns its complexity when the work genuinely spans machines.

What window size should I use? #

Large enough to amortise the per-window overhead and small enough to keep several windows in flight comfortably — 1,024 to 4,096 pixels square is the usual range for float32. Below about 512 the per-call cost starts to dominate; above 8,192 the memory advantage erodes.

Why did the destination come out larger than the source? #

Because a reprojection rotates the footprint, and the axis-aligned bounding box of a rotated rectangle is larger than the rectangle. Reprojecting a UTM tile to a conic frame can add 10 to 20 percent of cells, most of them nodata. Specifying the destination bounds explicitly — clipped to the study area — avoids paying for that.

How do I keep the band metadata? #

Copy it explicitly. rasterio carries the dataset profile but band descriptions and per-band tags are not part of it, so an hourly stack that loses its timestamps in the warp becomes an anonymous cube. Reading src.descriptions and src.tags(band) and writing them onto the destination costs two lines and preserves the only thing that makes the stack interpretable.

Is it worth compressing the output? #

Yes, with a lossless codec. LZW with a horizontal predictor typically halves a float32 raster at negligible read cost, and the saving compounds through every downstream read. What is not worth it is a lossy codec on data that feeds an arithmetic chain — the artefacts are small, systematic and impossible to distinguish from signal later.

Can a warp be resumed after a failure? #

Yes, if the destination is written window by window and the windows are addressable. Recording which destination windows have been written — a small sidecar or a per-window checksum — lets a rerun skip the completed ones, which on a multi-hour national warp turns a crash from a full restart into a few minutes of catch-up. Without that record the only safe action is to start again, because a partially written GeoTIFF is indistinguishable from a complete one.

Does the source need to be a Cloud-Optimised GeoTIFF? #

Not required, but it changes the economics substantially when the source is remote. A COG’s internal tiling and overviews mean a windowed read fetches kilobytes rather than the whole file, so the windowed strategy above is efficient over HTTP as well as over local disk. A striped, uncompressed GeoTIFF on object storage forces each window read to fetch whole rows, which is where most of the “windowing did not help” reports come from.

Related #

• Coordinate Reference Systems for Energy Projects — the parent workflow
• Solar Irradiance Raster Processing — the stacks this technique is usually applied to
• Resampling & Raster Kernel Quick Reference — choosing the kernel the warp uses
• Mosaicking Tiled GHI Rasters with rasterio.merge — the stage that usually precedes a warp