Exporting Cloud-Optimized GeoTIFFs with Rasterio

A Cloud-Optimized GeoTIFF (COG) is not a new file format — it is an ordinary GeoTIFF whose internal bytes are arranged so a client can read a single map tile over HTTP with two or three range requests instead of downloading the whole raster. For a drone survey that produces multi-gigabyte orthomosaics and float32 digital elevation models, this ordering discipline is the difference between a web map that streams instantly and one that stalls a browser for thirty seconds. The trap is that a raster can be a perfectly valid GeoTIFF and still fail every COG check: the pixels must be stored in internal tiles rather than strips, decimated overviews must be present, and — the part that catches most first attempts — the overview image file directories (IFDs) and their tile data must be laid out before the full-resolution pixel data so a reader hits the low-resolution pyramid first. This page shows how to produce that layout deterministically with rasterio, which parameters actually control the on-disk structure, and how to prove a file is a real COG before you publish it.

This work is the export stage of the broader DEM/DSM generation and raster export pipeline: point clouds become rasters, and rasters become web-servable COGs. If you are weighing the pure-Python approach here against the GDAL command-line path, read GDAL vs Rasterio for DEM and orthomosaic export for the trade-offs before committing a pipeline.

Audience and prerequisites. This guide targets Python 3.10+ on a 64-bit OS with enough RAM to hold one raster window at a time — not the whole mosaic. You should understand the difference between tiled and striped TIFF storage, what a raster overview (a decimated pyramid level) is, and the distinction between an 8-bit RGB orthomosaic and a single-band float32 elevation surface, because the correct compression and predictor differ sharply between the two. rasterio wraps GDAL directly, so the on-disk result is byte-for-byte what GDAL would produce; the value here is a scriptable, testable Python surface instead of shelling out to a CLI.

# Prerequisites

Install the raster stack in an isolated environment and pin the versions below. rasterio ships GDAL inside its binary wheel, so the single most common source of “works on my machine” COG discrepancies is a build server whose rasterio wheel bundles a different GDAL than the field workstation — GDAL 3.1 introduced the dedicated COG driver, and versions before 3.4 emit subtly different overview offsets. Confirm the bundled version with rasterio.gdal_version() on both machines.

Lock the set with uv pip compile or a conda lockfile. The rio-cogeo package is optional but pays for itself: it bundles cog_validate, the same structural check that a serving stack such as a tile server or a STAC catalogue will run against your output. If your DEMs carry a compound or geoid-based vertical CRS, confirm the CRS survives the write with a round-trip read (Step 6) — the COG driver preserves it, but a hand-rolled profile copy can drop it silently.

# Conceptual architecture

Producing a COG with rasterio is a six-stage flow. You open the source raster (an orthomosaic or a DEM straight out of the PDAL DSM/DTM stage), set a creation profile that forces tiled internal blocks at a fixed block size, build decimated overviews, choose a compression codec plus a matching predictor, write the result through the COG driver so the IFDs and overview tiles are ordered correctly, and finally validate the structure. The ordering matters: overviews must exist before the file is finalized, because the COG driver reshuffles the whole file so the reduced-resolution pyramid and its tile index sit ahead of the main image data.

Each stage has a narrow contract — a raster in, a raster with more structure out — so you can unit-test overview generation independently of compression and recombine both behind one command. The sections below implement them in order.

# Step 1: Open the source raster and capture its profile

Start by opening the source with rasterio and reading the profile dictionary, which carries the pixel data type, band count, transform (the affine mapping from pixel to ground coordinates), CRS, and NoData value. Everything downstream depends on preserving these unchanged: a COG that quietly loses its transform or CRS is worse than useless because it looks valid but geolocates nothing. Read them once and treat them as the contract for the output.

import rasterio

