Windowed Reads from Cloud Optimized GeoTIFF

A Cloud Optimized GeoTIFF lets you read a 512×512 patch from a 40 GB raster on object storage without downloading the file. This guide shows the windowed-read pattern with Rasterio, including reading by geographic coordinates rather than raw pixel offsets, and the PyProj reprojection that keeps your query window landing where you intend. It is for anyone working with satellite imagery, orthophotos, or DEMs too large to localize. It sits under Cloud-Native Geospatial Formats in Geospatial Data Ingestion & Processing Workflows.

Why This Approach / What Goes Wrong

A COG is internally tiled with overviews, and its TIFF header maps each tile to a byte range through the TileOffsets/TileByteCounts arrays. Rasterio — via GDAL's /vsicurl/ virtual filesystem — reads a Window by issuing HTTP range requests for just the tiles that overlap your area of interest: seconds and megabytes instead of minutes and gigabytes. The whole point is to never localize the file.

Three things defeat this, and they account for almost every "my windowed read is still slow/wrong" report.

There is a fourth, subtler cost that only shows up once the first three are fixed: the transfer is quantised by the internal block size, not by your window. A tile is the smallest unit GDAL can decompress, so a 300×300-pixel window over a 512×512-tiled COG still fetches and decodes four full tiles — roughly 1 MB of pixels to hand you 90,000. The read is correct, just less efficient than the numbers suggest, and the ratio gets worse the more small scattered windows you issue. Below the block size, asking for less costs the same; the lever is fetching fewer, larger, block-aligned windows rather than many tiny ones.

Layered on top of that is GDAL's own chunking. /vsicurl/ does not request exactly the byte ranges the driver asks for; it fetches in fixed chunks (CPL_VSIL_CURL_CHUNK_SIZE, 16 KB by default) and caches them. A 400 KB tile therefore becomes dozens of small HTTP requests unless the chunk size is raised or consecutive ranges are merged. On a high-latency link that request count, not the byte count, is what you are actually waiting for.

Data flow of a Cloud Optimized GeoTIFF windowed read An area-of-interest bounding box in EPSG:4326 is reprojected with PyProj into the raster CRS, then from_bounds maps it through the affine transform to a pixel Window. Rasterio via GDAL vsicurl issues HTTP range requests for only the internal tiles overlapping the Window, transferring megabytes from a Cloud Optimized GeoTIFF. A striped GeoTIFF has no tile table, so the whole file must stream sequentially — gigabytes. AOI bounding box EPSG:4326 · lon, lat PyProj reproject Transformer → raster CRS from_bounds() affine → pixel Window Rasterio via GDAL /vsicurl HTTP range requests for tiles overlapping the Window COG — internally tiled AOI Only overlapping tiles fetched · megabytes Striped GeoTIFF — contrast No tile table — whole file streams · gigabytes
The AOI is reprojected and mapped to a pixel Window; on a COG, Rasterio range-requests only the overlapping tiles, whereas a striped GeoTIFF forces a sequential read of the entire file.

Prerequisites

conda install -c conda-forge "rasterio=1.3.*" "pyproj=3.4.*"

Install from conda-forge rather than pip so GDAL, PROJ, and the Rasterio binding stay ABI-compatible — the same environment discipline used across Raster Data Handling with Rasterio.

Step-by-Step Implementation

1. Configure GDAL for efficient remote reads, then open the COG. The two options below stop GDAL from listing the whole bucket and from probing sidecar files it will never find over HTTP — both otherwise add a round trip per open.

import rasterio
from rasterio.env import Env

gdal_opts = {
    "GDAL_DISABLE_READDIR_ON_OPEN": "EMPTY_DIR",   # don't list the bucket on open
    "CPL_VSIL_CURL_ALLOWED_EXTENSIONS": ".tif",    # skip probing for .aux/.ovr sidecars
    "GDAL_HTTP_MULTIRANGE": "YES",                  # coalesce tile fetches into one request
    "GDAL_HTTP_MERGE_CONSECUTIVE_RANGES": "YES",    # adjacent tiles become one range
    "CPL_VSIL_CURL_CHUNK_SIZE": "1048576",          # 1 MB chunks instead of the 16 KB default
    "GDAL_CACHEMAX": "512",                         # MB of decoded block cache per process
    "GDAL_HTTP_MAX_RETRY": "3",                     # ride out transient 503s
    "GDAL_HTTP_RETRY_DELAY": "1",
}
cog_url = "https://example-bucket.s3.amazonaws.com/sentinel_ortho_cog.tif"

