Cropping a STAC Item to an AOI Without Full Download

To clip a STAC asset to a study area without transferring the whole scene, transform the area’s bounds into the asset’s CRS and read only that window:

import rasterio
from rasterio.warp import transform_bounds
from rasterio.windows import from_bounds

href = item.assets["B04"].href

with rasterio.Env(GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR"):
    with rasterio.open(href) as src:
        left, bottom, right, top = transform_bounds("EPSG:4326", src.crs, *aoi_bounds_wgs84,
                                                    densify_pts=21)
        window = from_bounds(left, bottom, right, top, transform=src.transform)
        arr = src.read(1, window=window, masked=True)
        clip_transform = src.window_transform(window)

This is the network-native version of the clipping workflow in Automated Image Clipping and Cropping.

Why This Arises in Remote Sensing Workflows

Study areas are small and scenes are large. A catchment, a district, a farm block or a park is typically a few kilometres across; a Sentinel-2 tile is 110 km on a side and a Landsat scene 185 km. Downloading the scene to extract a hundredth of it wastes almost all of the transfer, and at time-series scale that waste is the entire cost of the pipeline.

Cloud-optimized assets remove the need. Because the file is tiled and its offsets are indexed in the header, a windowed read touches only the tiles that intersect the window — the mechanism described in Reading a COG over S3 Without Downloading. The work in this page is therefore not about reading efficiently but about getting from a STAC item and a polygon to a correct window: transforming the bounds, rounding them safely, handling partial overlap, and keeping the output georeferenced.

What a windowed clip saves per scene A 12 by 12 kilometre study area occupies about one percent of a Sentinel-2 tile. Downloading the tile transfers around 130 megabytes per band; reading the intersecting tiles transfers about 2 megabytes. Across a 180-date time series and four bands, that is the difference between 94 gigabytes and 1.4 gigabytes. One 12 × 12 km study area inside a 110 km tile AOI ≈ 1.2% of the tile area download the tile 130 MB per band · 94 GB over 180 dates × 4 bands read the window 2 MB per band · 1.4 GB over the same series The saving is not marginal — it changes which pipelines are affordable to run at all. It also removes the local disk requirement, which is what lets the job run in a small container.

Environment & Setup

Package Version Why
pystac-client ≥0.7 Searching and resolving item assets
rasterio ≥1.3.0 Windowed remote reads, transform_bounds
shapely ≥2.0 AOI geometry handling
numpy ≥1.23 Array output 
pip install "pystac-client>=0.7" "rasterio>=1.3.0" "shapely>=2.0"

Complete Working Example

This function takes a STAC item, a list of asset keys and an AOI in WGS84, and returns one clipped array per asset with a shared transform — without transferring anything outside the window.

Outcomes of an AOI against an item footprint An AOI can sit inside the swath, overlap its edge, miss it entirely despite a bounding-box match, or be in a different CRS. Each needs a different response, and only the first works without explicit handling. Four AOI-versus-item outcomes to handle case what the window looks like what the code should do AOI inside the swath full requested size read and return AOI crosses the swath edge clipped, partly nodata read, report the valid fraction bbox matches, footprint misses zero width or height skip before reading AOI in a different CRS wildly wrong offsets transform bounds first, densified The third row is why the item geometry is tested and not only its bounding box. 
import numpy as np
import rasterio
from rasterio.warp import transform_bounds
from rasterio.windows import Window, from_bounds

GDAL_OPTS = {
    "GDAL_DISABLE_READDIR_ON_OPEN": "EMPTY_DIR",
    "GDAL_INGESTED_BYTES_AT_OPEN": "65536",
    "GDAL_HTTP_MULTIRANGE": "YES",
    "GDAL_HTTP_MERGE_CONSECUTIVE_RANGES": "YES",
    "VSI_CACHE": "TRUE",
}


def clip_item(
    item,
    asset_keys: list[str],
    aoi_bounds_wgs84: tuple[float, float, float, float],
    *,
    sign=None,
) -> dict:
    """Clip several assets of one STAC item to an AOI, reading only the window."""
    out: dict[str, object] = {}
    with rasterio.Env(**GDAL_OPTS):
        for key in asset_keys:
            href = item.assets[key].href
            if sign is not None:
                href = sign(href)                    # resolve tokens at read time

            with rasterio.open(href) as src:
                # Densified transformation: a projected box is not the transform of two corners
                left, bottom, right, top = transform_bounds(
                    "EPSG:4326", src.crs, *aoi_bounds_wgs84, densify_pts=21)

                window = from_bounds(left, bottom, right, top, transform=src.transform)
                window = window.round_offsets(op="floor").round_lengths(op="ceil")
                window = window.intersection(Window(0, 0, src.width, src.height))

                if window.width < 1 or window.height < 1:
                    out[key] = None                  # no overlap: report, do not read
                    continue

                out[key] = {
                    "array": src.read(1, window=window, masked=True),
                    "transform": src.window_transform(window),
                    "crs": src.crs,
                    "window": window,
                    "nodata": src.nodata,
                }
    return out


def write_clip(clip: dict, path: str) -> None:
    """Write one clipped array with the window's own transform."""
    arr = clip["array"]
    profile = {
        "driver": "COG", "height": arr.shape[0], "width": arr.shape[1], "count": 1,
        "dtype": arr.dtype, "crs": clip["crs"], "transform": clip["transform"],
        "nodata": clip["nodata"], "compress": "DEFLATE", "blocksize": 512,
    }
    with rasterio.open(path, "w", **profile) as dst:
        dst.write(arr.filled(clip["nodata"] if clip["nodata"] is not None else 0), 1)


if __name__ == "__main__":
    clips = clip_item(item, ["B04", "B08", "SCL"], (36.70, -1.45, 36.82, -1.33))
    for key, clip in clips.items():
        if clip is None:
            print(f"{key}: no overlap")
            continue
        print(key, clip["array"].shape, "valid:", int(clip["array"].count()))
        write_clip(clip, f"clip_{key}.tif")

Two lines do the safety work. round_offsets(floor) with round_lengths(ceil) guarantees the window covers the whole AOI rather than shaving a boundary row — the rounding argument from Clipping a Raster to a Bounding Box with Windowed Reads. And the intersection with the dataset window turns a partially-overlapping AOI into a smaller valid read instead of a negative offset error.

Variant Patterns

1. Search and clip in one pass

Search, filter, clip — in that order A STAC search bounded by the area of interest returns candidate items. Each item's own geometry is tested against the AOI to reject partial-tile misses that the bounding box let through. Only the survivors are opened, and each open reads a single window per asset. Every stage removes work from the next one STAC search bbox + datetime + cloud filter 214 geometry test item.geometry vs AOI costs nothing 168 windowed read one open per asset ~2 MB each clips 168 dates The geometry test matters because a bbox is a rectangle around an irregular footprint — it admits items whose actual swath misses the AOI entirely, and each of those is a wasted open. 
from shapely.geometry import shape, box

aoi = box(36.70, -1.45, 36.82, -1.33)

for item in search.items():
    if not shape(item.geometry).intersects(aoi):
        continue                      # bbox said maybe; the real footprint says no
    clips = clip_item(item, ["B04", "B08"], aoi.bounds)

2. Clipping to a polygon rather than a box

Reading the window is the expensive part; masking to the exact polygon afterwards is free by comparison.

from rasterio.features import geometry_mask

clip = clips["B04"]
inside = geometry_mask([aoi_geom_in_scene_crs],
                       out_shape=clip["array"].shape,
                       transform=clip["transform"],
                       invert=True)
polygon_clip = np.ma.masked_array(clip["array"], mask=clip["array"].mask | ~inside)

Always window first and mask second, for the reason set out in How to Clip Rasters to Irregular Polygon Boundaries.

3. An AOI spanning several tiles

from rasterio.merge import merge

pieces = []
for item in intersecting_items:
    clip = clip_item(item, ["B04"], aoi.bounds)["B04"]
    if clip is not None:
        pieces.append(clip)

# Merge the small clips, not the full tiles
mosaic, transform = merge([p["array"] for p in pieces], transform=[p["transform"] for p in pieces])

Clipping before merging keeps every transfer proportional to the AOI; merging first would assemble two full tiles to discard 99 percent of them, and the blending questions that follow are covered in Seamless Mosaicking and Edge Blending.

Verifying the Clip

Three assertions cover the failure modes that matter, and all are cheap:

import numpy as np
from rasterio.transform import array_bounds

clip = clips["B04"]
arr = clip["array"]

# 1. The clip covers the AOI: its bounds must contain the transformed AOI bounds
bounds = array_bounds(arr.shape[0], arr.shape[1], clip["transform"])
assert bounds[0] <= left and bounds[3] >= top, "window does not cover the AOI"

# 2. There is data: a fully masked clip means the AOI is outside the swath
assert arr.count() > 0, "clip is entirely nodata — check the item footprint"

# 3. The output georeferences correctly: pixel (0,0) maps to the window origin
assert clip["transform"].c == bounds[0] and clip["transform"].f == bounds[3]

The second is the one that catches partial-tile items: a scene whose bounding box overlaps the AOI but whose valid swath does not, which is common at the edge of an orbit.

Common Errors

WindowError: Bounds and transform are inconsistent

The AOI bounds were passed in degrees while the transform is in metres. Transform the bounds into the asset’s CRS first, with densify_pts set.

The clip is one row or column short

The window was rounded inward. Use round_offsets(op="floor") and round_lengths(op="ceil").

The output GeoTIFF is in the wrong place

The scene’s transform was written instead of the window’s. Use src.window_transform(window) for the clipped output.

Frequently Asked Questions

Q: Do I need to download the scene first? No. If the asset is a Cloud-Optimized GeoTIFF, GDAL fetches only the tiles intersecting the window. A 10 km study area inside a Sentinel-2 tile typically transfers a few megabytes instead of a few hundred.

Q: Should I use the item’s bbox to decide whether to read it? Yes, as a first filter — it costs nothing and rejects most non-intersecting items. But the bbox is a rectangle around a possibly irregular footprint, so a positive result still needs the actual read to confirm coverage.

Q: What if the AOI spans two tiles? Clip from each intersecting item and merge the results. Clipping first keeps every read small; merging first would mean assembling the full tiles before discarding most of them.