scverse / scverse/spatialdata

Raster bounding_box_query/polygon_query: 1-pixel positioning error when the box minimum has a fractional part in (0, 0.5]

Open
#1,203 0 comments 0 reactions 0 assignees View on GitHub

Nobody has claimed this yet.

bug 🚨 element: images 🌌 element: labels 🏷️ method: query needs: triage priority: high
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 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

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

Open the contributing guide

First steps

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. 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

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.