Raster Data Handling with Rasterio: Production Workflows & Pipelines

Rasterio is the pixel-processing tier of the Python geospatial stack, and this guide sits inside Mastering Core Geospatial Python Libraries as the counterpart to the vector-focused guides. Where GeoPandas DataFrames and Shapely geometry operations carry tabular geometry, Rasterio carries the grid: elevation models, multispectral imagery, land-cover classifications, and any other regularly-sampled surface. It wraps GDAL in a context-managed, NumPy-native API so you can read a sub-region of a terabyte-scale scene without loading the whole array, run band arithmetic as vectorized operations, mask against vector boundaries, and write standards-compliant output. This guide covers the full production arc — from environment setup and metadata inspection through band math, CRS alignment, and Cloud Optimized GeoTIFF export — with the memory and projection discipline that separates a reproducible pipeline from an ad-hoc script.

Windowed reads across raster bands A multi-band raster shown as stacked band layers and a tiled grid, with one small window highlighted to show that Rasterio reads a sub-region without loading the whole array. Read a window, not the whole array bands R · G · B · NIR window Window
Rasterio reads any band and any window independently — the key to processing rasters larger than memory.

Architecture & Data Structures

Every Rasterio session revolves around a DatasetReader (or DatasetWriter) obtained from rasterio.open(). The object is a thin, lazy handle over a GDAL dataset: opening it reads only the header, not the pixels. Three attributes define the spatial contract of any raster and you will touch them on almost every page:

import rasterio

# A Sentinel-2 surface-reflectance scene reprojected to UTM
scene_path = "data/sentinel2_utm33n.tif"

with rasterio.open(scene_path) as src:
    print(f"Driver:    {src.driver}")          # GTiff
    print(f"Size:      {src.width} x {src.height}")
    print(f"Bands:     {src.count}")            # 1-based band indices
    print(f"Dtype:     {src.dtypes[0]}")        # e.g. uint16
    print(f"CRS:       {src.crs}")              # EPSG:32633
    print(f"Transform: {src.transform}")        # Affine(10.0, 0.0, ...)
    print(f"NoData:    {src.nodata}")

Two properties of this model matter for everything downstream. First, band indices are 1-based — src.read(1) is the first band, and passing 0 raises IndexError. Second, the context manager guarantees the underlying file handle and GDAL resources are released on exit, which prevents OS-level file locks during batch processing and keeps GDAL's block cache from ballooning. For a deeper comparison of this API against the raw GDAL bindings, see Rasterio vs GDAL Python Bindings.

Pixel coordinates, world coordinates, and the half-pixel question

The Affine transform is the only thing standing between an array index and a place on Earth, and Rasterio exposes the conversion in both directions. src.index(x, y) turns projected coordinates into (row, col); src.xy(row, col) goes back. The detail that costs people an afternoon is that these are not exact inverses by default, because xy() returns the centre of the pixel while index() floors into whichever pixel contains the point:

import rasterio

with rasterio.open("data/dem_utm33n.tif") as src:
    row, col = src.index(692_450.0, 4_981_200.0)   # world -> array index (floored)
    cx, cy = src.xy(row, col)                      # array index -> PIXEL CENTRE
    ulx, uly = src.xy(row, col, offset="ul")       # array index -> upper-left corner

    print(row, col)          # e.g. 480 249
    print(cx - 692_450.0)    # a fraction of a pixel, not 0.0
    print(src.res)           # (10.0, 10.0) — so the offset is under 10 m

Round-tripping through xy() and index() is therefore stable, but comparing xy() output against your original coordinates is not — and a pipeline that snaps sample points by that comparison drifts by half a pixel in each direction. Two related helpers fall out of the same transform: rasterio.windows.from_bounds(*bounds, transform=src.transform) builds a Window from a map extent instead of array offsets, and src.window_transform(win) produces the transform for that sub-array so a clipped output stays georeferenced. Never construct the derived transform by adding offsets to the parent's c and f coefficients by hand.

Metadata that changes the numbers: scales, offsets and tags

