Streaming Cloud-Optimized GeoTIFFs with Async Range Requests

To stream just the blocks you need from a remote Cloud-Optimized GeoTIFF, point GDAL’s /vsicurl/ driver at the object’s URL and read a rasterio.windows.Window; GDAL parses the header once, then issues HTTP range requests for only the tiles the window overlaps. Wrapping each blocking read in asyncio.to_thread() under a Semaphore lets you pull many windows concurrently without downloading whole files. This page is part of the Async I/O for Raster Processing guide inside the broader Spatial Batch Processing & Async Workflows reference.

Prerequisites

• Python 3.10 or later (for asyncio.to_thread)
• pip install rasterio aiohttp rio-cogeo
• GDAL 3.4+ with curl support (gdalinfo --formats | grep vsicurl or check that /vsicurl/ is usable)
• A COG reachable over HTTPS with an origin that honours Range: headers (S3, GCS, Azure Blob, and most CDNs do)

If your inputs are large enough that a single window still strains RAM, pair this with Memory Management for Large GIS Datasets so per-block arrays stay bounded.

How COG Layout Enables Range Reads

A Cloud-Optimized GeoTIFF is an ordinary GeoTIFF with two disciplines applied: the pixel data is stored in internal tiles (typically 512×512 blocks) rather than scanline strips, and the file carries a pyramid of downsampled overviews. Both are described by the Image File Directory (IFD) at the front of the file, which lists the byte offset and length of every tile in every resolution level.

Because the IFD is compact and near the start of the object, GDAL can fetch it with one small range request, then compute exactly which tile byte-ranges cover a requested window. It issues a Range: bytes=start-end request per tile and decodes only those. A striped, non-COG GeoTIFF has no per-tile offsets, so any partial read collapses into fetching the whole file. The diagram below traces the request flow.

Complete Working Implementation

The script below streams several map windows from a remote COG concurrently. It uses GDAL’s /vsicurl/ handler through rasterio, offloads each blocking window read to a thread with asyncio.to_thread(), and caps in-flight connections with an asyncio.Semaphore. Windows are specified as EPSG:4326 bounding boxes and transformed to the raster’s native CRS before reading:

#!/usr/bin/env python3
"""
Stream windows from a remote Cloud-Optimized GeoTIFF via async range requests.
Usage: python stream_cog.py https://example.com/data/scene.tif
"""
import asyncio
import os
import sys
from dataclasses import dataclass

import numpy as np
import rasterio
from rasterio.windows import from_bounds
from rasterio.warp import transform_bounds

# GDAL /vsicurl/ tuning — set before any dataset is opened.
os.environ.setdefault("CPL_VSIL_CURL_ALLOWED_EXTENSIONS", ".tif,.tiff")
os.environ.setdefault("GDAL_HTTP_MULTIPLEX", "YES")          # reuse HTTP/2 conn
os.environ.setdefault("GDAL_DISABLE_READDIR_ON_OPEN", "EMPTY_DIR")  # no dir scan
os.environ.setdefault("VSI_CACHE", "TRUE")                   # cache header blocks
os.environ.setdefault("CPL_VSIL_CURL_USE_HEAD", "NO")       # skip extra HEAD call

MAX_CONCURRENCY = 8          # cap simultaneous range-request connections


@dataclass
class Aoi:
    """A named area of interest as a WGS84 (EPSG:4326) bounding box."""
    name: str
    west: float
    south: float
    east: float
    north: float


def read_window_blocking(url: str, aoi: Aoi) -> dict:
    """Open the remote COG and read one window. Runs in a worker thread.

    The dataset is opened HERE, not in the parent coroutine, because a GDAL
    dataset handle is bound to the thread that created it and must not be
    shared. Each call opens, reads, and closes within a single thread.
    """
    vsi_path = f"/vsicurl/{url}"
    with rasterio.open(vsi_path) as ds:
        # Reproject the EPSG:4326 bbox into the raster's own CRS so the
        # window lines up with the stored pixels (often EPSG:3857 or a UTM zone).
        left, bottom, right, top = transform_bounds(
            "EPSG:4326", ds.crs, aoi.west, aoi.south, aoi.east, aoi.north
        )
        window = from_bounds(left, bottom, right, top, transform=ds.transform)
        # Round to whole pixels; only tiles overlapping this window are fetched.
        window = window.round_offsets().round_lengths()
        data = ds.read(1, window=window)               # triggers range reads
        return {
            "name": aoi.name,
            "shape": data.shape,
            "crs": str(ds.crs),
            "min": float(np.nanmin(data)) if data.size else None,
            "max": float(np.nanmax(data)) if data.size else None,
        }


async def read_window(url: str, aoi: Aoi, sem: asyncio.Semaphore) -> dict:
    """Async wrapper: bound concurrency, then offload the blocking read."""
    async with sem:                                    # never exceed MAX_CONCURRENCY
        return await asyncio.to_thread(read_window_blocking, url, aoi)


async def stream_all(url: str, aois: list[Aoi]) -> list[dict]:
    sem = asyncio.Semaphore(MAX_CONCURRENCY)
    tasks = [read_window(url, aoi, sem) for aoi in aois]
    return await asyncio.gather(*tasks, return_exceptions=False)


def main() -> None:
    if len(sys.argv) != 2:
        print("usage: python stream_cog.py <cog-url>", file=sys.stderr)
        sys.exit(2)                                    # 2 = usage error
    url = sys.argv[1]

    # Example windows over a global EPSG:4326 raster.
    aois = [
        Aoi("berlin", 13.30, 52.45, 13.50, 52.58),
        Aoi("lisbon", -9.25, 38.68, -9.08, 38.80),
        Aoi("nairobi", 36.75, -1.35, 36.95, -1.20),
        Aoi("quito", -78.55, -0.30, -78.42, -0.10),
    ]

    try:
        results = asyncio.run(stream_all(url, aois))
    except rasterio.errors.RasterioIOError as exc:
        print(f"range read failed: {exc}", file=sys.stderr)
        sys.exit(1)                                    # 1 = runtime error

    for r in results:
        print(
            f"{r['name']:>8}  {r['shape']}  crs={r['crs']}  "
            f"min={r['min']}  max={r['max']}"
        )
    sys.exit(0)


