Extracting Image Chips Around Labelled Points

To read a fixed-size chip centred on a labelled point, convert the point’s coordinates to pixel indices with src.index, offset by half the chip size, and read the window with boundless=True so the shape is guaranteed:

import rasterio
from rasterio.windows import Window

SIZE = 256

with rasterio.open("s2_stack_36NYF.tif") as src:
    row, col = src.index(x, y)                       # map coords -> pixel indices
    win = Window(col - SIZE // 2, row - SIZE // 2, SIZE, SIZE)
    chip = src.read(window=win, boundless=True, fill_value=0)
    chip_transform = src.window_transform(win)       # keeps the chip georeferenced

That is the whole operation. Everything below is about the cases where it goes wrong: points outside the scene, points near the edge, points in the wrong CRS, and datasets of tens of thousands of points where read order starts to matter. The wider workflow this feeds is building training datasets from satellite imagery, inside Raster Machine Learning & Model Inference.

Why Point-Centred Chips Rather Than a Grid

Point labels are what field campaigns, crowd-sourced surveys and expert interpretation actually produce: a location and a class, with no polygon boundary. Chipping on a regular grid over such a dataset wastes almost every window, because labelled locations are sparse and clustered. Chipping around each point instead gives exactly as many samples as there are labels, each one guaranteed to contain its target.

From map coordinates to a centred window A label point in projected coordinates falls inside one pixel of the raster. The index lookup returns that pixel's row and column. Subtracting half the chip size from each gives the window origin, and the resulting window is centred on the labelled pixel with equal context on all four sides. src.index(x, y) then offset by half the chip point at (x, y) lands in pixel (row 2, col 2) src.index returns integer indices, floor-rounded window origin = (row - N//2, col - N//2) equal context on every side of the target

The other advantage is control over context. A grid chip may place the target at a corner where half its surroundings are missing; a centred chip always shows the model the same amount of neighbourhood in every direction, which makes the training distribution consistent.

Environment & Setup

Package Version pin Used for
rasterio >=1.3.0 index, Window, boundless reads, window_transform
geopandas >=0.14 Reading the point layer and reprojecting it
numpy >=1.23 Chip arrays and the validity mask 
pip install "rasterio>=1.3.0" "geopandas>=0.14" "numpy>=1.23"

Complete Working Example

This function takes a point layer and a scene, and returns chips with their transforms, skipping points that fall outside the scene entirely and flagging those that needed padding.

from dataclasses import dataclass

import geopandas as gpd
import numpy as np
import rasterio
from rasterio.windows import Window


@dataclass
class Chip:
    array: np.ndarray
    transform: object
    class_id: int
    padded: bool


def chips_from_points(raster_path: str, points_path: str, *,
                      size: int = 256, class_field: str = "class_id",
                      scale: float = 10_000.0) -> list[Chip]:
    points = gpd.read_file(points_path)
    out: list[Chip] = []

    with rasterio.open(raster_path) as src:
        if points.crs != src.crs:
            points = points.to_crs(src.crs)          # always before src.index

        fill = src.nodata if src.nodata is not None else 0
        half = size // 2

        # Sorting by row keeps successive reads close together in the file,
        # which matters enormously for a remote COG.
        points = points.assign(
            _row=[src.index(g.x, g.y)[0] for g in points.geometry],
            _col=[src.index(g.x, g.y)[1] for g in points.geometry],
        ).sort_values("_row")

        for row, col, cls in zip(points._row, points._col, points[class_field]):
            # Skip points that are not in this scene at all
            if not (0 <= row < src.height and 0 <= col < src.width):
                continue

            win = Window(col - half, row - half, size, size)
            padded = (col - half < 0 or row - half < 0
                      or col + half > src.width or row + half > src.height)

            arr = src.read(window=win, boundless=True, fill_value=fill)
            chip = arr.astype("float32") / scale
            chip[:, arr[0] == fill] = np.nan         # padding is not an observation

            out.append(Chip(chip, src.window_transform(win), int(cls), padded))

    return out


if __name__ == "__main__":
    chips = chips_from_points("s2_stack_36NYF.tif", "field_points.gpkg")
    padded = sum(c.padded for c in chips)
    print(f"{len(chips)} chips, {padded} needed padding at the scene edge")
    print("chip shape:", chips[0].array.shape)

Three lines carry most of the correctness. The to_crs call happens before any index lookup, because src.index interprets its arguments in the dataset’s CRS and will silently return nonsense indices for degrees fed to a UTM scene. The padded flag records which chips touched the edge, so they can be excluded or down-weighted later. And converting padding to NaN rather than leaving it as zero stops the model learning that “reflectance exactly zero” is a meaningful feature.

Variant Patterns

1. Jittered chips for augmentation

Perfectly centred targets teach a model that the answer is always in the middle. Offsetting the window by a random amount within a fraction of the chip breaks that, and it multiplies a small point dataset into a usable one.

Jittering the window around a point On the left one window is centred exactly on the label point, so the target always appears in the middle of the chip. On the right four windows are offset by random amounts up to a quarter of the chip size, so the same point appears at different positions and the model cannot rely on position alone. centred only jitter up to N/4 one sample per point four samples, target in four positions 
import numpy as np
from rasterio.windows import Window


def jittered_window(row: int, col: int, size: int, rng, max_shift: int | None = None) -> Window:
    """A window containing (row, col) but not necessarily centred on it."""
    max_shift = max_shift if max_shift is not None else size // 4
    dr = int(rng.integers(-max_shift, max_shift + 1))
    dc = int(rng.integers(-max_shift, max_shift + 1))
    return Window(col - size // 2 + dc, row - size // 2 + dr, size, size)

Keep the jitter well inside half the chip size, or the labelled pixel can fall outside the window entirely — a silent corruption that produces chips labelled for something they do not contain.

2. Reading the same points from a multi-date stack

For temporal models the chip is a cube: the same window across many dates. Reading date by date and stacking keeps memory bounded and lets missing dates be skipped rather than failing the whole point.

import numpy as np
import rasterio


def temporal_chip(paths: list[str], x: float, y: float, size: int = 64) -> np.ndarray:
    """Stack the same spatial window across dates into (time, bands, rows, cols)."""
    frames = []
    for path in paths:
        with rasterio.open(path) as src:
            row, col = src.index(x, y)
            win = rasterio.windows.Window(col - size // 2, row - size // 2, size, size)
            frames.append(src.read(window=win, boundless=True,
                                   fill_value=src.nodata or 0).astype("float32"))
    return np.stack(frames)

Every scene in the list must share a grid for this to be meaningful; if they do not, align them first with the approach in aligning two rasters with reproject_match.

3. Chipping directly from a remote COG

Nothing changes except the path. A windowed read against an HTTPS or S3 URL fetches only the overlapping byte ranges, which is what makes chipping from a public archive practical — the mechanics are covered in reading a COG over S3 without downloading.

Bytes moved per chip from a remote scene Downloading a one gigabyte scene to extract a single chip moves the whole file. A windowed read against the same Cloud-Optimized GeoTIFF issues one header request and a small number of range requests covering only the internal tiles the chip overlaps, moving a few hundred kilobytes. One 256 px chip from a 1 GB scene download then read 1,024,000 KB transferred windowed read ~340 KB transferred — header plus four internal tiles request 1: first 16 KB — header, tile offsets and sizes requests 2-5: the byte ranges of the four 512 px tiles the chip overlaps The ratio only holds if the file is internally tiled — a striped GeoTIFF fetches full-width rows instead. 
import os

import rasterio

os.environ["GDAL_DISABLE_READDIR_ON_OPEN"] = "EMPTY_DIR"   # no directory listing
os.environ["CPL_VSIL_CURL_ALLOWED_EXTENSIONS"] = ".tif"

with rasterio.open("https://example-bucket.s3.amazonaws.com/s2/36NYF_B08.tif") as src:
    row, col = src.index(x, y)
    chip = src.read(window=rasterio.windows.Window(col - 128, row - 128, 256, 256),
                    boundless=True, fill_value=0)

Common Errors

IndexError or wildly out-of-range indices from src.index

The points are in a different CRS from the raster. src.index does no transformation — it assumes its arguments are already in the dataset’s coordinate system. Reproject the layer first and assert points.crs == src.crs.

WindowError: Bounds and transform are inconsistent

A window with a negative offset was passed without boundless=True. Either pass boundless=True with a fill_value, or intersect the window with the dataset extent first.

Chips near the edge are the wrong shape

boundless was omitted, so rasterio clipped the window to the dataset. Any batched model will fail on the first such chip. Always read boundlessly and track the padding flag.

Frequently Asked Questions

Q: Why does my chip look shifted by one pixel? An even-sized chip has no centre pixel, so the point falls on a pixel corner. For a chip of size N the conventional offset is row minus N//2, which puts the point at index N//2 — the first pixel past the middle. Use an odd size when the exact centre pixel matters.

Q: How do I chip from a remote COG without downloading it? Open the HTTPS or S3 URL directly. A windowed read fetches only the byte ranges covering that window, so a 256 pixel chip from a 1 GB scene transfers a few hundred kilobytes. Sort the points by row before reading so requests hit nearby byte ranges.

Q: Should the point be at the chip centre or anywhere inside it? Centre it for a first dataset, then add jitter. A model trained only on perfectly centred targets can learn that the answer is always in the middle, which does not hold at inference time when windows tile the scene.