A DatasetReader also carries per-band metadata that silently rescales your analysis if you ignore it. src.scales and src.offsets hold the linear packing factors that many surface-reflectance and climate products use to store float measurements as int16 — a Sentinel-2 L2A band with scales=(0.0001,) means the raw 2431 you read is a reflectance of 0.2431, not 2431. Rasterio does not apply these automatically on read(); you must multiply. src.descriptions names each band, src.units records the physical unit where the writer bothered to set one, and src.tags() returns the dataset-level key/value metadata (acquisition date, processing baseline, sensor). src.tags(1) returns the same for band 1. Copying src.profile forward preserves the georeferencing but not these tags, so re-apply them with dst.update_tags(**src.tags()) when a derived product must stay traceable.

Environment Configuration & Dependency Resolution

Rasterio is a compiled extension linked against GDAL, which in turn links PROJ and GEOS. ABI compatibility across these C libraries is the single most common cause of a broken install, so pin the whole stack from one channel rather than mixing pip wheels with system packages.

# conda-forge ships GDAL/PROJ/GEOS binaries built against one another
conda create -n raster-pipeline python=3.12 \
    "rasterio>=1.3" "gdal>=3.6" numpy pyproj geopandas shapely rio-cogeo \
    -c conda-forge
conda activate raster-pipeline

Production considerations before you process anything at scale:

The pip wheels bundle their own GDAL and work well for pure-raster jobs, but the moment you also need GeoPandas or PostGIS drivers in the same environment, conda-forge's single-channel resolution saves hours of ABI debugging.

Which creation options you actually have is a property of the build

Rasterio declares a minimum GDAL, not an exact one, so two environments that both satisfy rasterio>=1.3 can offer different capabilities. Three differences bite in practice. The dedicated COG driver — which writes a compliant Cloud Optimized GeoTIFF in one step instead of the manual tile-and-overview recipe below — arrived in GDAL 3.1, so on anything older driver="COG" fails with an unhelpful driver-not-found error. The ZSTD and LERC compressors are compile-time options: a GDAL built without them accepts compress="zstd" in your profile and then errors at write time, which in a nightly batch means you discover it after the expensive read. And num_threads="ALL_CPUS" on a write is only honoured for compressors that support parallel encoding.

Rather than guessing, ask the driver what it supports before you commit a profile:

import rasterio

print(rasterio.__gdal_version__)          # e.g. 3.9.2 — the GDAL Rasterio is linked against

with rasterio.Env() as env:
    opts = env.drivers()                   # driver short name -> long name
    print("GTiff" in opts, "COG" in opts)  # True True on GDAL >= 3.1

from osgeo import gdal                      # same environment, same libgdal
spec = gdal.GetDriverByName("GTiff").GetMetadataItem("DMD_CREATIONOPTIONLIST")
print("ZSTD" in spec, "LERC" in spec)      # build-dependent, not version-dependent

Two GDAL environment variables are worth setting deliberately in any container image. GDAL_CACHEMAX (megabytes as a plain integer, or a percentage with %) sizes the block cache — the default is small enough that a windowed loop over a striped source re-reads the same compressed strips repeatedly. GDAL_NUM_THREADS=ALL_CPUS parallelises deflate compression on write and warping. Set both per-process in the image, and scope any per-job overrides inside a rasterio.Env() block so they do not leak into the next task in the same worker.

Vectorized Operations & Core Workflow

The canonical Rasterio workflow is: open, validate, read into NumPy, compute vectorized, write back. Because a band read returns a plain ndarray, all band math is ordinary NumPy — no per-pixel Python loops. The function below inspects a raster's spatial contract, fails loudly on a missing CRS, and returns the profile so callers can write derived products with identical georeferencing.

import rasterio
from rasterio.errors import CRSError


def inspect_raster(src_path: str) -> dict:
    """Validate georeferencing and report the array's numeric envelope."""
    with rasterio.open(src_path) as src:
        # A missing CRS silently corrupts every downstream spatial operation
        if not src.crs:
            raise CRSError(
                f"{src_path} has no CRS. Assign one on open with "
                "rasterio.open(path, crs=CRS.from_epsg(32633)) or supply a "
                ".prj / .aux.xml sidecar before processing."
            )

        print(f"{src.height} x {src.width}, {src.count} band(s)")
        print(f"CRS EPSG: {src.crs.to_epsg()}  |  pixel size: {src.res}")

        # Read one band; nodata is preserved as-is, not masked yet
        elevation = src.read(1)
        valid = elevation[elevation != (src.nodata if src.nodata is not None else elevation.min())]
        print(f"dtype {elevation.dtype}, valid range {valid.min()}–{valid.max()}")

        return src.profile.copy()


