Benchmarking COG Read Throughput from Object Storage

A useful raster benchmark measures three things per configuration — time to open, time to read a window, and bytes transferred — repeated enough times to see past the noise:

import time
import rasterio
from rasterio.windows import Window

with rasterio.Env(GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR"):
    t0 = time.perf_counter()
    with rasterio.open(url) as src:
        t_open = time.perf_counter() - t0
        t1 = time.perf_counter()
        arr = src.read(1, window=Window(2048, 2048, 512, 512))
        t_read = time.perf_counter() - t1

print(f"open {t_open*1000:.0f} ms · read {t_read*1000:.0f} ms · {arr.nbytes/1e6:.1f} MB decoded")

Measurement is what turns the guidance in Optimizing Pipeline Cost and Performance into decisions for your data, your region and your machine.

Why This Arises in Remote Sensing Workflows

Advice about COG reads is full of numbers — tile sizes, ingest sizes, thread counts — and all of them are contingent. They depend on the file’s internal layout, the distance between compute and storage, the object size distribution, and the shape of the windows the pipeline actually requests. A setting that halves the runtime of one pipeline can do nothing for another.

Benchmarking resolves that, and it is cheap: the whole exercise is a few hundred reads against one representative file. The difficulty is not running it but running it so the result means something. Object storage has a long latency tail, several caches warm up during the first requests, and a benchmark that changes two variables at once produces a number nobody can act on.

The payoff is concrete. Knowing that your reads cost 30 milliseconds and four requests each, rather than assuming it, tells you whether a job is I/O-bound before you scale the cluster — and it tells you which of the settings in Caching PROJ Data and GDAL Config in Containers are worth setting.

Where reads stop being latency-bound For small windows the time is almost constant, because it is dominated by request latency rather than by bytes. Past roughly one megabyte of useful data the curve becomes linear as bandwidth takes over. The knee is where batching stops helping and parallelism starts to. Median read time versus window size, same region 400 ms 0 latency-dominated time ≈ constant; more requests, not more bytes bandwidth-dominated time ∝ bytes; read less, or read in parallel the knee — around 1 MB here 64 KB 64 MB Left of the knee, concurrency is the lever. Right of it, transferring fewer bytes is. The knee moves with region, instance type and object size — which is exactly why it is measured, not assumed.

Environment & Setup

Package Version Why
rasterio ≥1.3.0 The reads being measured
numpy ≥1.23 Statistics over repeats
pandas ≥2.0 Tabulating the sweep
psutil (optional) ≥5.9 Process-level counters alongside timings 
pip install "rasterio>=1.3.0" "pandas>=2.0" "numpy>=1.23"

Complete Working Example

A benchmark harness that measures one configuration at a time, repeats it, and reports medians with spread.

Benchmark hygiene rules and what each prevents Turning caching off prevents a benchmark that measures the cache. Discarding the warm-up prevents measuring TLS and credential setup. Varying one axis prevents unattributable results, and recording the configuration is what makes two benchmarks comparable at all. Four rules that make numbers mean something rule prevents cost VSI_CACHE=FALSE while benchmarking measuring the cache, not the network none discard the first repeat measuring TLS and credential setup one read vary one axis at a time results with no attributable cause more runs record region, image, GDAL version benchmarks that cannot be compared a few fields report median and p90 a single sample from a heavy tail none A benchmark whose configuration is not recorded cannot be compared with the next one. 
import statistics
import time
from concurrent.futures import ThreadPoolExecutor
from dataclasses import asdict, dataclass

import pandas as pd
import rasterio
from rasterio.windows import Window

BASE_ENV = {
    "GDAL_DISABLE_READDIR_ON_OPEN": "EMPTY_DIR",
    "GDAL_HTTP_MULTIRANGE": "YES",
    "GDAL_HTTP_MERGE_CONSECUTIVE_RANGES": "YES",
    "VSI_CACHE": "FALSE",          # OFF for benchmarking: caching hides the real cost
}


@dataclass
class Result:
    label: str
    window_px: int
    concurrency: int
    open_ms: float
    read_ms: float
    total_ms: float
    mb_decoded: float


def one_read(url: str, size: int, offset: tuple[int, int] = (2048, 2048)) -> Result:
    """A single open-and-read, timed in two parts."""
    t0 = time.perf_counter()
    with rasterio.open(url) as src:
        t_open = time.perf_counter() - t0
        t1 = time.perf_counter()
        arr = src.read(1, window=Window(offset[0], offset[1], size, size))
        t_read = time.perf_counter() - t1
    return Result("", size, 1, t_open * 1e3, t_read * 1e3,
                  (t_open + t_read) * 1e3, arr.nbytes / 1e6)


def sweep_window_size(url: str, sizes=(256, 512, 1024, 2048, 4096), repeats: int = 11,
                      env: dict | None = None, label: str = "default") -> pd.DataFrame:
    """One axis at a time: window size, everything else fixed."""
    rows = []
    with rasterio.Env(**(env or BASE_ENV)):
        for size in sizes:
            timings = []
            for i in range(repeats):
                r = one_read(url, size)
                if i == 0:
                    continue                 # discard the warm-up
                timings.append(r)
            rows.append({
                "label": label,
                "window_px": size,
                "open_ms_median": round(statistics.median(t.open_ms for t in timings), 1),
                "read_ms_median": round(statistics.median(t.read_ms for t in timings), 1),
                "read_ms_p90": round(sorted(t.read_ms for t in timings)[int(0.9 * len(timings))], 1),
                "mb": round(timings[0].mb_decoded, 2),
                "mb_per_s": round(timings[0].mb_decoded /
                                  (statistics.median(t.read_ms for t in timings) / 1e3), 1),
            })
    return pd.DataFrame(rows)


def sweep_concurrency(url: str, size: int = 512, levels=(1, 2, 4, 8, 16, 32),
                      reads_per_level: int = 64) -> pd.DataFrame:
    """How much parallelism the storage backend actually gives you."""
    rows = []
    with rasterio.Env(**BASE_ENV):
        for n in levels:
            offsets = [(2048 + (i % 8) * 512, 2048 + (i // 8) * 512) for i in range(reads_per_level)]
            t0 = time.perf_counter()
            with ThreadPoolExecutor(max_workers=n) as pool:
                list(pool.map(lambda off: one_read(url, size, off), offsets))
            elapsed = time.perf_counter() - t0
            rows.append({"concurrency": n,
                         "reads": reads_per_level,
                         "elapsed_s": round(elapsed, 2),
                         "reads_per_s": round(reads_per_level / elapsed, 1)})
    return pd.DataFrame(rows)


if __name__ == "__main__":
    url = "https://example-bucket.s3.amazonaws.com/scene_cog.tif"

    print(sweep_window_size(url).to_string(index=False))
    print()
    print(sweep_concurrency(url).to_string(index=False))

    # Compare one setting change against the baseline, everything else fixed
    tuned = BASE_ENV | {"GDAL_INGESTED_BYTES_AT_OPEN": "65536"}
    print(sweep_window_size(url, env=tuned, label="ingest-64k").to_string(index=False))

VSI_CACHE=FALSE during benchmarking is deliberate and worth explaining: with caching on, the second read of the same window is free, which makes every configuration look excellent and tells you nothing. Turn it back on in production, where repeated access is real.

Variant Patterns

1. Measuring requests, not just time

Time conflates latency, bandwidth and request count. Counting requests separates them.

import rasterio

with rasterio.Env(CPL_CURL_VERBOSE="YES", CPL_DEBUG="ON"):
    with rasterio.open(url) as src:
        _ = src.profile
        arr = src.read(1, window=Window(2048, 2048, 512, 512))
# GDAL writes one line per HTTP request to stderr; count them with `2>&1 | grep -c '^> GET'`

Requests per useful megabyte is the metric that predicts cost most directly, because object storage bills for both. A configuration with slightly worse latency but half the requests is usually the cheaper one at scale — the accounting laid out in Reducing S3 Egress Costs in Raster Pipelines.

2. What to sweep, and in what order

Sweep the big levers first Region placement typically changes read time by a factor of three to five, tiling and alignment by two to four, GDAL open settings by up to two, and thread count by a factor of two before the backend throttles. Sweeping in that order avoids tuning a small lever while a large one is set wrong. Typical effect size, largest first region placement 3–5× — and it is the only one that also changes the bill per GB tiling and window alignment 2–4× — a property of the file, not the reader GDAL open settings 1.5–2× thread count up to 2×, then throttling Tuning threads while running cross-region is the classic misallocation of effort. Fix the top row first; the lower rows only matter once it is right.

3. Benchmarking the pipeline, not just the read

A read benchmark answers “how fast can I get bytes”. It does not answer “is my pipeline I/O-bound”, which needs the arithmetic timed alongside.

import time

t0 = time.perf_counter()
arr = src.read(1, window=win, masked=True)
t_read = time.perf_counter() - t0

t1 = time.perf_counter()
ndvi = (nir - arr) / (nir + arr)
t_compute = time.perf_counter() - t1

print(f"I/O {t_read:.3f}s · compute {t_compute:.3f}s · ratio {t_read / max(t_compute, 1e-9):.1f}")

A ratio well above one says the job is I/O-bound and that concurrency or region placement is the lever. Below one says the arithmetic dominates, and no amount of read tuning will help — the distinction that decides between the executors in Scaling Raster Processing with Dask.

Making the Numbers Trustworthy

Four rules keep a benchmark from misleading the team that acts on it.

Change one thing at a time. A run that changes both the ingest size and the thread count produces a number with no attributable cause, and the next person will change both back.

Discard the warm-up and report the spread. The first read of a session pays for DNS, TLS and credential resolution; the median of ten subsequent reads, with a p90 alongside, describes what a job will actually experience.

Benchmark where the job runs. A measurement from a laptop over a home connection is a measurement of that connection. Run the harness inside the same container image and region as production — the reproducibility argument from Containerizing Geospatial Python Environments.

Record the configuration with the result. Region, instance type, image tag, GDAL version and the environment variables in force all belong in the output table. A benchmark whose configuration is not recorded cannot be compared with the next one, which is the only use a benchmark really has.

Common Errors

Every configuration looks identical

Caching is on, so only the first read touches the network. Set VSI_CACHE=FALSE for the benchmark and vary the window offset between repeats.

Results vary by a factor of three between runs

Too few repeats against a heavily skewed latency distribution. Use at least ten repeats per configuration and compare medians, not single timings.

Throughput collapses above a certain concurrency

You have found the backend’s throttling point, which is a useful result. Set the pipeline’s concurrency below it rather than adding retries — the cap-then-back-off pattern from Parameterizing Prefect Flows for Multi-Tile Runs.

Frequently Asked Questions

Q: Why are my benchmark numbers so variable? Object storage latency is heavily skewed, and caches at several layers warm up. Run each configuration at least ten times, discard the first, and report the median with an interquartile range rather than a single number.

Q: Should I benchmark from my laptop? Only to compare relative effects. Absolute numbers from a laptop are dominated by the route to the bucket, so a setting that looks decisive locally may be irrelevant in-region, and vice versa.

Q: What is the single most useful metric? Requests per useful megabyte. It captures both the per-request overhead and the wasted transfer, and it predicts cost on object storage better than either latency or throughput alone.