def read_source_profile(src_path: str) -> dict:
    """Open a raster and return a copy of its profile plus useful metadata."""
    with rasterio.open(src_path) as src:
        profile = src.profile.copy()
        info = {
            "dtype": src.dtypes[0],       # e.g. 'uint8' (RGB) or 'float32' (DEM)
            "count": src.count,           # band count
            "crs": src.crs,               # coordinate reference system
            "transform": src.transform,   # affine pixel -> ground mapping
            "nodata": src.nodata,         # may be None; set it explicitly downstream
        }
    return profile, info

The dtype returned here decides your compression and predictor later: an 8-bit uint8 RGB orthomosaic and a float32 DEM want different codecs. If the source came straight from the PDAL rasterization step, its NoData will already encode the interpolation voids; if you are still seeing gaps, resolve them first with fixing holes and voids in drone-derived DSMs before exporting, because a COG faithfully preserves whatever NoData pattern it is handed.

# Step 2: Set a tiled creation profile

The first hard COG requirement is internal tiling. A default GeoTIFF stores pixels in horizontal strips, which forces a reader to fetch a full-width band of the image to get any single map tile. COGs store pixels in square internal blocks (typically 512×512) so a reader fetches exactly the block it needs. In rasterio this is tiled=True with blockxsize and blockysize set to the same power-of-two value. Set the interleave to pixel for RGB imagery so a viewer reads all bands of a tile in one contiguous run.

def make_tiled_profile(profile: dict, blocksize: int = 512) -> dict:
    """Return a profile that forces square internal tiling."""
    profile = profile.copy()
    profile.update(
        driver="GTiff",
        tiled=True,
        blockxsize=blocksize,   # must equal blockysize; power of two
        blockysize=blocksize,
        interleave="pixel",     # 'pixel' for RGB; 'band' for single-band DEM
        BIGTIFF="IF_SAFER",     # auto-switch to BigTIFF past ~4 GB
    )
    return profile

BIGTIFF="IF_SAFER" is the safety valve for large surveys: a classic TIFF caps at 4 GB of addressable offsets, and a full-resolution city orthomosaic plus its overviews routinely exceeds that. IF_SAFER lets GDAL promote the container to BigTIFF automatically when the projected size gets close, avoiding a mid-write TIFFAppendToStrip: Maximum TIFF file size exceeded abort.

# Step 3: Build decimated overviews

The second COG requirement is overviews — a pyramid of progressively downsampled copies of the raster embedded in the same file. Without them, a client zoomed out to a whole-survey view would have to read every full-resolution pixel and downsample on the fly. You build overviews at decimation factors (each level halves resolution) and pick a resampling method appropriate to the data: average for continuous surfaces like a DEM or a photographic orthomosaic, nearest for categorical or palettized rasters where inventing intermediate values would be wrong.

from rasterio.enums import Resampling

def add_overviews(path: str, factors=(2, 4, 8, 16), method="average") -> None:
    """Build embedded decimated overviews in place."""
    resampling = Resampling[method]      # 'average', 'nearest', 'bilinear', ...
    with rasterio.open(path, "r+") as dst:
        dst.build_overviews(factors, resampling)
        dst.update_tags(ns="rio_overview", resampling=method)

Choose the number of levels so the coarsest overview is roughly one tile wide — for a 20000-pixel-wide mosaic, factors (2, 4, 8, 16, 32) bring the top level down to about 625 pixels, near a single 512 block. The setting compression and overviews for orthomosaic COGs page works through how far to take the pyramid and why average versus nearest changes both file size and visual fidelity.

# Step 4: Choose compression and a matching predictor

Compression is where an orthomosaic and a DEM diverge most. DEFLATE and ZSTD are lossless and safe for both; ZSTD compresses faster and usually smaller but requires a GDAL built with ZSTD support. The multiplier is the predictor: predictor 2 (horizontal differencing) helps integer rasters where adjacent pixels are similar, and predictor 3 (floating-point differencing) is specifically for float DEMs — using predictor 2 on float32 data produces larger files or corrupt output. LZW is the universally supported fallback. For pure RGB visual orthomosaics where a small amount of loss is acceptable, JPEG compression at quality 85 shrinks files dramatically, but never apply it to a DEM.

