Writing Cloud-Optimized GeoTIFFs from Interpolated Sensor Grids

A Cloud-Optimized GeoTIFF (COG) is a standard GeoTIFF whose internal bytes are arranged so an HTTP client can read just the pixels it needs over a byte-range request, without downloading the whole file. For interpolated environmental surfaces served to web maps, tile caches, or partner pipelines, that layout is the difference between a 200 ms window read and a multi-megabyte full-file fetch. The recipe is specific: internal tiling at a fixed block size, a pyramid of decimated overviews, appropriate compression, and — critically — a header and tile ordering the COG driver produces but a hand-rolled GeoTIFF does not. Starting from a georeferenced DataArray produced in the gridded raster fields with xarray and rioxarray workflow, you can write and validate a COG in Python with rioxarray and rasterio in one pass.

What Makes a GeoTIFF “Cloud-Optimized”

Three properties, all internal to the file, distinguish a COG from an ordinary GeoTIFF. Internal tiling stores pixels in square blocks (conventionally 512×512) rather than in scanline strips, so a reader fetching a small map window touches only the overlapping tiles. Internal overviews are progressively decimated copies of the raster — half resolution, quarter, eighth — so a zoomed-out view reads a tiny coarse pyramid level instead of decimating the full-resolution grid on the fly. Header-first layout places the metadata and the image file directory (IFD) entries at the front of the file and orders overviews before full-resolution data, so a client learns the tile offsets from the first few kilobytes and issues precise range requests for everything after.

For interpolated sensor grids these properties matter more than they do for satellite imagery, because the surfaces are smooth and compress extremely well: a 400×400 PM2.5 field can be a few hundred kilobytes, and the overhead of building overviews is trivial while the payoff at read time is large. The one place smoothness bites is overview resampling — average resampling preserves the field magnitude when zoomed out, whereas nearest (correct for classified rasters) would make a continuous field look blocky at coarse zoom.

Production-Ready Implementation

The function below takes a CRS-aware DataArray — one that already carries a CRS, a nodata value, and a valid north-up transform — and writes a valid COG. It uses rioxarray’s to_raster with driver="COG", which builds the overviews and arranges the layout in a single call, and it validates the result before returning so a broken write fails loudly rather than silently shipping a non-COG file.

# python 3.11 · numpy==1.26.4 · xarray==2024.3.0 · rioxarray==0.15.0 · rasterio==1.3.9 · pyproj==3.6.1
from __future__ import annotations
from pathlib import Path
import xarray as xr
import rioxarray  # noqa: F401 — registers the .rio accessor
import rasterio
from rasterio.enums import Resampling


def write_cog_from_grid(
    da: xr.DataArray,
    path: str | Path,
    *,
    blocksize: int = 512,
    compress: str = "DEFLATE",
    predictor: int = 2,
    overview_resampling: str = "average",
    overview_levels: tuple[int, ...] = (2, 4, 8, 16),
) -> Path:
    """
    Write a CRS-aware DataArray to a valid Cloud-Optimized GeoTIFF.

    Parameters
    ----------
    da : xr.DataArray
        Georeferenced 2-D field with dims (y, x). Must already have
        ``da.rio.crs`` and ``da.rio.nodata`` set and a north-up transform.
    path : str | Path
        Destination ``.tif`` path.
    blocksize : int
        Internal tile size in pixels. 512 is the standard COG block; use 256
        only for small rasters. Must be a multiple of 16.
    compress : str
        'DEFLATE' (with predictor=2) for continuous float fields, 'LZW' for
        classified integer rasters. Never lossy JPEG on measurement data.
    predictor : int
        2 = horizontal differencing, which improves compression on smooth
        surfaces. Set to 1 for classified/integer rasters.
    overview_resampling : str
        'average' for continuous concentration/temperature fields, 'nearest'
        for integer class rasters. Averaging a sentinel nodata corrupts values.
    overview_levels : tuple[int, ...]
        Decimation factors for the overview pyramid.

    Returns
    -------
    Path
        The written COG path (validated before return).

    Raises
    ------
    ValueError
        If the input lacks a CRS/transform or the written file fails COG checks.

    Time complexity:  O(N) in pixels for the base write, plus ~O(N/3) for the
                      full overview pyramid (geometric series over 4x decimation).
    Space complexity: O(blocksize^2) streamed per tile, not the whole grid.
    """
    path = Path(path)
    if da.rio.crs is None:
        raise ValueError("DataArray has no CRS; call rio.write_crs first.")
    if da.rio.transform().e >= 0:
        raise ValueError("Transform is not north-up (e >= 0); sort y descending.")

    da.rio.to_raster(
        path,
        driver="COG",
        blocksize=blocksize,
        compress=compress,
        predictor=predictor,
        overview_resampling=Resampling[overview_resampling],
        overviews="AUTO",
        overview_levels=list(overview_levels),
        num_threads="ALL_CPUS",
        BIGTIFF="IF_SAFER",
    )

    _assert_valid_cog(path, blocksize)
    return path


