Clipping a Raster to a Bounding Box with Windows

To clip a GeoTIFF to a rectangular bounding box without reading the whole image, convert the box to a pixel window with rasterio.windows.from_bounds, then pass that window to src.read. Only the overlapping blocks are fetched:

import rasterio
from rasterio.windows import from_bounds

left, bottom, right, top = 500000, 4600000, 510000, 4610000  # box in raster CRS

with rasterio.open("scene.tif") as src:
    window = from_bounds(left, bottom, right, top, src.transform)
    clip = src.read(1, window=window)
    transform = src.window_transform(window)

This is the rectangular, memory-lean path inside Automated Image Clipping and Cropping. It is deliberately distinct from How to Clip Rasters to Irregular Polygon Boundaries: a polygon clip masks pixels outside an arbitrary shape, whereas a bounding-box window is an axis-aligned rectangle that GDAL can satisfy by reading a contiguous block. When you only need a rectangular extent, the window approach is far cheaper.

Why This Arises in Remote Sensing Workflows

Satellite scenes are large — a single Sentinel-2 or Landsat tile is hundreds of megabytes — but analysis is usually confined to a study area a fraction of that size. Reading the full array to then slice it wastes memory and, for remote assets, transfers pixels you immediately discard. A windowed read flips that: you translate the study-area bounds into a pixel offset and size, and GDAL fetches only the blocks that intersect it.

The mechanism rests on the affine transform. Every raster carries a geotransform mapping pixel row/column to world coordinates. windows.from_bounds inverts that mapping: given four world coordinates and the transform, it returns the fractional pixel window covering them. Round it to whole pixels, clamp it to the raster extent, and you have a Window object that read and write both understand.

The saving is largest for internally tiled rasters, which includes every well-formed Cloud-Optimized GeoTIFF. Such files store pixels in independent blocks — typically 256×256 or 512×512 — and GDAL can seek directly to the blocks a window overlaps. A 200 MB scene from which you want a 2 km study box might touch only a handful of blocks, so both memory use and, for remote assets, network transfer collapse to a few megabytes. On a striped (non-tiled) raster the same call still works, but GDAL must read whole rows, so the win is smaller; that difference is one more reason to convert archives to tiled COGs before running clip-heavy pipelines.

This topic lives in Automated Image Clipping and Cropping under Satellite Processing Workflows & Index Pipelines. It shares its core mechanic — windowed I/O — with Optimizing Rasterio Window Reads for Memory Efficiency, which goes deeper on block alignment and chunked iteration; here the focus is turning a geographic box into a correct, georeferenced clip.

From World Bounds to a Pixel Window

The diagram shows how four bounding-box coordinates become a pixel window, and why only the intersecting blocks of a tiled raster are read.

Bounding box to pixel window over a tiled raster A geographic bounding box is converted by from_bounds into a pixel window. Overlaid on a tiled raster, only the blocks intersecting the window are read from disk or network; the rest are skipped. Tiled raster (256×256 blocks) bbox from_bounds(...) col_off, row_off width, height (px) src.read(window=window) reads only shaded blocks window_transform → clip origin

The rounded window can extend past the raster if the box overhangs the scene edge. Intersecting it with Window(0, 0, src.width, src.height) clamps it, so the read never asks for out-of-bounds pixels.

Environment & Setup

Package Minimum version Why required
rasterio 1.3.0 windows.from_bounds, window_transform, windowed read/write
numpy 1.24 Array container for the clipped block 
pip install "rasterio>=1.3.0" "numpy>=1.24"

Complete Working Example

This function clips any GeoTIFF (local or a remote COG) to a bounding box, reprojecting the box into the raster CRS if needed, clamping the window to the raster extent, and writing a correctly georeferenced clip.

From map coordinates to a safe pixel window The inverse transform turns the bounding box corners into fractional pixel coordinates. Rounding outward guarantees the window covers the whole request rather than shaving a partial row. Intersecting with the dataset window clamps a box that extends past the scene edge, which is what stops a negative offset reaching read. Three steps, and the one that prevents a crash bbox in CRS units (499980, 9890200, 512340, 9901000) invert fractional pixels col 0.0 … 1236.0 row 47.6 … 1127.4 round out whole pixels col 0 … 1236 row 47 … 1128 window.intersection(dataset_window) clamps a box that runs off the scene — without it, read() gets a negative offset Rounding inward instead of outward is the classic source of a one-pixel gap along tile seams. 
import rasterio
from rasterio.windows import from_bounds, Window
from rasterio.warp import transform_bounds

def clip_to_bbox(
    src_path: str,
    bbox: tuple[float, float, float, float],
    out_path: str,
    bbox_crs: str | None = None,
) -> None:
    """
    Clip a raster to a bounding box using a windowed read.

    bbox : (left, bottom, right, top)
    bbox_crs : CRS of bbox, e.g. "EPSG:4326". If None, bbox is assumed
               to already be in the raster CRS.
    """
    with rasterio.open(src_path) as src:
        left, bottom, right, top = bbox

        # Reproject the box into the raster CRS when they differ
        if bbox_crs is not None and src.crs.to_string() != bbox_crs:
            left, bottom, right, top = transform_bounds(
                bbox_crs, src.crs, left, bottom, right, top
            )

        # World bounds -> fractional pixel window, then snap to whole pixels
        window = from_bounds(left, bottom, right, top, src.transform)
        window = window.round_offsets().round_lengths()

        # Clamp to the raster so we never read outside the array
        full = Window(0, 0, src.width, src.height)
        window = window.intersection(full)
        if window.width <= 0 or window.height <= 0:
            raise ValueError("Bounding box does not overlap the raster extent.")

        data = src.read(window=window)                 # only these blocks are read
        out_transform = src.window_transform(window)   # shift origin to clip corner

        profile = src.profile.copy()
        profile.update(
            height=int(window.height),
            width=int(window.width),
            transform=out_transform,
        )
        with rasterio.open(out_path, "w", **profile) as dst:
            dst.write(data)


