Forestry GIS Python · Forestry · Ecology

Reading Cloud Optimized GeoTIFFs with rasterio Windows

A windowed read is the operation that makes cloud-native geospatial work: given a URL and a study area, fetch the bytes covering the study area and nothing else. It is three lines of rasterio when it goes right and a full scene transfer when it goes wrong, and the difference is entirely in how the window is derived and how GDAL is configured. This guide covers that narrow step inside Cloud-Native Raster Workflows for Forestry, part of Ecological GIS Data Foundations in Python.

When to use this approach

Read pattern Bytes moved Use when
src.read() on a remote COG The whole scene Never, unless you genuinely need every pixel
src.read(window=...) The overlapping internal tiles A study area smaller than the scene — the normal case
src.read(out_shape=...) The matching overview level A preview, a quicklook, or a coarse regional summary
rasterio.mask.mask(..., crop=True) The geometry’s bounding window You want the clip and the mask in one step
Download, then read locally The whole scene, once You will read the same scene many times and disk is cheaper than repeat requests

Minimal reproducible example

The core is rasterio.windows.from_bounds, which converts map coordinates into a pixel window. Everything difficult is on either side of it: the bounds have to be in the raster’s CRS, and the window has to be snapped so the read lands on whole internal tiles.

import geopandas as gpd
import rasterio
from rasterio.windows import from_bounds

GDAL_ENV = {
    "GDAL_DISABLE_READDIR_ON_OPEN": "EMPTY_DIR",
    "CPL_VSIL_CURL_ALLOWED_EXTENSIONS": ".tif,.TIF,.tiff",
    "GDAL_CACHEMAX": "512",
    "VSI_CACHE": "TRUE",
}


def read_stand_window(cog_url: str, stands_path: str, band: int = 1):
    """Read only the pixels covering a stand layer from a remote COG."""
    stands = gpd.read_file(stands_path)
    with rasterio.Env(**GDAL_ENV), rasterio.open(cog_url) as src:
        # The bounds MUST be expressed in the raster's own CRS, not the vector's.
        minx, miny, maxx, maxy = stands.to_crs(src.crs).total_bounds
        window = from_bounds(minx, miny, maxx, maxy, transform=src.transform)
        window = window.round_offsets().round_lengths()
        # Clip to the raster so a partly-overlapping stand layer does not error.
        window = window.intersection(
            rasterio.windows.Window(0, 0, src.width, src.height))
        data = src.read(band, window=window)
        transform = src.window_transform(window)
        profile = src.profile | {
            "height": data.shape[0], "width": data.shape[1],
            "transform": transform, "count": 1,
        }
    return data, profile

The to_crs(src.crs) call is the line that most often goes missing. Passing latitude and longitude bounds to a raster in UTM does not raise — from_bounds happily computes a window somewhere near the origin, the read succeeds, and you get a small array of NoData that looks like a genuinely empty stand.

Snapping a request window to internal tile boundaries A grid of internal COG tiles is drawn with a study-area rectangle crossing several of them at arbitrary offsets. The requested window is shown as a solid outline that cuts through tiles; the enlarged window that snaps outward to whole tile edges is shown as a shaded region. The shading beyond the requested rectangle is labelled as pixels the server sends anyway, because a tile is the smallest unit that can be transferred. requested shaded = whole tiles the server must send Tile-aligned read Why alignment matters • a tile is the smallest transferable unit • a window off by one pixel touches an extra row and column of tiles • at 512 px tiles that is a large fraction of a small request The fix round the window outward, then slice the array you wanted from what you received

Read at the right resolution

For a regional overview you rarely need full resolution. Passing out_shape makes GDAL serve the request from the nearest overview level, which is both faster and far smaller — a 10× decimation reads roughly 1% of the bytes.

