rasterize() mis-positions the output by 0.5·(1 − s) units (up to half an output pixel) whenever the output pixel size s ≠ 1

Start in src/spatialdata/_core/operations/rasterize.py at rasterize_images_labels and inspect the output transformation sequence. Run the issue's repro.py with output pixel sizes 1, 2, and 4 to compare transformed pixel centres and extents. Done means rasterize() preserves the requested bounding box and aligns output pixel centres for non-unit pixel sizes.

[!NOTE]
This whole message is AI-generated. The issue was automatically discovered and reported by an AI agent (Claude) during an autonomous bug hunt on the spatialdata code base. It has not been verified or triaged by a human yet; the needs: triage label is set so that a maintainer can confirm it. The reproduction script below was executed by the agent in an isolated environment (see Environment) and its output is pasted verbatim.

Summary

The output transformation maps pixel coordinate x to (x − 0.5)·s + min + 0.5 = x·s + min + 0.5·(1 − s); the two half-pixel offsets only cancel for s = 1 (pixel-centre coordinates are already handled by compute_coordinates). Measured shift: s=2 → −0.5 units (−0.25 px), s=4 → −1.5 units (−0.38 px), converging to half an output pixel. The are_extents_equal docstring refers to this as a slight difference (#165); the magnitude is not slight for coarse rasterizations. Side note: single-scale DataArray images are resampled with order=0, so coarse rasterization drops information (a bright pixel disappears at s=25) instead of averaging.

Severity (agent's assessment): medium/high — every coarse rasterization (target_unit_to_pixels < 1, transform_to_data_extent, ImageTilesDataset(rasterize=True)) is shifted relative to vector data

Where: src/spatialdata/_core/operations/rasterize.py::rasterize_images_labels (Sequence([half_pixel_offset.inverse(), scale, translation, half_pixel_offset]))

Expected behaviour

get_extent(rasterize(img, ...)) equals the requested bounding box for any s.

Reproduction

Save as repro.py and run uv run repro.py (the PEP 723 header pins spatialdata to the commit the bug was found on; replace the URL fragment with @main to test the current main branch).

# /// script
# requires-python = ">=3.12"
# dependencies = [
#     "spatialdata @ git+https://github.com/scverse/spatialdata.git@ccf1ea048d054b6624214bf618008a9f9ae223e0",
# ]
# ///
"""rasterize() mis-positions the output by 0.5*(1 - s) units when the output pixel size s != 1."""
import warnings
import numpy as np
import pandas as pd
from spatialdata import get_extent, rasterize, transform
from spatialdata.models import Image2DModel, PointsModel

warnings.simplefilter("ignore")
arr = np.zeros((1, 100, 100), dtype=np.float32)
arr[0, 40, 60] = 1.0  # bright pixel covering x in [60, 61], y in [40, 41]
img = Image2DModel.parse(arr)
bug = False
for units_per_pixel in [1.0, 2.0, 4.0]:
    out = rasterize(img, axes=("x", "y"), min_coordinate=[0, 0], max_coordinate=[100, 100], target_coordinate_system="global", target_unit_to_pixels=1 / units_per_pixel)
    o = np.asarray(out.data.compute())[0]
    ys, xs = np.nonzero(o)
    centre_px = pd.DataFrame({"x": xs.astype(float) + 0.5, "y": ys.astype(float) + 0.5})
    centre_units = transform(PointsModel.parse(centre_px, transformations={"global": out.attrs["transform"]["global"]}), to_coordinate_system="global").compute()
    expected_x = (xs[0] + 0.5) * units_per_pixel
    shift = centre_units.x.iloc[0] - expected_x
    ext = get_extent(out)["x"]
    print(f"output pixel size s={units_per_pixel}: bright output pixel x={xs[0]}, its centre maps to x={centre_units.x.iloc[0]:.2f} (expected {expected_x:.2f}), shift={shift:+.2f} units; extent x=({ext[0]:.1f}, {ext[1]:.1f}) (expected (0, 100))")
    bug |= abs(shift) > 1e-9
print("VERDICT:", "BUG REPRODUCED" if bug else "NOT REPRODUCED")

Observed output

output pixel size s=1.0: bright output pixel x=60, its centre maps to x=60.50 (expected 60.50), shift=+0.00 units; extent x=(0.0, 100.0) (expected (0, 100))
output pixel size s=2.0: bright output pixel x=30, its centre maps to x=60.50 (expected 61.00), shift=-0.50 units; extent x=(-0.5, 99.5) (expected (0, 100))
output pixel size s=4.0: bright output pixel x=15, its centre maps to x=60.50 (expected 62.00), shift=-1.50 units; extent x=(-1.5, 98.5) (expected (0, 100))
VERDICT: BUG REPRODUCED

Possible fix direction (unverified)

Use Sequence([scale, translation]) for the output transformation; consider order=1/coarsening for images. Related: #165, #166.

Environment

uv run repro.py with the PEP 723 metadata in the script (fresh, isolated environment; spatialdata built from main @ ccf1ea0 (2026-08-28); Python 3.13, latest releases of the dependencies at run time: pandas 3.0, anndata 0.13, zarr 3.3, dask 2026.8, numpy 2.5, geopandas 1.1, shapely 2.1). macOS (arm64). Also reproduced in a second environment with pandas 2.3.3 / anndata 0.12.11 / numpy 2.4.4 / zarr 3.2.1.

Possibly related issues

#165, #166

Automatically generated; discovered by an AI agent (Claude) and not yet reviewed by a human.