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

Inizia in src/spatialdata/_core/operations/rasterize.py, su rasterize_images_labels, e ispeziona la sequenza delle trasformazioni dell'output. Esegui il repro.py della issue con dimensioni dei pixel di output pari a 1, 2 e 4 per confrontare i centri e le estensioni dei pixel trasformati. Il lavoro è completato quando rasterize() preserva il bounding box richiesto e allinea i centri dei pixel di output per dimensioni dei pixel diverse da 1.

[!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.