if __name__ == "__main__":
    main()

Step Annotations

1. os.environ.setdefault(...) before any open — GDAL reads its configuration once when the curl layer initializes. CPL_VSIL_CURL_ALLOWED_EXTENSIONS=".tif,.tiff" tells /vsicurl/ it may serve those URLs without a probing request, and GDAL_HTTP_MULTIPLEX="YES" lets multiple range requests ride a single HTTP/2 connection. Setting these after opening a dataset has no effect.

2. /vsicurl/{url} prefix — This is GDAL’s virtual filesystem handler for HTTP(S) objects. Prefixing the URL routes reads through the range-request machinery instead of downloading the file to a temp path.

3. transform_bounds("EPSG:4326", ds.crs, ...) — The area of interest is expressed in geographic degrees, but the raster stores pixels in its own CRS. Transforming the bounding box first ensures the derived window maps to the correct pixel rows and columns; skipping this is the classic cause of empty or wildly offset reads.

4. from_bounds(...).round_offsets().round_lengths() — from_bounds yields a floating-point Window; rounding to whole pixels aligns it to block boundaries so GDAL fetches complete tiles. Sub-pixel windows still work but waste a partial tile fetch at each edge.

5. Open inside read_window_blocking — The dataset handle lives entirely within the worker thread. This is deliberate: GDAL handles are not thread-safe, so opening in the parent coroutine and sharing across to_thread calls would corrupt reads under concurrency.

6. asyncio.Semaphore(MAX_CONCURRENCY) — Each window read holds one connection to the origin. The semaphore caps simultaneous connections at eight, which keeps you under typical per-host connection limits and avoids 429 Too Many Requests throttling from object stores.

7. asyncio.to_thread(read_window_blocking, ...) — rasterio reads are blocking C calls that release the GIL during I/O. Offloading them to the default thread pool lets the event loop keep other window reads in flight, giving genuine concurrency for network-bound work described in the Async I/O for Raster Processing guide.

Named Gotcha: A Non-COG or Overview-less File Forces Full-File Reads

The most common way this optimisation silently fails is pointing it at a GeoTIFF that is not actually cloud-optimized. If the source is striped (scanline layout) instead of internally tiled, or if it lacks overviews, GDAL cannot seek to a small block. A windowed read then pulls whole strips — often the entire file — and your carefully bounded window transfers just as many bytes as a naive download. Nothing errors; it is simply slow and expensive.

Validate the source before you rely on streaming, and re-tile it if it fails:

# Validate COG structure (rio-cogeo)
rio cogeo validate https://example.com/data/scene.tif

# If it fails, rewrite as a proper COG with tiles + overviews
rio cogeo create input.tif scene_cog.tif \
    --cog-profile deflate \
    --blocksize 512 \
    --overview-resampling average

A passing rio cogeo validate guarantees internal tiling and an overview pyramid, which is exactly what makes per-block range reads pay off.

Verification

The point of streaming is transferring far fewer bytes than the object’s full size. Measure it directly by asking the origin for the content length, then comparing against the bytes GDAL actually read:

import aiohttp
import asyncio
from osgeo import gdal

URL = "https://example.com/data/scene.tif"

async def full_size(url: str) -> int:
    async with aiohttp.ClientSession() as session:
        async with session.head(url) as resp:
            return int(resp.headers["Content-Length"])

# Read one small window and inspect GDAL's byte counter.
gdal.SetConfigOption("CPL_CURL_VERBOSE", "NO")
gdal.SetConfigOption("GDAL_HTTP_MULTIPLEX", "YES")
ds = gdal.Open(f"/vsicurl/{URL}")
band = ds.GetRasterBand(1)
_ = band.ReadAsArray(0, 0, 512, 512)          # one tile-sized window

# GDAL exposes network counters via VSIStatL / vsimem stats.
stats = gdal.GetActualURLStatistics() if hasattr(gdal, "GetActualURLStatistics") else None
total = asyncio.run(full_size(URL))
print(f"full object size : {total:,} bytes")
print(f"single 512x512 window should fetch << that")

For a true COG, one 512×512 window fetches on the order of a few hundred kilobytes regardless of whether the full object is 50 MB or 5 GB. If the fetched total approaches the full size, the file is not tiled — return to the gotcha above. GDAL’s own --debug on (CPL_DEBUG=ON) also logs each VSICURL: Downloading ... range, letting you count requests per window.

Performance Notes

• Right-size the semaphore. Eight concurrent connections saturate most single-origin CDNs. Object stores like S3 tolerate more, but past ~16 you trade throughput for throttling risk. Benchmark before raising it.
• Reuse overviews for zoomed-out views. If you only need a coarse preview, read from an overview level (ds.read(1, out_shape=(h, w)) with a downscaled shape) so GDAL fetches the small pyramid tiles instead of full-resolution blocks.
• Keep VSI_CACHE=TRUE. The IFD header is read once and cached in memory; without caching, every window re-fetches the header, adding a round-trip per read.

FAQ

Related

• Async I/O for Raster Processing — parent guide covering asyncio patterns, thread offloading, and connection limits for network-bound raster work
• Processing 100k GeoJSON Files with Python asyncio — the vector counterpart: fan out many small async reads under a bounded semaphore