def _assert_valid_cog(path: Path, blocksize: int) -> None:
    """Reopen with rasterio and confirm the COG-defining properties hold."""
    with rasterio.open(path) as src:
        prof = src.profile
        if not prof.get("tiled", False):
            raise ValueError(f"{path} is not internally tiled.")
        if prof.get("blockxsize") != blocksize:
            raise ValueError(
                f"{path} block size {prof.get('blockxsize')} != {blocksize}."
            )
        # band 1 must carry an overview pyramid for zoomed-out reads
        if not src.overviews(1):
            raise ValueError(f"{path} has no internal overviews.")
        if src.crs is None:
            raise ValueError(f"{path} lost its CRS on write.")

The _assert_valid_cog gate is deliberately strict: a tiled GeoTIFF is not automatically a COG, so checking tiled, blockxsize, and the presence of overviews catches the most common ways a “COG” write silently degrades to a plain GeoTIFF.

Parameter Tuning Guide

Block size, compression, and overview depth scale with grid size and the read pattern. The defaults suit a typical city-scale interpolated field; adjust as the raster grows or the consumer changes.

Grid size (pixels) Block size Compression Predictor Overview levels Notes
< 256 × 256 256 DEFLATE 2 (2, 4) Small field; a 512 tile would cover most of it
300 × 300 (dense urban) 512 DEFLATE 2 (2, 4, 8) Standard air-quality mesh output
400 × 400 (regional) 512 DEFLATE 2 (2, 4, 8, 16) Deeper pyramid for wider zoom range
Integer class raster 512 LZW 1 (2, 4, 8) — nearest Averaging would invent fake class codes
> 2000 × 2000 (continental) 512 DEFLATE 2 (2, 4, 8, 16, 32) BIGTIFF=IF_SAFER; overviews essential

Overview depth rule of thumb: add levels until the coarsest overview is roughly a single block — around 512 pixels on its longest side. A 4000-pixel raster needs decimation to /8 (500 px), so (2, 4, 8) at minimum; adding /16 and /32 costs little and smooths very-zoomed-out reads.

Compression on smooth fields: DEFLATE with predictor=2 typically beats LZW by 15–30% on continuous interpolated surfaces because horizontal differencing turns smooth gradients into small residuals. For integer class rasters the predictor hurts, so pair LZW with predictor=1.

Verification and Testing

Validate structurally with rasterio, not by eye. The test below writes a COG from a synthetic georeferenced field and asserts the three COG-defining properties — tiling, overviews, and a preserved CRS — plus round-trip fidelity of the finite values.

# python 3.11 · numpy==1.26.4 · xarray==2024.3.0 · rioxarray==0.15.0 · rasterio==1.3.9
import numpy as np
import xarray as xr
import rioxarray  # noqa: F401
import rasterio
import pytest


def _synthetic_field(n: int = 400) -> xr.DataArray:
    x = 500000 + (np.arange(n) + 0.5) * 50.0          # UTM easting, 50 m cells
    y = 5200000 - (np.arange(n) + 0.5) * 50.0          # descending: north first
    vals = (20 + 5 * np.sin(np.linspace(0, 4, n))[:, None]
               + 5 * np.cos(np.linspace(0, 4, n))[None, :]).astype("float32")
    da = xr.DataArray(vals, dims=("y", "x"), coords={"y": y, "x": x}, name="pm25")
    da = da.rio.write_crs("EPSG:32633")
    da = da.rio.write_nodata(np.nan, encoded=True)
    return da


def test_cog_is_valid(tmp_path):
    da = _synthetic_field()
    out = write_cog_from_grid(da, tmp_path / "field.tif")

    with rasterio.open(out) as src:
        assert src.profile["tiled"] is True
        assert src.profile["blockxsize"] == 512
        assert src.overviews(1)                       # non-empty overview list
        assert src.crs.to_epsg() == 32633
        # round-trip: base resolution values survive within float32 tolerance
        read_back = src.read(1)
        assert np.allclose(read_back, da.values, atol=1e-4, equal_nan=True)


def test_missing_crs_rejected(tmp_path):
    da = _synthetic_field().rio.write_crs(None)
    with pytest.raises(ValueError, match="no CRS"):
        write_cog_from_grid(da, tmp_path / "bad.tif")

If you have the GDAL command-line tools available, rio cogeo validate field.tif from the rio-cogeo plugin is the authoritative external check; the rasterio-only assertions above cover the same properties without an extra dependency and run inside your test suite.

Gotchas

A tiled GTiff is not a COG. Writing with driver="GTiff" plus tiled=True produces internal tiles but neither internal overviews nor the header-first IFD ordering the spec requires. Validators reject it, and clients fall back to full-file reads. Always use driver="COG", which arranges the layout correctly.

average overviews blend sentinel nodata. If nodata is an integer sentinel like -9999, average resampling mixes it into neighboring real values, painting a halo of wrong numbers around masked regions at coarse zoom. Use np.nan nodata (which resampling skips) for float fields, or nearest resampling if you must keep an integer sentinel.

Overbuilding overviews on tiny rasters wastes space. A 200×200 field does not need a five-level pyramid; the overviews would be larger overhead than benefit. Match overview depth to grid size using the rule of thumb above — stop when the coarsest level is about one block.

Re-tiling an existing plain GeoTIFF loses the interpolation attributes. If you convert a previously written GeoTIFF to COG with a fresh rasterio write, carry the source tags() forward explicitly; the provenance attributes (interpolation method, station count, timestamp) do not migrate automatically and an unattributed raster is unauditable.