# profile = inspect_raster("data/dem_utm33n.tif")

Two habits keep this loop reliable. Always carry src.profile forward when writing (profile.update(dtype="float32", count=1)) so the output inherits the driver, nodata, transform, and CRS instead of silently dropping them. And when a band participates in arithmetic, cast to float32 first — integer DEMs and uint16 reflectance overflow or truncate the moment you subtract or divide. The band-selection and memory-safe reading patterns are covered in depth in Reading Multi-Band TIFFs with Rasterio.

Where a single-process windowed loop stops scaling

It helps to know the arithmetic before you hit the wall. A full Sentinel-2 10 m tile is 10,980 × 10,980 pixels; at uint16 that is 241 MB per band, so a twelve-band read materialises 2.9 GB before you have done any maths — and casting to float32 for an index doubles it again. Windowed reads keep the resident set at roughly chunk² × bands × itemsize, which is why a 1024-pixel chunk over four bands stays under 35 MB regardless of how large the scene is.

What the loop does not solve is wall-clock time on a catalogue of scenes. A single Python process spends most of a windowed NDVI run inside GDAL's decompressor, one block at a time, and that is a single core. Scaling out has three rules that are easy to get wrong:

The practical crossover: below roughly a dozen scenes, a single windowed loop is simpler and fast enough. Above that, parallelise across files with a process pool rather than across windows within a file, because per-scene work is embarrassingly parallel and needs no shared state. When the same scenes must also be aligned along a time axis, the labelled-array route in xarray & rioxarray raster cubes replaces the loop with a chunked, Dask-scheduled graph.

Geometry / Data Processing Details

Two processing patterns dominate real raster work: clipping a raster to vector boundaries, and computing spectral indices over tiles too large for RAM. Both stay vectorized and both hinge on getting the geometry–pixel relationship right.

Masking a raster to vector boundaries

rasterio.mask.mask burns polygons into the pixel grid, keeping pixels inside the geometry and setting the rest to nodata. It expects GeoJSON-like mappings in the raster's CRS, so reproject the vector first — this is the most common masking bug. Pairing it with Shapely geometry operations lets you validate and repair boundaries before they touch the raster.

import geopandas as gpd
import rasterio
from rasterio.mask import mask
from shapely.geometry import mapping


def clip_raster_to_catchment(raster_path, catchment_path, out_path):
    catchment = gpd.read_file(catchment_path)
    catchment = catchment[catchment.geometry.is_valid]
    if catchment.empty:
        raise ValueError("No valid catchment geometries to mask with.")

    with rasterio.open(raster_path) as src:
        # Reproject vector to the raster CRS — mask() does NOT do this for you
        if catchment.crs and not catchment.crs.equals(src.crs):
            catchment = catchment.to_crs(src.crs)

        shapes = [mapping(geom) for geom in catchment.geometry]
        clipped, clipped_transform = mask(
            src, shapes, crop=True, all_touched=True, nodata=src.nodata
        )

        out_profile = src.profile.copy()
        out_profile.update(
            height=clipped.shape[1],
            width=clipped.shape[2],
            transform=clipped_transform,
            compress="deflate",
            tiled=True,
        )
        with rasterio.open(out_path, "w", **out_profile) as dst:
            dst.write(clipped)

Use all_touched=True to keep any pixel the polygon boundary crosses (otherwise thin features lose edge pixels). When a catchment falls entirely outside the raster extent, mask raises ValueError: Input shapes do not overlap raster — catch it and skip. For multipart features, catchment.explode(ignore_index=True) first so each part masks independently. Loading the vector boundaries themselves is a GeoPandas job.

The same catchment polygon burned into a pixel grid under both rasterisation rules Two identical six-by-five pixel grids each carry the same diamond-shaped catchment boundary. On the left, all_touched is False, so a pixel is kept only when its centre — shown as a small dot — falls inside the polygon; eight pixels survive and the top and bottom rows are lost entirely even though the boundary passes through them. On the right, all_touched is True, so every pixel the boundary crosses is kept, giving twenty pixels: the same eight core pixels plus twelve edge pixels in a lighter tint. The difference is entirely at the boundary, which is why narrow features vanish under the default rule. One catchment, two rasterisation rules — the difference is all edge all_touched=False all_touched=True centre inside boundary only catchment pixel centre 8 pixels kept — centre-in-polygon test 20 pixels kept — every pixel the boundary crosses A feature narrower than one pixel has no centres to catch, so the default rule returns an empty mask.
The default centre-in-polygon rule discards every pixel the boundary merely clips, which is why thin catchments and narrow corridors come back empty until all_touched=True widens the burn to the full boundary.

