Windowed COG Reads over HTTP with Rasterio

One-sentence answer: set GDAL’s HTTP options once per process, reproject the field boundary into the raster’s CRS, turn its bounds into a pixel window rounded outward, and read that window — two to four range requests instead of a 100 MB download.

Context

A Sentinel-2 band asset is roughly 60–120 MB and covers 110 km square. A field covers a few hundred metres. Downloading the file to read 0.02% of it is the default behaviour of most code, and it is why so many imagery pipelines are described as “running overnight”. The cloud-optimized layout exists precisely so a reader can fetch a small header, learn where the internal tiles are, and request only the byte ranges it needs — but only if the client is configured to take advantage of it.

This guide is the read half of satellite imagery APIs and STAC catalogues, applied to the scenes accepted in searching Sentinel-2 scenes with pystac-client. The same mechanics apply to your own published imagery — see cloud-optimized storage for field imagery.

What a windowed read asks the object store for A sequence from the reader to the object store. First a small header range request of about sixteen kilobytes returns the tile offsets. Then three range requests fetch the internal tiles that overlap the field window, about three hundred kilobytes each. A comparison line shows the whole-file alternative at over one hundred megabytes. rasterio + GDAL /vsicurl/ holds a 64 MB per-file block cache object store B08.tif · 104 MB 512 px internal tiles 5 overview levels GET Range: bytes=0-16383 — header and tile index 16 KB — offsets of every tile 3 × GET Range — only the tiles the window overlaps ≈ 900 KB total, 60–200 ms whole-file alternative: 104 MB, 8–20 s, and the disk to put it on Without GDAL_DISABLE_READDIR_ON_OPEN, an extra prefix listing precedes every one of these opens.

Prerequisites

Beyond the parent topic’s stack: rasterio 1.3.* built against GDAL 3.6+, and network access to the object store. Read the GDAL settings below before the first rasterio.open in the process — some are only consulted at open time.

Step-by-step

1. Configure GDAL in a context manager so the settings are explicit and scoped.

From a boundary in degrees to a window in pixels Five steps: start from the field boundary in EPSG:4326, reproject it into the raster's own coordinate reference system, convert its bounds to a fractional pixel window, round that window outward with a small pad so boundary pixels are not clipped away, and read it in two to four HTTP range requests. Boundary EPSG:4326 polygon Reproject into the raster CRS a few coordinates Bounds to window fractional offsets Round outward plus a 2 px pad Read 2–4 range requests Reproject the polygon, never the remote raster: one is a handful of coordinates, the other resamples every pixel.

2. Reproject the boundary into the raster’s CRS.

3. Compute the window from bounds and round outward, so no field pixel is clipped away.

4. Read each band onto one grid, resampling the 20 m bands explicitly rather than relying on broadcasting.

5. Time the read and assert it.

PYTHON 
import time

import geopandas as gpd
import numpy as np
import rasterio
from rasterio.enums import Resampling
from rasterio.windows import from_bounds, Window
from shapely.geometry import shape

GDAL_HTTP = {
    "GDAL_DISABLE_READDIR_ON_OPEN": "EMPTY_DIR",
    "CPL_VSIL_CURL_ALLOWED_EXTENSIONS": ".tif,.TIF,.tiff",
    "GDAL_HTTP_MULTIPLEX": "YES",
    "GDAL_HTTP_VERSION": "2",
    "VSI_CACHE": "TRUE",
    "VSI_CACHE_SIZE": str(64 * 1024 * 1024),
    "GDAL_CACHEMAX": "512",
}


def read_field_window(href: str, boundary_geojson: dict, *, pad_px: int = 2):
    """Read one band's field window from a remote COG. Returns (array, transform, crs)."""
    with rasterio.Env(**GDAL_HTTP):
        with rasterio.open(href) as src:
            assert src.crs is not None, f"{href} has no CRS"
            proj = (gpd.GeoSeries([shape(boundary_geojson)], crs="EPSG:4326")
                    .to_crs(src.crs).iloc[0])
            assert proj.intersects(
                shape({"type": "Polygon", "coordinates": [[
                    (src.bounds.left, src.bounds.bottom), (src.bounds.right, src.bounds.bottom),
                    (src.bounds.right, src.bounds.top), (src.bounds.left, src.bounds.top),
                    (src.bounds.left, src.bounds.bottom)]]})), (
                "the field does not intersect this raster — wrong tile or wrong CRS")

            win = from_bounds(*proj.bounds, transform=src.transform)
            # Round outward and pad, so boundary pixels are never clipped off.
            win = Window(
                col_off=int(np.floor(win.col_off)) - pad_px,
                row_off=int(np.floor(win.row_off)) - pad_px,
                width=int(np.ceil(win.width)) + 2 * pad_px,
                height=int(np.ceil(win.height)) + 2 * pad_px,
            ).intersection(Window(0, 0, src.width, src.height))

            arr = src.read(1, window=win)
            transform = src.window_transform(win)
            return arr, transform, src.crs


