Raster bounding_box_query/polygon_query: 1-pixel positioning error when the box minimum has a fractional part in (0, 0.5]
Nobody has claimed this yet.
- Dominant language
- Python
- Stars
- 394
- Forks
- 95
- Avg merge
- 4d 3h
- Merged PRs (30d)
- 7
Description
[!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 thespatialdatacode base. It has not been verified or triaged by a human yet; theneeds: triagelabel 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
Raster coordinates are pixel centres (0.5, 1.5, ...). image.sel(x=slice(2.3, 7.8)) keeps pixel 2 (centre 2.5 >= 2.3) as the first pixel, but the translation written into the returned transformation is ceil(2.3) = 3. The first kept pixel index is ceil(min - 0.5), not ceil(min), so for roughly half of all non-integer query boxes the returned raster is shifted by one pixel. Affects DataArray and DataTree inputs, and polygon_query on rasters (which uses the polygon bounds).
Severity (agent's assessment): high — silent mis-registration of the cropped raster by one full pixel relative to all other elements
Where: src/spatialdata/_core/query/_utils.py::_create_slices_and_translation (translation_vectors[i, j] = np.ceil(max(min_values[i, j], 0)))
Expected behaviour
The translation equals the index of the first selected pixel (e.g. 2 for xmin=2.3), so get_extent() of the crop matches the pixels that were kept.
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",
# ]
# ///
"""bounding_box_query on rasters: the translation is off by one pixel for fractional box minima in (k, k+0.5]."""
import warnings
import numpy as np
from spatialdata import bounding_box_query, get_extent
from spatialdata.models import Image2DModel
from spatialdata.transformations import get_transformation
warnings.simplefilter("ignore")
# pixel value encodes 10*y + x, so the first column value of row 0 == x index of the first selected pixel
img = Image2DModel.parse(np.arange(100, dtype=np.float64).reshape(1, 10, 10))
bug = False
for xmin in [2.0, 2.3, 2.5, 2.6, 3.0]:
res = bounding_box_query(img, axes=("x", "y"), min_coordinate=[xmin, 0], max_coordinate=[7.8, 10], target_coordinate_system="global")
first_pixel = int(res.data[0, 0, 0].compute())
tx = get_transformation(res, "global").to_affine_matrix(input_axes=("x", "y"), output_axes=("x", "y"))[0, 2]
ext = get_extent(res, coordinate_system="global")["x"]
ok = tx == first_pixel
print(f"xmin={xmin}: first selected pixel x-index={first_pixel}, translation x={tx}, extent x=({ext[0]:.1f}, {ext[1]:.1f}) -> {'OK' if ok else f'MISMATCH ({tx - first_pixel:+.0f} px)'}")
bug |= not ok
print("VERDICT:", "BUG REPRODUCED" if bug else "NOT REPRODUCED")
Observed output
xmin=2.0: first selected pixel x-index=2, translation x=2.0, extent x=(2.0, 8.0) -> OK
xmin=2.3: first selected pixel x-index=2, translation x=3.0, extent x=(3.0, 9.0) -> MISMATCH (+1 px)
xmin=2.5: first selected pixel x-index=2, translation x=3.0, extent x=(3.0, 9.0) -> MISMATCH (+1 px)
xmin=2.6: first selected pixel x-index=3, translation x=3.0, extent x=(3.0, 8.0) -> OK
xmin=3.0: first selected pixel x-index=3, translation x=3.0, extent x=(3.0, 8.0) -> OK
VERDICT: BUG REPRODUCED
Possible fix direction (unverified)
Derive the translation from the data actually selected, e.g. result.coords[axis].values[0] - 0.5 after image.sel(...) (or np.ceil(max(min - 0.5, 0))), and add a regression test with non-integer box bounds.
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.
Automatically generated; discovered by an AI agent (Claude) and not yet reviewed by a human.
Contributor guide
First steps
- Read the whole issue, then the project's contributing guide.
- Comment on the issue to say you are picking it up — it saves two people doing the same work.
- Fork the repository and make your change on a branch.
- Open a pull request that references the issue number.
Research direction
Start by running the provided repro.py with uv, then inspect src/spatialdata/_core/query/_utils.py, especially _create_slices_and_translation. Add a regression test covering fractional minimum bounds in (k, k+0.5] for raster bounding_box_query and polygon_query. Done means the translation equals the first selected pixel index and the returned extent matches the selected raster pixels.
Written by the indexing model from the issue text.
Assessment
- Tech stack
- python
- Domain
- data
- Issue type
- Bug
- Difficulty
- 3/5
- Estimated time
- 1-2 days
- Activity status
- Active
- Clarity
- Clearly specified
- Newbie friendliness
- 76/100