Windowed band math for large scenes

Spectral indices like NDVI are a subtraction and a division — trivial arithmetic, but a full Sentinel-2 tile stack will not fit in memory. Iterate over windows aligned to the file's native tiling and write each result block as you go.

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


def ndvi_windowed(scene_path, out_path, chunk=1024):
    with rasterio.open(scene_path) as src:
        # Sentinel-2 band order here: 4 = Red, 8 = NIR (1-based indices)
        red_idx, nir_idx = 4, 8

        profile = src.profile.copy()
        profile.update(count=1, dtype="float32", nodata=-9999.0)

        with rasterio.open(out_path, "w", **profile) as dst:
            for row in range(0, src.height, chunk):
                for col in range(0, src.width, chunk):
                    win = Window(
                        col, row,
                        min(chunk, src.width - col),
                        min(chunk, src.height - row),
                    )
                    red = src.read(red_idx, window=win).astype("float32")
                    nir = src.read(nir_idx, window=win).astype("float32")

                    denom = nir + red
                    ndvi = np.where(denom == 0, -9999.0, (nir - red) / denom)
                    dst.write(ndvi.astype("float32"), 1, window=win)

Guarding the 0/0 case with np.where avoids RuntimeWarning: invalid value encountered in divide and keeps nodata explicit. Align chunk to src.block_shapes[0] so each read maps onto native TIFF tiles and you never re-read the same block. Band numbering is sensor-specific — confirm it against the product spec before trusting an index.

Crossing the raster–vector boundary in both directions

rasterio.features is the bridge between the grid and the feature table, and it goes both ways. rasterize() burns attribute values from geometries into a new array on a transform you supply — the standard way to build a zone raster whose pixel values are parcel identifiers, so later analysis is a NumPy bincount rather than a per-polygon mask. shapes() does the inverse, walking a classified array and yielding GeoJSON polygons for each connected run of equal values.

import geopandas as gpd
import numpy as np
import rasterio
from rasterio.features import rasterize, shapes
from shapely.geometry import shape

with rasterio.open("data/landcover_utm33n.tif") as src:
    grid_transform, grid_shape = src.transform, (src.height, src.width)
    landcover = src.read(1)
    raster_crs = src.crs

parcels = gpd.read_file("data/urban_parcels.gpkg").to_crs(raster_crs)

# 1. Vector -> raster: burn a small integer key, never a float id
parcels["zone_key"] = np.arange(1, len(parcels) + 1, dtype="int32")
zone_raster = rasterize(
    ((geom, key) for geom, key in zip(parcels.geometry, parcels.zone_key)),
    out_shape=grid_shape,
    transform=grid_transform,
    fill=0,                 # 0 means "no parcel here"
    dtype="int32",
    all_touched=False,      # centre-in-polygon, matching mask() semantics
)

# 2. Raster -> vector: polygonise one land-cover class only
water = (landcover == 5).astype("uint8")
water_polys = [
    shape(geom)
    for geom, value in shapes(water, mask=water.astype(bool), transform=grid_transform)
    if value == 1
]
water_gdf = gpd.GeoDataFrame({"class_id": 5}, geometry=water_polys, crs=raster_crs)
print(len(water_gdf), "water polygons")

Three things go wrong here reliably. The geometries must already be in the raster's CRS — rasterize takes a transform, not a projection, so a mismatched vector burns into empty space and returns an all-fill array with no error. The burn value must fit the dtype you declare: passing a float parcel id into an int32 output truncates it, and passing a value above 255 into uint8 wraps around. And shapes() on an unclassified continuous raster produces one polygon per distinct float value, which on a DEM means millions of single-pixel squares — always classify or threshold first, and pass a boolean mask so background is skipped entirely.

The output of the first direction is exactly what per-zone summaries consume; the statistics themselves belong to zonal statistics and raster sampling. The output of the second is an ordinary vector layer, so run it through the topology validation and repair checks before publishing — polygonised rasters have stair-stepped boundaries and frequently share edges only approximately.