def compression_options(dtype: str, lossy_rgb: bool = False) -> dict:
    """Pick a codec + predictor from the raster dtype."""
    if dtype.startswith("float"):
        return {"compress": "ZSTD", "predictor": 3, "zstd_level": 9}
    if lossy_rgb and dtype == "uint8":
        return {"compress": "JPEG", "quality": 85, "photometric": "YCBCR"}
    # integer, lossless default
    return {"compress": "DEFLATE", "predictor": 2, "zlevel": 6}

Predictor 3 on floating-point elevation is not optional cosmetics — it commonly halves DEM size because neighbouring terrain heights differ by centimetres, and differencing turns those into long runs of near-zero bytes that DEFLATE and ZSTD crush. The full decision table across RGB, float, and palettized inputs lives in setting compression and overviews for orthomosaic COGs.

# Step 5: Write through the COG driver

The third and least obvious COG requirement is IFD and tile ordering: the overview directories and their tile data must sit at the front of the file, ahead of the full-resolution image, so a reader fetches the header and pyramid index in the first range request. You cannot get this ordering by writing a normal GTiff — you must write through GDAL’s dedicated COG driver, which reorganizes the whole file on close. The COG driver builds overviews for you if you pass overview_resampling, so the cleanest path is to open the source, set the COG creation options, and copy the data across in one pass.

import rasterio
from rasterio.shutil import copy as rio_copy

def write_cog(src_path: str, dst_path: str, *, blocksize: int = 512,
              compress: str = "DEFLATE", predictor: int = 2,
              overview_resampling: str = "average") -> None:
    """Emit a valid COG via the GDAL COG driver (correct IFD ordering)."""
    cog_profile = {
        "driver": "COG",
        "blocksize": blocksize,               # COG driver uses 'blocksize', not block[xy]size
        "compress": compress,
        "predictor": predictor,
        "overview_resampling": overview_resampling,
        "BIGTIFF": "IF_SAFER",
    }
    # rio_copy streams windowed data; it never loads the whole raster into RAM.
    rio_copy(src_path, dst_path, **cog_profile)

Two ordering notes that trip people up: the COG driver takes blocksize (one key), not the blockxsize/blockysize pair the GTiff driver uses, and it manages overviews internally — if you pre-built overviews in Step 3 on a GTiff and then re-emit through the COG driver, the driver rebuilds them from the source resolution, so you do not double up. When a file that “should” be a COG fails validation anyway, the specific reasons and fixes are catalogued in fixing COG validation errors in GDAL.

# Step 6: Validate the output structure

A file that opens in QGIS is not necessarily a COG. Validation asserts three things: the pixels are internally tiled, decimated overviews are present, and the byte layout passes the cloud-optimized structural check. rasterio exposes tiling and overview state directly, and rio-cogeo’s cog_validate runs the same layout test a serving stack uses.

import rasterio
from rio_cogeo.cogeo import cog_validate

def verify_cog(path: str) -> None:
    with rasterio.open(path) as ds:
        assert ds.profile["tiled"] is True, "raster is striped, not tiled"
        # overviews() returns the decimation factors present for band 1
        factors = ds.overviews(1)
        assert len(factors) >= 3, f"too few overviews: {factors}"
        assert ds.crs is not None, "CRS was dropped on write"
        assert ds.transform is not rasterio.Affine.identity(), "no geotransform"
    is_valid, errors, warnings = cog_validate(path)
    assert is_valid, f"not a valid COG: {errors}"
    print(f"valid COG · overviews={factors} · warnings={warnings}")

