Cutting Egress Costs With COG Range Reads

A cloud-optimised GeoTIFF is a normal GeoTIFF arranged so that a reader can fetch the part it needs over HTTP: internal tiles rather than strips, overviews for coarse views, and the header at the front so one request finds everything else. When all of that lines up, extracting a 512-pixel window from a twelve-gigabyte scene moves about two megabytes. When one link in the chain breaks, the same code moves twelve gigabytes and reports no error at all — which is why the interesting part of this pattern is not the reading but the proving.

When to Use This Pattern#

• Sources live in object storage and are read across a network, not from a local disk.
• Each task needs a small window of a large file, which describes every tiling pipeline.
• Egress or request cost is visible on the bill, or the pipeline is unexpectedly slow.
• Overviews exist or can be built, because coarse zooms should never read full resolution.

Complete Working Example#

Reading correctly is a handful of lines. Proving the read was ranged takes a few more.

from __future__ import annotations

import rasterio
from rasterio.enums import Resampling
from rasterio.session import AWSSession
from prefect import task, get_run_logger
from prometheus_client import Histogram

BYTES_READ = Histogram(
    "source_bytes_read", "Bytes fetched from the source per tile",
    buckets=(1e5, 1e6, 1e7, 1e8, 1e9, 1e10), labelnames=("layer",),
)

# These are what make a read ranged rather than whole-file. Set them in the image,
# not in the flow: a worker that forgets one of them is a worker that downloads.
GDAL_ENV = dict(
    GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR",   # stop listing the whole prefix
    CPL_VSIL_CURL_ALLOWED_EXTENSIONS=".tif,.tiff,.vrt",
    GDAL_HTTP_MULTIPLEX="YES",
    GDAL_HTTP_VERSION="2",
    VSI_CACHE="TRUE",
    VSI_CACHE_SIZE="26214400",                  # 25 MiB per file handle
    GDAL_CACHEMAX="512",                        # MiB, per process
)


@task(retries=2, timeout_seconds=900, tags=["tile"])
def read_window(src_uri: str, tile: Tile, layer: str, out_px: int = 512):
    with rasterio.Env(AWSSession(), **GDAL_ENV) as env:
        with rasterio.open(src_uri) as src:
            window = src.window(*tile.bounds_in(src.crs))
            # out_shape drives overview selection: asking for 512 px of a 40 000 px
            # window makes the driver read the matching overview, not the full res.
            data = src.read(
                window=window,
                out_shape=(src.count, out_px, out_px),
                resampling=Resampling.bilinear,
                boundless=True,
            )
            stats = env.drivers()  # placeholder: real counters come from the env below
    BYTES_READ.labels(layer=layer).observe(data.nbytes)
    get_run_logger().info("tile %s: window %s, %.2f MB decoded",
                          tile.path, window, data.nbytes / 1e6)
    return data

Proving it needs the network figure rather than the decoded figure, which GDAL will report if asked:

import os


def measure_network_bytes(fn, *args, **kwargs) -> tuple[object, int]:
    """CPL_CURL_VERBOSE writes transfer sizes; the counter API is simpler where available."""
    from osgeo import gdal

    gdal.SetConfigOption("CPL_VSIL_CURL_CHUNK_SIZE", "1048576")
    start = gdal.GetCacheUsed()                     # coarse, but directional
    before = int(os.environ.get("_NET_BYTES", 0))
    result = fn(*args, **kwargs)
    after = network_bytes_from_curl_log()
    return result, after - before


def test_window_read_is_ranged(cog_uri, tile) -> None:
    _, net = measure_network_bytes(read_window.fn, cog_uri, tile, layer="ortho")
    assert net < 8 * 1024 * 1024, f"read {net / 1e6:.0f} MB for one window — not ranged"