CRS Alignment & Projection Pipeline

Raster reprojection resamples pixels onto a new grid, so it needs a source CRS, a target CRS, and a freshly computed transform. Getting the CRS objects right is where most errors hide, and the rules echo those in Coordinate Systems with PyProj:

import rasterio
from rasterio.crs import CRS
from rasterio.warp import calculate_default_transform, reproject, Resampling


def reproject_scene(src_path, out_path, dst_epsg=25832):
    """Reproject to an appropriate projected CRS (ETRS89 / UTM 32N here)."""
    dst_crs = CRS.from_epsg(dst_epsg)
    with rasterio.open(src_path) as src:
        transform, width, height = calculate_default_transform(
            src.crs, dst_crs, src.width, src.height, *src.bounds
        )
        profile = src.profile.copy()
        profile.update(crs=dst_crs, transform=transform, width=width, height=height)

        with rasterio.open(out_path, "w", **profile) as dst:
            for band in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, band),
                    destination=rasterio.band(dst, band),
                    src_transform=src.transform, src_crs=src.crs,
                    dst_transform=transform, dst_crs=dst_crs,
                    # bilinear for continuous data; nearest for classes
                    resampling=Resampling.bilinear,
                )

Choose the resampling kernel by data semantics: Resampling.bilinear or cubic for continuous surfaces (elevation, reflectance), and Resampling.nearest for categorical rasters (land cover, classified masks) so class codes are never averaged into nonsense values. When the reprojection also has to survive limited RAM, the windowed strategy in Reprojecting Large Datasets Without Memory Errors applies directly.

Three failure modes cluster around this call. First, nodata has to be declared on both sides: reproject() accepts src_nodata and dst_nodata, and if the source band has nodata=None the warper treats fill as real data and smears it into the valid area along every edge — a one-to-two pixel halo that grows with the resampling kernel width. Copying src.profile carries the nodata forward only when the source actually set it; check src.nodata is not None and pass the value explicitly rather than assuming.

Second, calculate_default_transform picks a resolution for you, derived from the source pixel size warped into the target CRS. That is convenient and almost never what you want when several scenes must align: two adjacent tiles reprojected independently end up on grids offset by a fraction of a pixel, and every later stack, difference or zonal summary is silently interpolated. Pass resolution=(10.0, 10.0) (or explicit dst_width/dst_height) so the whole batch lands on one grid, and reproject a categorical mask with Resampling.nearest even when the paired continuous band uses bilinear.

Third, a UTM zone is a bad target CRS for anything spanning a zone boundary. A catchment straddling the 12°E meridian reprojected into zone 32N pushes the eastern half far outside the zone's valid extent, where scale error grows quickly; the same applies to any scene crossing the antimeridian, where calculate_default_transform can return a bounding box spanning the whole globe. Pick a regional projected CRS or an equal-area projection for those cases — the selection logic is in choosing a UTM zone automatically in Python.

When you only need aligned pixels for one downstream read and not a file on disk, WarpedVRT presents a reprojected view of a dataset without materialising it:

import rasterio
from rasterio.vrt import WarpedVRT
from rasterio.enums import Resampling

with rasterio.open("data/sentinel2_utm33n.tif") as src:
    with WarpedVRT(
        src,
        crs="EPSG:25832",
        resampling=Resampling.bilinear,
        src_nodata=src.nodata,
        nodata=src.nodata,
    ) as vrt:
        # Reads are warped on the fly, window by window — nothing is written
        print(vrt.crs, vrt.width, vrt.height)
        patch = vrt.read(1, window=rasterio.windows.Window(0, 0, 1024, 1024))

This is the right tool when a pipeline reads a small fraction of a large scene in a foreign CRS: you pay the warp only for the windows you touch. It is the wrong tool when every pixel will be read more than once, because a VRT re-warps on each access instead of caching a materialised grid.