def read_aligned_stack(hrefs: dict[str, str], boundary_geojson: dict, target_band: str = "B08"):
    """Read several bands onto the target band's grid, resampling coarser ones explicitly."""
    ref, ref_transform, ref_crs = read_field_window(hrefs[target_band], boundary_geojson)
    out = {target_band: ref}
    with rasterio.Env(**GDAL_HTTP):
        for band, href in hrefs.items():
            if band == target_band:
                continue
            with rasterio.open(href) as src:
                proj = (gpd.GeoSeries([shape(boundary_geojson)], crs="EPSG:4326")
                        .to_crs(src.crs).iloc[0])
                win = from_bounds(*proj.bounds, transform=src.transform)
                arr = src.read(1, window=win, out_shape=ref.shape,
                               resampling=Resampling.bilinear)
            out[band] = arr
    shapes = {b: a.shape for b, a in out.items()}
    assert len(set(shapes.values())) == 1, f"bands are not aligned after resampling: {shapes}"
    return out, ref_transform, ref_crs

Inline verification — assert the latency, because correctness alone will not reveal a striped file:

PYTHON 
start = time.perf_counter()
arr, transform, crs = read_field_window(item.assets["nir"].href, boundary)
elapsed = time.perf_counter() - start

print(f"{arr.shape} pixels in {elapsed * 1000:.0f} ms, CRS {crs}")
assert arr.size > 0, "empty window — the field does not overlap this raster"
assert elapsed < 2.0, (
    f"window read took {elapsed:.1f}s — the file is probably striped, "
    "or GDAL_DISABLE_READDIR_ON_OPEN is unset")

Settings worth being deliberate about

Setting Default here Why it matters
GDAL_DISABLE_READDIR_ON_OPEN EMPTY_DIR The single highest-impact setting. Without it, opening one asset lists the entire bucket prefix; in a prefix with thousands of objects that listing dominates the request
VSI_CACHE_SIZE 64 MB Per-file block cache. Reading many windows from one open dataset reuses cached blocks; reading one window per open never touches it
GDAL_CACHEMAX 512 MB Process-wide block cache. Raising it helps repeated reads over the same tile and does nothing for a single pass over many scenes
GDAL_HTTP_MULTIPLEX / VERSION YES / 2 HTTP/2 multiplexing lets several range requests share one connection, which is most of the benefit of concurrency for small reads
pad_px 2 Pixels of margin around the window. Enough to survive rounding at the boundary; a larger pad reads whole extra tiles for nothing
Concurrent reads 16 I/O bound, so more threads help until the object store’s per-prefix limit answers with 503s that retries then amplify

Set these once at process start-up, not per call. They are read when a dataset is opened, so a worker that configures them inside a loop pays for the first open in the loop having none of them.

Gotchas and edge cases

  • Rounding a window inward silently trims the field. from_bounds returns fractional offsets; truncating them drops the partial pixels at the boundary, which is where mixed edge pixels live. Round outward and pad, then mask afterwards — the reverse order loses data you cannot recover.
Why the rounding direction decides what the edge of the field means Two panels contrasting inward and outward rounding of a pixel window. Rounding inward truncates the partial pixels at the boundary, which are the mixed canopy and soil pixels, and the loss is invisible. Rounding outward with a pad keeps them, and masking afterwards marks what falls outside. Rounding the window inward from_bounds returns fractional offsets Truncating drops the partial edge pixels. Those are exactly the mixed canopy-and-soil pixels. The loss is invisible — the array is simply smaller. Boundary statistics are quietly biased. Rounding outward, then masking Floor the offsets, ceil the size, add a pad Every field pixel is inside the array. Clip to the polygon afterwards, not before. Nodata marks what falls outside the boundary. Nothing is lost that cannot be recovered.

  • Band resolutions differ within one scene. Sentinel-2 red and NIR are 10 m; red edge and the classification band are 20 m. Reading both without an explicit out_shape produces arrays of different sizes, and any arithmetic on them either raises or, worse, broadcasts into something meaningless. Resample once, deliberately, as above.

  • rasterio.Env settings are read at open time. Setting them after opening a dataset changes nothing for that dataset. Wrap the open, not just the read.

  • Bilinear resampling of a categorical band is wrong. The classification band holds class codes; interpolating them produces class 6.5, which is not a class. Use Resampling.nearest for anything categorical — the same rule as in clipping rasters to field boundaries.

  • Cross-region reads cost egress on every range request. They work perfectly and quietly generate a bill. Run the pipeline in the archive’s region, or copy the windows you need once and read them locally afterwards.

  • The per-file cache is per open. Reading the same asset for twenty fields opens it twenty times unless you keep the dataset open and loop the windows inside. For a field group sharing a tile, open once and read many.

  • A signed URL that expires mid-run looks like a corrupt file. Assets behind time-limited credentials return a 403 partway through a long read, and GDAL commonly surfaces that as an unreadable dataset rather than as an authorisation error. Refresh the signature before each batch rather than once at the start of a job, and when a read fails on a file that worked an hour ago, check the credential’s lifetime before suspecting the object.

This guide is part of Satellite Imagery APIs & STAC Catalogues for Field Monitoring — see there for search, reflectance scaling and manifest registration.

Parent section Satellite Imagery APIs & STAC Catalogues for Field Monitoring Related Searching Sentinel-2 Scenes with pystac-client