with Env(**gdal_opts):
    with rasterio.open(cog_url) as src:
        print(src.profile["driver"], src.shape, src.crs)  # GTiff (20000, 20000) EPSG:32633
        print(src.block_shapes[0])                        # (512, 512) — internally tiled
        print(src.overviews(1))                           # [2, 4, 8, 16] — pyramid present

If block_shapes reports a full-width strip like (1, 20000) or overviews(1) is empty, the file is not a COG — jump to the debugging section before going further.

Two of those options are worth understanding rather than pasting. CPL_VSIL_CURL_CHUNK_SIZE sets the granularity at which /vsicurl/ fetches and caches; the 16 KB default is tuned for reading a header, not for pulling 400 KB tiles, so raising it to 1 MB typically cuts the request count for a single window read by an order of magnitude. Push it much higher and you start over-fetching neighbours you never asked for, which hurts when windows are sparse. GDAL_CACHEMAX is the decoded-block cache, in megabytes per process — raise it when you read many overlapping windows from the same file, leave it alone when each window is visited once, and remember that it is per process, so sixteen worker processes at 512 MB will try to use 8 GB.

HTTP round trips exchanged during one windowed read A request sequence between Rasterio driving GDAL vsicurl and an S3 bucket. The first round trip is a ranged GET of the opening sixteen kilobytes, which returns the TIFF header carrying the TileOffsets and TileByteCounts arrays. The bucket listing and the probes for aux and ovr sidecar files are suppressed by GDAL_DISABLE_READDIR_ON_OPEN and CPL_VSIL_CURL_ALLOWED_EXTENSIONS, so they are never sent. The second round trip is a multi-range GET that coalesces the four overlapping tile byte ranges into a single request and returns roughly three megabytes, leaving the rest of the forty gigabyte file in the bucket. What actually crosses the wire on one windowed read Rasterio · GDAL /vsicurl S3 object storage GET Range: bytes=0-16383 TIFF header — TileOffsets + TileByteCounts · 16 KB bucket listing and .aux.xml / .ovr probes — never sent GDAL_DISABLE_READDIR_ON_OPEN · CPL_VSIL_CURL_ALLOWED_EXTENSIONS GET Range: bytes=a-b, c-d, e-f, g-h four tile ranges coalesced by GDAL_HTTP_MULTIRANGE=YES 4 × 512 × 512 tiles ≈ 3 MB — the AOI patch 2 round trips · roughly 3 MB transferred the remaining 39.99 GB never leaves the bucket
The three GDAL options collapse the exchange to two round trips: one ranged GET for the tile offset table, one multi-range GET for the tiles the window touches.

2. Convert a geographic bounding box to a pixel window. Reproject the AOI bounds into the raster's CRS first, then let from_bounds derive the pixel Window from the affine transform. Pass always_xy=True so PyProj keeps lon/lat (x/y) order and does not silently swap the axes — the classic PROJ 6+ axis-order trap.

from rasterio.windows import from_bounds
from pyproj import Transformer

# Area of interest in EPSG:4326 (lon, lat): a slice of central Berlin
aoi_4326 = (13.35, 52.50, 13.42, 52.54)   # (min_lon, min_lat, max_lon, max_lat)

with Env(**gdal_opts), rasterio.open(cog_url) as src:
    # Reproject AOI corners into the raster CRS (here EPSG:32633, UTM 33N — a metric grid)
    to_raster = Transformer.from_crs("EPSG:4326", src.crs, always_xy=True)
    xmin, ymin = to_raster.transform(aoi_4326[0], aoi_4326[1])
    xmax, ymax = to_raster.transform(aoi_4326[2], aoi_4326[3])

    window = from_bounds(xmin, ymin, xmax, ymax, transform=src.transform)
    aoi_patch = src.read(1, window=window)          # only overlapping tiles transfer
    patch_transform = src.window_transform(window)  # georeferencing for the patch

3. Read a decimated overview when full resolution is more than you need. For a wide-area preview, pass out_shape so GDAL pulls a coarser pyramid level and moves a fraction of the bytes — invaluable for thumbnails or when the AOI spans thousands of pixels.