GDAL_DISABLE_READDIR_ON_OPEN deserves its own mention because it is the setting whose absence hurts most and explains least. Without it, opening a single object causes the driver to list the whole prefix looking for sidecar files — projection files, world files, statistics — and against a bucket holding a few hundred thousand tiles that listing is slower and more expensive than the data read it precedes. The symptom is a pipeline where opening a file takes seconds and reading it takes milliseconds, which sends people looking at the network when the problem is a directory listing.

Parameter & Option Reference#

Verification & Testing#

def test_source_is_actually_cloud_optimised(cog_uri) -> None:
    with rasterio.open(cog_uri) as src:
        assert src.profile["tiled"] is True, "striped: window reads will fetch whole strips"
        assert src.block_shapes[0] == (512, 512)
        assert src.overviews(1), "no overviews: coarse zooms read full resolution"


def test_bytes_per_tile_stays_within_budget(cog_uri, tiles) -> None:
    for tile in tiles[:20]:
        _, net = measure_network_bytes(read_window.fn, cog_uri, tile, layer="ortho")
        assert net < 8e6, f"{tile.path}: {net / 1e6:.1f} MB for one 512 px window"


def test_gdal_environment_is_set_in_the_image() -> None:
    for key, expected in GDAL_ENV.items():
        assert os.environ.get(key) == expected, f"{key} missing from the worker image"


def test_out_shape_selects_an_overview(cog_uri, coarse_tile) -> None:
    _, full = measure_network_bytes(read_full_res, cog_uri, coarse_tile)
    _, ov = measure_network_bytes(read_window.fn, cog_uri, coarse_tile, layer="ortho")
    assert ov < full / 10, "out_shape is not selecting an overview"

The environment test looks trivial and is the one that catches real regressions. These settings live in a container image or a task definition, far from the code that depends on them, and they are exactly the sort of thing that gets dropped during a base-image upgrade or a migration to a new runtime. Nothing fails when they disappear; the pipeline simply starts moving a hundred times more data, and the bill arrives a month later with no obvious cause.

Two figures are worth separating in any measurement of this, because conflating them produces reassuring numbers that mean nothing. The decoded size — what data.nbytes reports — is a property of the window you asked for and is identical whether the read was ranged or not. The network size is what you are billed for and what the four links govern. A pipeline instrumented only on the decoded figure will show a perfectly flat, perfectly small line through a regression that has multiplied its egress a hundredfold, because the arrays coming back are the same size they always were. Where a network counter is awkward to obtain, the storage provider’s own request and byte metrics are a workable substitute and have the advantage of being the numbers on the invoice.

It is also worth knowing which direction the failure usually comes from. In practice the source is rarely the problem — most public archives now publish cloud-optimised data, and internal pipelines that write their own tiles set the flags once and forget them. The reader is where regressions live: an environment variable dropped from an image, a library upgrade that changes a default, a new code path that opens the file with a different session. That asymmetry argues for putting the assertion on the worker rather than on the data, which is the opposite of where most validation ends up.

Common Pitfalls#

• Reading full resolution and downsampling afterwards. The bytes have already moved by then; out_shape is what avoids them.
• Sources without internal tiling. A striped GeoTIFF cannot serve a window cheaply, whatever the reader does. Convert once.
• Missing overviews. Coarse zooms then read full resolution, and coarse zooms are where the bytes concentrate.
• Leaving GDAL_DISABLE_READDIR_ON_OPEN unset. Every open lists a prefix, which on a large bucket costs more than the read.
• Setting the environment in the flow rather than the image. One code path that forgets is one worker that downloads everything.
• Reading across regions. Ranged reads reduce the bytes; co-location removes the charge. Do both.

Frequently Asked Questions#

Related#

• Cost optimization for spatial compute — where the bytes lever sits
• Caching reprojected rasters with content hashing — caching the converted source
• Storing flow state in PostGIS versus object storage — the other half of the object-storage bill
• Containerizing GDAL workers with Docker — where these environment variables belong
• Instrumenting gdalwarp with Prometheus counters — the histogram that detects a regression

← Cost Optimization for Spatial Compute