def read_overview(cog_url: str, bounds, decimation: int = 8):
    """Read a decimated window, served from the COG's overview pyramid.

    `bounds` must already be expressed in the raster's own CRS.
    """
    with rasterio.Env(**GDAL_ENV), rasterio.open(cog_url) as src:
        window = from_bounds(*bounds, transform=src.transform).round_offsets().round_lengths()
        out_h = max(1, int(window.height) // decimation)
        out_w = max(1, int(window.width) // decimation)
        data = src.read(
            1, window=window, out_shape=(out_h, out_w),
            resampling=rasterio.enums.Resampling.average,
        )
        transform = src.window_transform(window) * rasterio.Affine.scale(
            window.width / out_w, window.height / out_h)
    return data, transform

Use Resampling.average when decimating a continuous surface such as reflectance or height, and Resampling.nearest for a classified raster, where averaging class codes produces values that mean nothing.

Parameter reference

Setting Where Recommended Why
GDAL_DISABLE_READDIR_ON_OPEN env EMPTY_DIR Stops a directory listing on every open — frequently the dominant cost
CPL_VSIL_CURL_ALLOWED_EXTENSIONS env .tif,.TIF,.tiff Prevents speculative requests for sidecar files that do not exist
GDAL_CACHEMAX env 512 (MB) Holds decoded tiles so overlapping windows are not re-fetched
VSI_CACHE / VSI_CACHE_SIZE env TRUE / 128 MB Caches raw byte ranges, which helps when reading several bands of one scene
GDAL_HTTP_MULTIPLEX env YES Allows HTTP/2 multiplexing where the host supports it
window src.read rounded outward Snaps the read to whole internal tiles
out_shape src.read decimated Serves from overviews rather than full resolution
boundless src.read False True silently pads outside the raster; prefer an explicit intersection

Expected output and verification

The correctness question is whether the window you read covers the geometry you asked for. Compare the georeferenced bounds of the returned array against the geometry’s bounds in the same CRS:

import rasterio
from rasterio.transform import array_bounds


def verify_window(data, profile, wanted_bounds) -> dict:
    """Confirm the returned array actually covers the requested extent."""
    got = array_bounds(profile["height"], profile["width"], profile["transform"])
    covers = (got[0] <= wanted_bounds[0] + 1e-6 and got[1] <= wanted_bounds[1] + 1e-6
              and got[2] >= wanted_bounds[2] - 1e-6 and got[3] >= wanted_bounds[3] - 1e-6)
    return {
        "shape": tuple(data.shape),
        "got_bounds": got,
        "covers_request": bool(covers),
        "all_nodata": bool((data == profile.get("nodata")).all()),
    }


report = verify_window(data, profile, wanted_bounds)
assert report["covers_request"], "window does not cover the request — check the CRS of the bounds"
assert not report["all_nodata"], "everything is NoData — the window is probably in the wrong CRS"

The all_nodata assertion is the one that catches the missing to_crs. A window computed from degrees against a metre-based transform lands in a corner of the scene that is almost always empty, so the read succeeds and returns a perfectly valid array of nothing.

To confirm the read is actually windowed rather than pulling the whole file, measure it. GDAL’s curl layer reports bytes transferred when asked:

import time

import rasterio


def measure_read(cog_url: str, bounds) -> dict:
    """Time a windowed read and report the array size, as a proxy for bytes moved."""
    t0 = time.perf_counter()
    with rasterio.Env(**GDAL_ENV, CPL_CURL_VERBOSE="NO"), rasterio.open(cog_url) as src:
        window = from_bounds(*bounds, transform=src.transform).round_offsets().round_lengths()
        data = src.read(1, window=window)
        full_px = src.width * src.height
    dt = time.perf_counter() - t0
    read_px = int(data.shape[0]) * int(data.shape[1])
    return {"seconds": round(dt, 2), "read_px": read_px, "scene_px": full_px,
            "fraction": round(read_px / full_px, 4)}

A windowed read of a small stand from a Sentinel-2 tile should complete in a second or two over a normal connection. If it takes a minute, the file is almost certainly not a COG — a striped GeoTIFF cannot serve a window, so GDAL falls back to reading everything.

Cost of the same request against three raster layouts Three horizontal bars compare the same stand-sized request. Reading a window from a tiled Cloud Optimized GeoTIFF transfers a small amount and finishes quickly. Reading the identical window from a striped GeoTIFF transfers the whole scene because a strip cannot be partially served. Reading the window from the overview pyramid transfers less again, at reduced resolution. COG, windowed Striped GeoTIFF, same window COG, overview 1:8 14 MB · 1.6 s 690 MB · 74 s 0.4 MB · 0.7 s A slow "windowed" read almost always means the file is not tiled — check with gdalinfo before blaming the network. The request sequence a windowed read actually issues A timeline of requests. A Cloud Optimized GeoTIFF read issues a small header request, then a tile index request, then a batch of range requests for the nine overlapping tiles, which can be issued concurrently. A striped GeoTIFF read issues a single request for the entire file, because a strip spans the full width and cannot be partially served. The total bytes are annotated for both. Cloud Optimized GeoTIFF header index 9 tile ranges, issued concurrently 14 MB total striped GeoTIFF one request for the whole file 690 MB The COG path is more requests and two orders of magnitude fewer bytes.

Common pitfalls

  • Bounds in the wrong CRS. The most common failure and the quietest: no exception, a valid array, and every pixel NoData.
  • Reading band by band without caching. Four bands of one scene means four sets of range requests unless VSI_CACHE is on; with it, the header and shared tiles are fetched once.
  • Using boundless=True to avoid an error. It pads outside the raster with NoData, so a study area that only half-overlaps a scene returns an array that looks complete and is half fabricated. Intersect explicitly instead.
  • Assuming every remote GeoTIFF is a COG. Many public archives still serve striped files. gdalinfo reports Block=...x1 for a striped raster and Block=512x512 for a tiled one.
  • Opening inside a loop without rasterio.Env. The GDAL settings are per-environment; opening outside one silently uses defaults and the directory listing returns.

Frequently Asked Questions

How do I check whether a remote file is really a COG?

gdalinfo /vsicurl/<url> reports the block size and whether overviews exist. Block=512x512 with an Overviews: line is a COG; Block=...x1 means strips and no windowed read is possible. The rio cogeo validate command from rio-cogeo gives a definitive answer including layout order.

Should I read with rasterio.mask instead of a window?

rasterio.mask.mask(src, shapes, crop=True) does the window and the geometry mask in one call and is the better choice when you want both. It reads the geometry’s bounding window, so the byte cost is the same; the difference is only whether you get the pixels outside the polygon masked for you.

Why is my first read slow and later reads fast?

The header and tile index are fetched on open and then cached, so the first request to a scene carries that overhead. Within one rasterio.Env and one process the cache persists; across processes it does not. Reading several windows from the same scene inside a single environment is markedly cheaper than opening it repeatedly.

Up: Cloud-Native Raster Workflows · Ecological GIS Data Foundations in Python