with Env(**gdal_opts), rasterio.open(cog_url) as src:
    # Ask for the same window at ~1/8 resolution — reads an overview, not full-res tiles
    preview = src.read(
        1,
        window=window,
        out_shape=(int(window.height // 8), int(window.width // 8)),
    )
    print(preview.shape)   # e.g. (18, 39) — a cheap decimated view of the AOI
Overview levels and the bytes each one moves Five stacked bars represent the pyramid levels stored inside a Cloud Optimized GeoTIFF for the same area of interest. Full resolution reads a 4096 by 4096 pixel window and moves about 33.5 megabytes; the 2x, 4x, 8x and 16x decimated levels read 2048, 1024, 512 and 256 pixel windows and move 8.4 megabytes, 2.1 megabytes, 524 kilobytes and 131 kilobytes respectively. The 8x level is highlighted as the one an out_shape of window height and width divided by eight selects. Which pyramid level answers the read overview level window read at that level bytes moved 1:1 full res 4096 × 4096 px ≈ 33.5 MB 2× decimated 2048 × 2048 px ≈ 8.4 MB 4× decimated 1024 × 1024 px ≈ 2.1 MB 8× decimated 512 × 512 px out_shape=(h // 8, w // 8) one range request against a smaller level ≈ 524 KB 16× decimated 256 × 256 px ≈ 131 KB The pyramid ships inside the COG — a coarser out_shape costs one range request, never a full-resolution read plus resample.
Asking for a smaller out_shape makes GDAL serve the request from a pre-built pyramid level, so an eight-fold decimation moves roughly one sixty-fourth of the bytes.

4. Snap the window to whole internal blocks so no tile is decoded twice. from_bounds returns a floating-point window that almost never aligns with the 512-pixel grid. Rounding it out to full blocks costs you a few extra pixels and saves partial-tile handling, which matters when you loop over many adjacent windows: unaligned windows make neighbouring reads re-fetch the same boundary tiles.

from rasterio.windows import from_bounds, round_window_to_full_blocks, intersection, Window

with Env(**gdal_opts), rasterio.open(cog_url) as src:
    raw = from_bounds(xmin, ymin, xmax, ymax, transform=src.transform)

    # Expand outward to the internal tile grid, then clip to the raster extent
    aligned = round_window_to_full_blocks(raw, src.block_shapes)
    aligned = intersection(aligned, Window(0, 0, src.width, src.height))

    print(raw)      # Window(col_off=8134.2, row_off=3901.7, width=312.4, height=148.9)
    print(aligned)  # Window(col_off=7680.0, row_off=3584.0, width=1024.0, height=1024.0)
    block_patch = src.read(1, window=aligned)

The intersection call is not optional. round_window_to_full_blocks happily returns offsets past the raster edge for an area of interest that overhangs the image, and read on an out-of-range window raises WindowError rather than clipping for you.

5. Read past the edge deliberately with a boundless read. When the area of interest genuinely straddles the raster boundary — a tile-aligned processing grid, or an AOI drawn before anyone checked the footprint — boundless=True returns an array of exactly the requested shape, padding the outside with fill_value. Without it you get a smaller array than you asked for, and every downstream index is off by the difference.

import numpy as np

with Env(**gdal_opts), rasterio.open(cog_url) as src:
    overhanging = Window(col_off=src.width - 200, row_off=0, width=512, height=512)
    padded = src.read(1, window=overhanging, boundless=True, fill_value=0)

    print(padded.shape)                    # (512, 512) — full shape, not (512, 200)
    print(int((padded == 0).sum()))        # 159744 — the padded region outside the raster

Mask the fill before computing statistics; a fill_value of 0 is a legitimate reflectance or elevation value, so treat the padded cells as nodata explicitly with np.ma.masked_equal or by reading src.read_masks.

6. Persist the patch as its own small GeoTIFF, carrying the correct transform. Write patch_transform, not the source transform, or the patch is georeferenced to the parent raster's origin.

with Env(**gdal_opts), rasterio.open(cog_url) as src:
    profile = src.profile.copy()

profile.update(
    width=aoi_patch.shape[1],
    height=aoi_patch.shape[0],
    transform=patch_transform,
    count=1,
)

with rasterio.open("aoi_patch.tif", "w", **profile) as dst:
    dst.write(aoi_patch, 1)

Multi-band imagery follows the identical pattern with src.read(window=window) returning a (bands, rows, cols) stack — see reading multi-band TIFFs with Rasterio for band-indexing details.

7. Fan out over many windows or many scenes with a thread pool. The dominant cost of a remote windowed read is network latency, not CPU, and GDAL releases the GIL while waiting on HTTP — so threads scale this workload well. The one hard rule is that a dataset handle is not thread-safe: open a fresh handle inside each worker rather than sharing one across the pool.

from concurrent.futures import ThreadPoolExecutor

sample_windows = [Window(c, r, 512, 512) for r in range(0, 4096, 512)
                                          for c in range(0, 4096, 512)]

def read_window(win):
    # One dataset handle per call — sharing a handle across threads corrupts reads
    with Env(**gdal_opts), rasterio.open(cog_url) as src:
        return win, src.read(1, window=win)

with ThreadPoolExecutor(max_workers=8) as pool:
    patches = dict(pool.map(read_window, sample_windows))

print(len(patches), "windows read concurrently")   # 64 windows read concurrently

Eight to sixteen workers is the useful range against a single object store. Beyond that the bucket begins returning 503 SlowDown, and because each worker carries its own GDAL block cache, memory grows linearly with the pool. If the windows all come from the same file, reopening per call is cheap once the header is in the /vsicurl/ cache; if they come from hundreds of different scenes, the per-open header fetch dominates and you are better off batching by scene.

When the windows feed a larger tabular workflow rather than a per-patch computation, the same fan-out belongs at the task level instead — see Scaling with Dask-GeoPandas for the partitioned equivalent.

Verification

Confirm the window is the expected size and the patch is georeferenced where you asked — bracketing the reprojected AOI corners from step 2.

import rasterio

print("Patch shape:", aoi_patch.shape)          # Patch shape: (148, 312)
assert aoi_patch.size > 0, "Empty window — bounds may miss the raster extent or wrong CRS"

with rasterio.open("aoi_patch.tif") as check:
    left, bottom, right, top = check.bounds
    print("Patch bounds (raster CRS):", round(left), round(bottom), round(right), round(top))
    assert check.crs == rasterio.crs.CRS.from_epsg(32633)
    # left/bottom/right/top should bracket (xmin, ymin, xmax, ymax) from step 2

The stronger check is on bytes, not shape: a read that returns the right pixels having transferred the whole file is still a failure. Turn on GDAL's curl tracing for one read and count the requests.

import rasterio

with rasterio.Env(CPL_CURL_VERBOSE="YES", CPL_DEBUG="ON", **gdal_opts):
    with rasterio.open(cog_url) as src:
        _ = src.read(1, window=aligned)
# stderr shows one ranged GET for the header and a small number of tile ranges.
# A single unbounded "GET /sentinel_ortho_cog.tif" with no Range header means
# the file is striped, or the server ignored Range and returned 200 instead of 206.

A 200 OK where you expected 206 Partial Content is the diagnostic for a server or CDN that has not been configured to honour range requests — the client is correct, the storage layer is not.

Edge Cases & Debugging

Frequently Asked Questions

Is it ever faster to just download the whole COG? Yes, once your windows cover a large enough fraction of the file or you revisit it many times. The crossover is roughly when total windowed transfer approaches the file size: a few hundred scattered 512-pixel patches from a 40 GB scene is overwhelmingly a win for windowed reads, while a hundred windows that between them touch most tiles of a 200 MB scene is not. Latency counts too — a hundred sequential windows at 80 ms round trip each is eight seconds of pure waiting, which a single bulk download may beat.

Should I read in the raster's native CRS or reproject on the fly? Read natively and reproject the small result. WarpedVRT can present the dataset in any CRS and lets you window in target coordinates, which is convenient, but the warp is recomputed on every read and it resamples pixels you may then resample again downstream. Native reads keep the pixel values untouched and defer exactly one resampling step to the end. Reserve on-the-fly warping for the case where many differently-projected sources must be read on a single common grid.

How do I pick the right overview level instead of hardcoding a divisor? Derive it from the ratio between the window size and the pixels you actually need to display, then let out_shape do the selection: GDAL picks the closest level at or above your requested resolution. Hardcoding // 8 breaks when a scene ships a different overview factor list — check src.overviews(1) and note that a COG may have been written with factors [2, 4, 8, 16, 32] or with none at all, a choice made at encode time as covered in resampling and overviews when writing COGs.

Does this work against a plain web server, or do I need S3? Any HTTP server that supports Range requests and reports Accept-Ranges: bytes works — /vsicurl/ cares about the protocol, not the vendor. What breaks it is a CDN or proxy configured to buffer and re-serve whole objects, which silently converts every ranged request into a full download. Verify with a curl -I for Accept-Ranges and a ranged curl -r 0-1023 that returns 206.

Why is my first window read slow and the rest fast? The first read pays for the header and tile-offset table, which on a large multi-band COG can be several megabytes of TileOffsets/TileByteCounts, plus TLS handshake and DNS. Subsequent reads hit the /vsicurl/ cache for the header and reuse the connection. If every read is slow, the dataset handle is being reopened per read inside a loop — hoist the rasterio.open outside the loop when reads are sequential in one thread.

Can I compute area or distance from a windowed patch directly? Only if the raster is in a suitable projected CRS. Pixel dimensions come from the affine transform, so a patch from a geographic (degree-based) raster has cells whose ground size varies with latitude, and a Web Mercator raster has cells whose ground size varies even more. Reproject the patch — or better, read from a source already in a UTM zone or national grid, chosen the way Choosing a UTM Zone Automatically in Python describes — before turning pixel counts into square metres.

Related