When a synthetic raster is larger than RAM, the only way to write it as a valid cloud-optimized GeoTIFF is to generate and write it in block-aligned windows and build overviews without ever materializing the full array, so peak memory stays flat whether the grid is 10,000 or 1,000,000 pixels on a side.
Part of Memory-Optimized Large-Grid Spatial Generation: that page covers keeping generation within a memory envelope in general; this page drills into the output stage — writing directly to a cloud-optimized GeoTIFF (COG) with windowed, tiled writes and internal overviews — where the failure is not the generation algorithm but the naive write() that tries to hold the whole band in memory and either OOM-kills the job or produces a striped, non-COG file that no range-request reader can use efficiently.
A GeoTIFF is COG-valid only when its data is organized into internal tiles (not scanline strips) and it carries a pyramid of reduced-resolution overviews, with the image file directory laid out so a client can fetch any tile or overview level with a single HTTP range request. The memory problem and the layout problem are linked: the natural way to satisfy “tiled + overviews” is to build the array in memory and let the driver reorganize it on close, and that array is exactly what does not fit.
Consider the arithmetic. A single-band float32 grid of N×N pixels needs 4N2 bytes. At N=100,000 that is 40 GB for one band before any copy — and the standard “build then write” path needs at least two such buffers live at once (the source array and the driver’s tiled reorganization), plus a third to compute each overview level by decimating the full-resolution array. Peak memory therefore scales as O(N2), which is the property that makes it fail: doubling the grid edge quadruples the requirement.
The fix inverts the dependency. Instead of the writer pulling from a complete array, the generator pushes one block-aligned window at a time into a file already opened with a tiled profile, so only one window — O(blocksize2), a few megabytes — is ever resident. Overviews are then computed by reading back the just-written tiles at successively coarser decimation, again one window at a time, never the whole level. Peak memory becomes independent of N and depends only on the block size and the number of bands. This windowed, push-based pattern is the same tiling discipline that governs tuning Dask chunk size for spatial tiling — the block size that is efficient for the writer is usually the chunk size that is efficient for the parallel generator feeding it.
The one hard constraint: writes must be block-aligned. A window whose bounds do not fall on the internal tile grid forces the driver to read-modify-write partially filled tiles, which both reintroduces memory pressure and can corrupt the tiled layout that COG validity depends on.
The whole-array path scales peak memory with grid area and OOMs; the windowed path keeps one block resident, pushes it into a tiled file, then decimates written tiles into an overview pyramid — memory bounded by block size, not N.
import numpy as np
# The whole-array buffer for a modest 60k x 60k float32 grid:
edge, bands =60_000,1
gib = edge * edge * bands * np.dtype("float32").itemsize /2**30print(f"{gib:.1f} GiB for one live copy; build-then-write needs 2-3x")# 13.4 GiB# Even if it fits, a default strip-organized GeoTIFF is NOT a valid COG:# $ rio cogeo validate strips.tif# strips.tif is NOT a valid cloud optimized GeoTIFF.# The file is greater than 512xH or 512xW and is not tiled.# The file is greater than 512xH or 512xW and has no overviews.
A file that is not tiled and has no overviews forces a client to download the entire image to read any window — the opposite of what a COG is for.
Open the file once with a tiled COG-shaped profile, generate and write one block-aligned window at a time, then build overviews by decimation and finalize as a COG.
Set tiled=True with power-of-two blockxsize/blockysize (256 or 512), a compression the reader supports, and BIGTIFF so the file can exceed 4 GB. Use a metric CRS when the grid represents area-true density so pixel size is uniform; EPSG:4326 degree pixels distort with latitude.
python
import rasterio
from rasterio.windows import Window
BLOCK =512defopen_streaming_cog(path:str, width:int, height:int, transform, crs="EPSG:6933"):"""Open a tiled GeoTIFF sized for block-aligned windowed writes."""
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,"zlevel":6,"BIGTIFF":"YES",# required past ~4 GB}return rasterio.open(path,"w",**profile)
Iterate the internal block windows the driver reports, generate each block’s data on demand, and write it. Only one block is ever resident. Seed the generator per block from a base seed and the block index so the raster is reproducible and each block is independent — the same seeding contract the parallel grid generators use, so a windowed write and a Dask-backed write produce byte-identical output.
python
import numpy as np
import rasterio
from rasterio.windows import Window
defstream_blocks(dst, generate_block,*, base_seed:int)->None:"""Push generator output into the COG one internal block at a time."""for _, window in dst.block_windows(1):# driver's block-aligned windows# deterministic per-block stream: independent, reproducible, order-free
seed =(base_seed *1_000_003+ window.row_off *73_856_093+ window.col_off *19_349_663)&0xFFFFFFFF
rng = np.random.default_rng(seed)
block = generate_block(window, rng).astype("float32")# (h, w) for THIS windowassert block.shape ==(window.height, window.width)# block-aligned guard
dst.write(block,1, window=window)# peak memory ~ one block
After the full-resolution data is written, build decimated overviews and finalize. rio-cogeo handles the overview pass and the IFD reordering that makes the file COG-valid; it reads tiles back rather than the whole array.
python
from rasterio.enums import Resampling
from rio_cogeo.cogeo import cog_translate
from rio_cogeo.profiles import cog_profiles
deffinalize_cog(src_path:str, cog_path:str)->None:"""Add internal overviews and reorder to a valid COG, streaming tile-by-tile."""
cog_translate(
src_path, cog_path,
cog_profiles.get("deflate"),
overview_level=5,# /2 /4 /8 /16 /32
overview_resampling=Resampling.average,# area-preserving for density grids
in_memory=False,# CRITICAL: never buffer the whole raster
web_optimized=False,)
Setting in_memory=False is the load-bearing flag: with it, the overview build streams tiles from disk and peak memory stays bounded by the block size; without it, the finalize step reintroduces the very O(N2) buffer the windowed write eliminated. The overview resampling should match the semantics of the grid — average for a density surface feeding density mapping and heat generation, nearest for a categorical class raster.
Gate on two invariants: the output is a structurally valid COG, and peak RSS did not scale with grid area. Both belong in CI so a regression that reintroduces a whole-array buffer is caught before release.
python
import tracemalloc
import numpy as np
from rio_cogeo.cogeo import cog_validate
deftest_streamed_cog_is_valid_and_bounded(cog_path, block=512, bands=1):# 1. Structural COG validity: tiled + overviews + correct IFD order.
is_valid, errors, _warnings = cog_validate(cog_path)assert is_valid,f"not a valid COG: {errors}"# 2. Peak memory during a full read-back stays near one block, not the whole grid.import rasterio
tracemalloc.start()with rasterio.open(cog_path)as src:for _, window in src.block_windows(1):
_ = src.read(1, window=window)# one block resident at a time
_cur, peak = tracemalloc.get_traced_memory()
tracemalloc.stop()
block_bytes = block * block * bands * np.dtype("float32").itemsize
assert peak <8* block_bytes,f"peak {peak/2**20:.1f} MiB scales with grid, not block"
The COG-validity assertion catches a missing overview level or an untiled write; the peak-memory bound catches the subtler regression where the file is valid but some stage still buffered the full array. Together they pin both the layout and the memory contract.
A final row or column of partial blocks. When the grid edge is not a multiple of the block size, the last row and column of windows are narrower or shorter than BLOCK. Writing a full BLOCK-shaped array into a partial window raises or silently pads with garbage, so always shape each block to (window.height, window.width) as the reproducer’s guard does — never assume a square block at the margins.
Antimeridian-spanning extent breaks the transform. A grid whose geographic extent crosses ±180° longitude in EPSG:4326 produces an affine transform with a discontinuity, and windowed writes near the seam map to the wrong world coordinates so downstream range-request reads fetch the wrong tiles. Generate in an equal-area projected CRS such as EPSG:6933 (or a UTM zone) whose extent does not wrap, and only reproject on read if a client needs EPSG:4326.
Compression predictor mismatched to dtype. The horizontal-differencing predictor (predictor=2) accelerates deflate on smooth integer grids but is defined for floating point as predictor=3; using 2 on float32 data can inflate the file or, with some GDAL builds, produce tiles a strict reader rejects. Match the predictor to the dtype — 2 for integer density counts, 3 for float32 surfaces — and re-run cog_validate after any compression change.
Why must writes be aligned to the internal block grid?
A tiled GeoTIFF stores data in fixed-size tiles, and the driver can write a whole tile in one pass only when the write window matches a tile boundary. A misaligned window forces a read-modify-write of partially filled tiles, which both reintroduces memory pressure and can leave the tiled layout inconsistent, breaking COG validity. Iterating the driver's own block windows guarantees alignment for free.
What blocksize should I choose for a large synthetic grid?
Use a power of two, typically 256 or 512. Smaller tiles reduce per-window memory but multiply the number of range requests a client makes and the metadata overhead; larger tiles cut request count but raise resident memory and waste bandwidth when a reader needs only a small area. 512 is a good default for dense synthetic surfaces, and it should usually match the chunk size of any parallel generator feeding the writer.
Why is in_memory=False so important when building overviews?
The overview build is where a memory regression hides. With in_memory=False the finalizer reads full-resolution tiles from disk, decimates them, and writes each overview tile without ever holding the whole raster, so peak memory stays bounded by the block size. With in_memory=True it buffers the entire image to build the pyramid, which reintroduces the order-N-squared allocation the windowed write was designed to avoid.
Can I write the blocks out of order or in parallel?
You can write blocks in any order as long as each write is block-aligned and no two writers touch the same tile concurrently, because tiles are independent regions of the file. Seed each block deterministically from its window offset rather than a running counter so the output is identical regardless of write order, which is what lets a serial windowed write and a parallel tiled generator produce byte-for-byte the same COG.
