A synthetic grid that fits comfortably in a notebook at 1,000 × 1,000 cells becomes an out-of-memory crash at 200,000 × 200,000 — the resolution a continental land-cover, elevation, or population surface actually demands. This page is part of Spatial Distribution & Pattern Generation: where the other generators in this area decide the statistics of what a grid contains, memory-optimized generation is the engineering that lets you produce it at real scale without the peak resident set exceeding the machine. The central discipline is simple to state and hard to hold: peak memory must stay flat as the grid grows. A single float64 array covering a continent at 30-metre resolution is tens of terabytes — it was never going to be a NumPy array in RAM. The generator has to produce, validate, and persist the grid in bounded windows, touching only a working-set slice at any instant.
The generator holds one block-aligned window in memory, writes it out of core to a tiled GeoTIFF, and releases it — so peak RAM is set by the window, not the grid.
The failure this page solves is the materialize-then-write anti-pattern — building the entire synthetic grid as one dense in-memory array and only afterward writing it to disk. It works flawlessly in development at small sizes and then dies the first time it meets production resolution. The crash is the visible symptom; the underlying defects are subtler and each has its own fix.
Whole-array allocation.np.zeros((h, w), dtype="float64") for a continental grid asks for tens of terabytes in one call. Even at 30-metre resolution over a single country, an intermediate float64 array can exceed physical RAM by orders of magnitude, and the allocation fails before a single value is generated.
Hidden intermediate copies. Even when the target fits, a chain of NumPy operations that each returns a new array can transiently hold three or four full copies. Peak memory is the sum of the live intermediates, not the size of the result, so a grid that “fits” still OOMs mid-pipeline.
Block-misaligned writes. Writing arbitrary rectangular windows into a tiled raster forces the format driver to read, modify, and rewrite whole internal blocks — read-modify-write amplification that both slows the write and inflates memory as partial blocks are buffered.
Accumulate-then-persist. Collecting every generated tile in a Python list to concatenate at the end recreates the whole-array problem one level up: the list holds the entire grid before the single terminal write, so the concatenation is the OOM point.
The remedy is a single principle applied everywhere: never hold more than a bounded working set. Generate into windows sized to the output’s internal blocks, validate each window in isolation, write it out of core, and release it. Peak memory becomes a function of window size and worker count — both of which you choose — and is independent of the grid’s total extent. This is the same reproducible-tiling discipline that governs async execution for large grids; that page focuses on concurrency and determinism, while this one focuses on the memory envelope.
Large-grid generation leans on the chunked-array and windowed-raster stacks. Pin majors so chunking heuristics and GeoTIFF block behaviour stay stable across environments:
Two environment decisions dominate the memory envelope. First, align your generation window to the output raster’s internal block size. A cloud-optimized GeoTIFF is internally tiled (commonly 256 × 256 or 512 × 512); if your write window is not a whole number of blocks on both axes, GDAL must read the surrounding blocks, merge, and rewrite them, which both slows the write and quietly raises peak memory. Second, set the dtype deliberately. A grid stored as float64 costs eight bytes a cell; the same field as float32, int16, or a scaled uint8 with a nodata value costs half, a quarter, or an eighth — and the working-set window shrinks proportionally. Choose the narrowest dtype the attribute’s precision allows and carry it end to end, rather than generating in float64 and downcasting at the write.
The mechanism that makes an unbounded grid tractable is the working set: at any instant the process holds only the windows currently being generated, validated, or written, never the whole grid. Everything else is either not yet produced or already flushed to the out-of-core sink. Three cooperating ideas realize this.
Windowed (tiled) processing partitions the target extent into a grid of windows and processes one — or a bounded pool of a few — at a time. Each window is a self-contained generation task with its own coordinate offset, so the synthetic values it produces depend only on its position, not on any other window’s state. Because windows are independent, memory is bounded by window_cells × dtype_bytes × concurrent_windows, a number you set directly and that does not grow with the extent.
Chunked arrays are the in-memory realization of that idea. Rather than a dense NumPy array, a Dask array is a lazy grid of chunks; operations build a task graph and materialize only the chunks needed for the block currently being written. Set the chunk shape equal to the output block shape and the array never holds more than a few blocks of live data even as the logical array spans a continent. This is what lets an ordinary array expression — add noise, apply a spatial gradient, threshold into classes — run over a grid far larger than RAM without a rewrite; the same chunk-tuning trade-offs govern the density mapping and heat generation surfaces that are produced the same way.
Out-of-core writes close the loop. Each finished window is written straight into its block position in the output raster with a windowed write and then dropped from memory. The grid accumulates on disk, not in RAM. When the sink is a cloud-optimized format, the write also builds the internal tiling and overview pyramid that let a consumer later read any sub-region without downloading the whole file — the subject of the drill-down on streaming large grids to cloud-optimized GeoTIFF, which covers the block ordering and overview-generation passes that keep the write single-pass.
The one subtlety is window interdependence. Purely local generators — per-cell noise, elevation from a coordinate function — need no context beyond the window. Generators with spatial extent — a smoothing kernel, a diffusion, a distance transform — need a halo of neighbouring cells, or the window edges show seams. Read or generate a buffer margin of at least the kernel radius, compute, then trim the halo before writing, exactly as a tiled tessellation trims its buffered seams. The same discipline underpins the polygon tessellation algorithms when their cells are produced tile by tile.
The pipeline runs as four auditable stages. Each block assumes the pinned toolchain and a metric CRS.
Step 1 — Define a block-aligned tiling of the extent. Compute window offsets that are whole multiples of the output block size so every write lands on block boundaries.
python
import rasterio
from rasterio.windows import Window
BLOCK =512# match the output GeoTIFF internal tile sizedefblock_windows(width:int, height:int, block:int= BLOCK):for row_off inrange(0, height, block):for col_off inrange(0, width, block):
w =min(block, width - col_off)
h =min(block, height - row_off)yield Window(col_off, row_off, w, h)# block-aligned by construction
Step 2 — Generate each window deterministically into a narrow dtype. Derive the RNG seed from the window’s integer offsets plus a pipeline version hash, never from wall-clock time, so a window’s values are identical no matter when or in what order it runs.
python
import numpy as np
VERSION =7# bump when the generation logic changesdefgenerate_window(win: Window, transform)-> np.ndarray:# Deterministic per-window seed: stable across runs, machines, and ordering.
seed =(VERSION *73856093)^(win.col_off *19349663)^(win.row_off *83492791)
rng = np.random.default_rng(seed &0xFFFFFFFF)# Example field: a smooth spatial gradient plus bounded noise, in float32.
rows, cols = np.mgrid[0:win.height,0:win.width]
xs = transform.c +(win.col_off + cols)* transform.a
ys = transform.f +(win.row_off + rows)* transform.e
field =(0.001* ys).astype("float32")
field += rng.normal(0.0,0.05, size=field.shape).astype("float32")return field # one window only; released after the write below
Step 3 — Stream each window straight to a cloud-optimized GeoTIFF. Open the sink once with tiling and compression enabled, then write windows in place; peak memory holds only the active window plus GDAL’s block buffer.
python
defwrite_grid(path:str, width:int, height:int, crs, transform)->None:
profile ={"driver":"GTiff","dtype":"float32","count":1,"width": width,"height": height,"crs": crs,"transform": transform,"tiled":True,"blockxsize": BLOCK,"blockysize": BLOCK,"compress":"deflate","predictor":2,"BIGTIFF":"YES",}with rasterio.open(path,"w",**profile)as dst:for win in block_windows(width, height):
data = generate_window(win, transform)# bounded working set
dst.write(data, window=win, indexes=1)# out-of-core windowed writedel data # release before next window
Step 4 — Build overviews in a single trailing pass. After the full-resolution write, add a decimated overview pyramid so consumers can read coarse zoom levels cheaply. GDAL streams the overview build block by block, so it too stays memory-bounded.
A large-grid generator needs two kinds of gate: one proving the memory envelope stayed bounded, and one proving the content is correct and reproducible despite being produced piecewise.
python
import numpy as np
import rasterio
import tracemalloc
defgate_memory_envelope(build_fn, budget_mb:float)->None:"""Assert peak Python allocation stays under a flat budget during a build."""
tracemalloc.start()
build_fn()# runs the full windowed generation
_, peak = tracemalloc.get_traced_memory()
tracemalloc.stop()assert peak /1e6< budget_mb,f"peak {peak/1e6:.0f} MB exceeded {budget_mb} MB"defgate_seams_and_determinism(path_a:str, path_b:str)->None:with rasterio.open(path_a)as a, rasterio.open(path_b)as b:# 1. Two independent runs are byte-identical (deterministic per-window seeds).for win in a.block_windows(1):
wa = a.read(1, window=win[1])
wb = b.read(1, window=win[1])assert np.array_equal(wa, wb),"non-deterministic window output"# 2. No seam discontinuity: cross-block gradient stays within tolerance.
full = a.read(1)
col_jump = np.abs(np.diff(full, axis=1))[:, BLOCK -1::BLOCK]assert col_jump.max()<1.0,"visible seam at a block boundary"
The memory-envelope gate is the load-bearing check unique to this page: it runs the whole build and asserts peak allocation stays under a flat budget that does not scale with grid size — the direct test that the working-set discipline held. Wire both gates as build-failing checks through the CI/CD integration pipeline, running the memory gate on a scaled-down but same-shaped extent so CI never itself needs the full grid.
Once memory is bounded, throughput becomes the constraint, and the knobs interact.
Window size is a memory-versus-overhead trade. Larger windows amortize per-window scheduling and GDAL open costs but raise the resident working set; smaller windows shrink peak memory but multiply overhead. Size the window to a small integer number of output blocks and profile, rather than guessing.
Bound the worker pool explicitly. Peak memory is window_bytes × concurrent_windows, so an unbounded pool defeats the whole design. Cap concurrency to a fixed pool sized to your RAM budget; the async execution for large grids layer covers the scheduler that keeps that pool saturated without oversubscription.
Compression trades CPU for I/O and disk. Deflate with a horizontal predictor shrinks smooth synthetic fields dramatically, cutting write bandwidth and final size at a modest CPU cost — usually a net win when the sink is remote object storage.
Halo cost scales with kernel radius. A generator needing a neighbour halo reads or regenerates extra cells per window; keep the kernel radius small relative to the window, or the halo overhead dominates. For large kernels, widen the window so the halo is a small fraction of it.
Avoid the concatenation trap. Never collect finished windows to concatenate at the end — that reintroduces the whole-array allocation this design exists to remove. Each window’s only destination is the out-of-core sink.
My generator OOMs at production resolution but runs fine in development. What changed?
You are materializing the whole grid as one in-memory array before writing it, which only fits at small development sizes. Switch to windowed generation: split the extent into block-aligned windows, generate one bounded window at a time, write it out of core to the raster, and release it. Peak memory then depends on window size and worker count, both of which you set, not on the total grid extent.
Peak memory is far higher than the size of my output array. Why?
Chained NumPy operations that each return a new array can hold several full copies live at once, so peak memory is the sum of the intermediates, not the size of the result. Use chunked Dask arrays whose chunk shape matches the output block size so only a few blocks materialize at a time, and prefer in-place operations. The peak then tracks the working-set window rather than the whole expression.
Writes are slow and memory spikes even though each window is small.
Your write windows are almost certainly not aligned to the raster's internal block size, so GDAL must read the surrounding blocks, merge your partial data, and rewrite them — read-modify-write amplification. Make every window a whole number of internal tiles on both axes by using the same block size for the tiling scheduler and the GeoTIFF blockxsize and blockysize, and the writes become clean in-place block writes.
Visible seams appear at the boundaries between windows.
A generator with spatial extent — a smoothing kernel, diffusion, or distance transform — needs neighbouring context that a bare window does not provide, so its edges disagree with the next window. Generate a halo margin of at least the kernel radius around each window, compute over the padded window, then trim the halo before writing. Purely local per-cell generators need no halo and never show seams.
Two runs of the same configuration produce different grids.
The per-window RNG is drawing from wall-clock time or a global counter that depends on execution order, so concurrent windows desynchronize. Derive each window's seed deterministically from its integer offsets and a pipeline version hash, so a window yields identical values regardless of when or in what order it runs. Add a CI test that builds the grid twice and asserts the two rasters are byte-identical.