The overviews(1) call returns the decimation factors GDAL actually wrote (e.g. [2, 4, 8, 16]); an empty list is the single most common validation failure and means the overviews never made it into the file. The cog_validate return splits errors (structural failures that make the file not a COG) from warnings (block-size or overview-count advisories that still serve correctly). Treat errors as a hard gate in CI and log warnings for review.

# Parameter deep-dive

Every option trades file size against read performance or write time. Tune them against a representative tile, not the whole survey.

# Verification and output inspection

Beyond the structural gate in Step 6, prove that the pixels survived the round trip. Read a window from the source and the same window from the COG and assert they match within the codec’s tolerance — exact for lossless DEFLATE/ZSTD, approximate for JPEG. Also assert the NoData value carried across, because a dropped NoData turns transparent voids into opaque black in a web viewer.

import numpy as np
import rasterio
from rasterio.windows import Window

def assert_roundtrip(src_path: str, cog_path: str, lossless: bool = True) -> None:
    win = Window(0, 0, 512, 512)
    with rasterio.open(src_path) as a, rasterio.open(cog_path) as b:
        assert a.nodata == b.nodata, f"nodata drift: {a.nodata} -> {b.nodata}"
        sa = a.read(1, window=win)
        sb = b.read(1, window=win)
        if lossless:
            assert np.array_equal(sa, sb), "lossless pixels changed on write"
        else:
            # JPEG: assert visually close, not identical
            assert np.abs(sa.astype("int32") - sb).mean() < 3.0, "excess JPEG loss"
    print("round-trip verified")

Run this on a tile that contains both data and NoData so the mask path is exercised. In CI, wrap the whole export in a context manager and call the validation as a non-zero exit gate so a malformed COG never reaches a tile server. When the export itself blows up on memory rather than structure — a “cannot allocate memory” abort on a very large mosaic — the diagnosis and windowed-write fix live in decoding GDAL “cannot allocate memory” on export.

# Troubleshooting

cog_validate reports “The file is greater than 512xN or Nx512, it is not tiled”. Your raster is stored in strips, not internal tiles. You wrote it with a plain GTiff driver without tiled=True, or an intermediate step re-striped it. Re-emit through driver="COG" (Step 5), which always produces tiled output, and re-run the validator.

The COG opens fine locally but tiles load slowly over HTTP. Overviews are missing, so a zoomed-out request reads full-resolution pixels. Check ds.overviews(1) — if it returns an empty list, no pyramid was written. Rebuild through the COG driver with overview_resampling set, and confirm the coarsest level is about one block wide.

A float32 DEM COG is barely smaller than the uncompressed source. You applied predictor 2 (or none) to floating-point data. Horizontal integer differencing does nothing useful on IEEE floats. Switch to predictor=3, the floating-point predictor, and DEM size typically drops by roughly half because neighbouring elevations differ by only centimetres.

Transparent void areas render as solid black in the web map. The NoData value was not preserved on write, so the viewer has no mask to make voids transparent. Set nodata explicitly in the creation profile and assert a.nodata == b.nodata after writing (Verification section); for RGB, consider an internal alpha mask band instead.

TIFFAppendToStrip: Maximum TIFF file size exceeded on a large city mosaic. The output crossed the 4 GB classic-TIFF ceiling. Add BIGTIFF="IF_SAFER" to the creation options so GDAL promotes the container to BigTIFF automatically once the projected size approaches the limit.

CPLE_NotSupported: driver COG does not support creation option PREDICTOR (or ZSTD). The bundled GDAL is too old for that option, or was built without ZSTD. Check rasterio.gdal_version(); you need GDAL ≥ 3.4 for the full COG option set, and a ZSTD-enabled build to use compress="ZSTD". Fall back to DEFLATE if you cannot upgrade.

# Related

• GDAL vs Rasterio for DEM and Orthomosaic Export
• Fixing COG Validation Errors in GDAL
• Setting Compression and Overviews for Orthomosaic COGs
• Generating DSM and DTM from Point Clouds with PDAL

← DEM/DSM Generation & Raster Export Automation