Raster reprojection and COG delivery pipeline A left-to-right pipeline: a source GeoTIFF in native UTM feeds a warp step, then a tiled compressed write and a COG validation, all inside a native-UTM analysis zone where distances and areas stay valid. A dashed divider marks a single reprojection to EPSG:3857 that produces the web tile server output, labelled display only. Reproject once, at the end — analysis in UTM, display in Web Mercator Source GeoTIFF native UTM EPSG:32633 Warp calc_transform + resampling Tiled write deflate · 512px overviews COG validate rio cogeo validate Web tile server EPSG:3857 HTTP range reads reproject → 3857 Native UTM (EPSG:32633) — metric-safe NDVI · slope · zonal statistics computed here Web Mercator display only
Every metric operation runs in the scene's native UTM zone; the pipeline reprojects to Web Mercator exactly once, to produce the display tiles — never for measurement.

Production Export & Integration

The de-facto delivery format for analysis-ready rasters is the Cloud Optimized GeoTIFF (COG): an internally tiled, overview-bearing GeoTIFF that HTTP range requests can read partially. Writing one is a matter of the right creation options — tiling, compression, and internal overviews.

import rasterio
from rasterio.enums import Resampling
from rasterio.shutil import copy as rio_copy

def write_cog(src_path, cog_path):
    with rasterio.open(src_path) as src:
        profile = src.profile.copy()
        profile.update(
            driver="GTiff", tiled=True, blockxsize=512, blockysize=512,
            compress="deflate", predictor=2, interleave="pixel",
        )
        with rasterio.open(cog_path, "w", **profile) as dst:
            dst.write(src.read())
            # Internal overviews power fast zoomed-out reads
            dst.build_overviews([2, 4, 8, 16], Resampling.average)
            dst.update_tags(ns="rio_overview", resampling="average")

A short export checklist before shipping:

Once published to object storage, the same COG feeds partial reads over HTTP — see Windowed Reads from Cloud Optimized GeoTIFF in the ingestion and processing workflows section. For catalogs of rasters queried alongside vector footprints, load the extents into PostGIS and let spatial analysis and advanced query techniques drive selection before you touch pixels.

Windows / Platform Edge Cases & Debugging

Most Rasterio failures are environment failures, and they concentrate on Windows and in cloud runners where the GDAL/PROJ data paths are not where the wheel expects them.

Frequently Asked Questions

When should I stop using Rasterio and move to xarray? The moment the third dimension becomes meaningful. Rasterio's model is one file, one grid, bands addressed by integer; as soon as you are keeping a dictionary of dates to arrays, or writing if band == 8 in more than one place, you are hand-rolling what a labelled cube gives you for free. The crossover is discussed concretely in xarray vs Rasterio for time-series rasters. For single-scene reads, masking, and writing, Rasterio stays the lighter and more predictable option.

Do I need to reproject a raster before clipping it to a vector layer? No — reproject the vector instead. rasterio.mask.mask needs the geometries in the raster's CRS, and reprojecting a few hundred polygons is instantaneous, while warping a multi-gigabyte scene resamples every pixel and costs you accuracy you cannot get back. Only reproject the raster when the output itself has to be in a different CRS.

Why is my windowed loop slower than reading the whole file? Almost always because the windows do not match the file's internal layout. A striped GeoTIFF has block_shapes like (1, 20000) — one row per block — so a square 1024 × 1024 window forces GDAL to decompress 1024 full-width strips and throw most of each away. Check src.block_shapes[0] first; if the source is striped and you need square windows repeatedly, rewrite it once as a tiled file and read from that.

Is src.read() lazy? rasterio.open() is lazy — it reads the header only. read() is not: it allocates and fills a NumPy array immediately. There is no lazy array in Rasterio's model, which is precisely the gap WarpedVRT, windowed loops, and the Dask-backed cubes in xarray & rioxarray raster cubes exist to fill.

Should I store analysis outputs as COG or as NetCDF/Zarr? COG when the consumer is a map client or another GDAL-based tool reading spatial windows over HTTP, because range requests plus internal overviews are exactly what browsers and tile servers need. NetCDF or Zarr when the consumer slices along time or variable — a COG has no time axis, so a 200-date series becomes 200 files and every temporal query becomes a file listing. The trade-offs across the whole family are in cloud-native geospatial formats.

How do I know a derived raster is still correctly georeferenced? Assert on the transform and CRS, not on the array. After any write, reopen the output and check dst.crs == src.crs, dst.transform == expected_transform, and that dst.bounds overlaps the input bounds by the amount you intended. A raster that is shifted by one pixel looks completely normal in a viewer and only surfaces as a systematic bias in zonal statistics months later.

Related