if __name__ == "__main__":
    # A WGS84 study box clipped from a UTM-projected scene
    clip_to_bbox(
        src_path="S2_scene_utm.tif",
        bbox=(11.20, 46.05, 11.35, 46.15),  # lon/lat
        out_path="study_area_clip.tif",
        bbox_crs="EPSG:4326",
    )

The two lines that make the clip correct rather than merely cropped are window.intersection(full), which guards against a box that hangs off the scene, and src.window_transform(window), which moves the affine origin so the output carries real-world coordinates for its own top-left corner. Copying src.profile and overriding only height, width, and transform preserves the CRS, dtype, band count, nodata value, and compression of the source, so the clip is a faithful subset rather than a re-encoded approximation.

round_offsets().round_lengths() deserves a note: from_bounds returns a window with fractional offsets because a study-area corner rarely lands exactly on a pixel boundary. Rounding the offsets and lengths to whole pixels snaps the window to the integer grid the array is stored on, which is required before reading — a fractional window has no meaning to read. The trade-off is a sub-pixel shift of at most half a cell, which is immaterial for a rectangular extract and far cheaper than resampling the raster to align it exactly to the box.

Variant Patterns

1. Clipping every band of a multi-band scene

Every case your clipping function will meet A box fully inside returns a normal array. A box overlapping the edge returns a smaller array unless you pad it. A box entirely outside produces an empty window, which should be detected and reported rather than read. A box in a different CRS must be transformed first or it will silently land somewhere else on Earth. Handle all four, or the fourth will find you in production fully inside array as requested overlapping smaller array, or pad with fill disjoint width or height 0 — raise, do not read wrong CRS transform_bounds first Blue is the request, pale blue is the scene. Only the first case works without explicit handling. Assert that the clipped array is non-empty before writing it — an empty COG is a silent pipeline failure.

Omitting the band index in read returns all bands as a 3-D array; the same window applies to every band, and the profile already knows the band count.

import rasterio
from rasterio.windows import from_bounds

with rasterio.open("multiband.tif") as src:
    win = from_bounds(*bbox, src.transform).round_offsets().round_lengths()
    stack = src.read(window=win)             # shape (bands, rows, cols)
    print(stack.shape, src.count)

2. Clipping a remote COG over HTTPS

Because a windowed read only fetches overlapping blocks, the same call against a Cloud-Optimized GeoTIFF transfers just the clip over the network. See How to Read COG Headers Without Downloading Full Files for the header-first validation that pairs with this.

import rasterio
from rasterio.env import Env
from rasterio.windows import from_bounds

url = "https://example-bucket.s3.amazonaws.com/scene_cog.tif"
with Env(GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR"):
    with rasterio.open(url) as src:
        win = from_bounds(*bbox, src.transform).round_offsets().round_lengths()
        subset = src.read(1, window=win)     # only the clip travels over HTTP

3. Snapping the window to the block grid

Aligning the window to the internal tile grid means GDAL reads whole blocks with no partial-block overhead. Round offsets down and lengths up to the block size for the tightest read.

from rasterio.windows import Window

def snap_to_blocks(window: Window, block: int = 256) -> Window:
    col = (int(window.col_off) // block) * block
    row = (int(window.row_off) // block) * block
    w = (int(window.width + window.col_off - col + block - 1) // block) * block
    h = (int(window.height + window.row_off - row + block - 1) // block) * block
    return Window(col, row, w, h)

Common Errors

The clip is offset or lands in the wrong place

The bounding box and raster are in different coordinate systems. Reproject the box with rasterio.warp.transform_bounds into the dataset CRS before calling from_bounds, which assumes the bounds are already in the raster’s coordinate system.

rasterio.errors.WindowError or an empty array

The window fell entirely outside the raster, or its rounded width/height collapsed to zero. Intersect the window with the full-raster window and check window.width > 0 and window.height > 0 before reading.

The written clip has the original raster’s coordinates

You reused src.transform in the output profile. Replace it with src.window_transform(window) so the affine origin points at the clip’s corner, not the parent scene’s corner.

Frequently Asked Questions

Q: Does a windowed read load the whole raster into memory first? No. When you pass window= to src.read, GDAL reads only the blocks overlapping that window. For a tiled GeoTIFF or a remote COG this means memory and network transfer scale with the clip size, not the full image size.

Q: How do I keep the clip georeferenced correctly? Compute the output transform with src.window_transform(window). This shifts the affine origin to the window’s top-left corner so the written clip carries the correct coordinates instead of the original raster origin.

Q: What if my bounding box is in a different CRS than the raster? Reproject the bounding box coordinates into the raster CRS with rasterio.warp.transform_bounds before building the window. from_bounds assumes the bounds are already expressed in the dataset coordinate system.