Bounding box query on points

spatialdata/_core/_spatial_query.py

@@ -0,0 +1,87 @@
+from __future__ import annotations
+
+from dataclasses import dataclass
+from typing import TYPE_CHECKING
+
+import numpy as np
+import pyarrow as pa
+
+from spatialdata._core.coordinate_system import CoordinateSystem, _get_spatial_axes
+
+if TYPE_CHECKING:
+ pass
+
+
+@dataclass(frozen=True)
+class BaseSpatialRequest:
+ """Base class for spatial queries."""
+
+ coordinate_system: CoordinateSystem
+
+ def __post_init__(self) -> None:
+ # validate the coordinate system
+ spatial_axes = _get_spatial_axes(self.coordinate_system)
+ if len(spatial_axes) == 0:
+ raise ValueError("No spatial axes in the requested coordinate system")
+
+
+@dataclass(frozen=True)
+class BoundingBoxRequest(BaseSpatialRequest):
+ """Query with an axis-aligned bounding box.
+
+ Attributes
+ ----------
+ coordinate_system : CoordinateSystem
+ The coordinate system the coordinates are expressed in.
+ min_coordinate : np.ndarray
+ The coordinate of the lower left hand corner (i.e., minimum values)
+ of the bounding box.
+ max_coordiate : np.ndarray
+ The coordinate of the upper right hand corner (i.e., maximum values)
+ of the bounding box
+ """
+
+ min_coordinate: np.ndarray # type: ignore[type-arg]
+ max_coordinate: np.ndarray # type: ignore[type-arg]
+
+
+def _bounding_box_query_points(points: pa.Table, request: BoundingBoxRequest) -> pa.Table:
+ """Perform a spatial bounding box query on a points element.
+
+ Parameters
+ ----------
+ points : pa.Table
+ The points element to perform the query on.
+ request : BoundingBoxRequest
+ The request for the query.
+
+ Returns
+ -------
+ query_result : pa.Table
+ The points contained within the specified bounding box.
+ """
+ spatial_axes = _get_spatial_axes(request.coordinate_system)
+
+ for axis_index, axis_name in enumerate(spatial_axes):
+ # filter by lower bound
+ min_value = request.min_coordinate[axis_index]
+ points = points.filter(pa.compute.greater(points[axis_name], min_value))
+
+ # filter by upper bound
+ max_value = request.max_coordinate[axis_index]
+ points = points.filter(pa.compute.less(points[axis_name], max_value))
+
+ return points
+
+
+def _bounding_box_query_points_dict(
+ points_dict: dict[str, pa.Table], request: BoundingBoxRequest
+) -> dict[str, pa.Table]:
+ requested_points = {}
+ for points_name, points_data in points_dict.items():
+ points = _bounding_box_query_points(points_data, request)
+ if len(points) > 0:
+ # do not include elements with no data
+ requested_points[points_name] = points
+
+ return requested_points

spatialdata/_core/_spatialdata.py

@@ -14,6 +14,11 @@
 from ome_zarr.types import JSONDict
 from spatial_image import SpatialImage
 
+from spatialdata._core._spatial_query import (
+ BaseSpatialRequest,
+ BoundingBoxRequest,
+ _bounding_box_query_points_dict,
+)
 from spatialdata._core.coordinate_system import CoordinateSystem
 from spatialdata._core.core_utils import SpatialElement, get_dims, get_transform
 from spatialdata._core.models import (
@@ -142,6 +147,12 @@ def __init__(
  Table_s.validate(table)
  self._table = table
 
+ self._query = QueryManager(self)
+
+ @property
+ def query(self) -> QueryManager:
+ return self._query
+
  def _add_image_in_memory(self, name: str, image: Union[SpatialImage, MultiscaleSpatialImage]) -> None:
  if name in self._images:
  raise ValueError(f"Image {name} already exists in the dataset.")
@@ -618,3 +629,34 @@ def _gen_elements(self) -> Generator[SpatialElement, None, None]:
  for element_type in ["images", "labels", "points", "polygons", "shapes"]:
  d = getattr(SpatialData, element_type).fget(self)
  yield from d.values()
+
+
+class QueryManager:
+ """Perform queries on SpatialData objects"""
+
+ def __init__(self, sdata: SpatialData):
+ self._sdata = sdata
+
+ def bounding_box(self, request: BoundingBoxRequest) -> SpatialData:
+ """Perform a bounding box query on the SpatialData object.
+
+ Parameters
+ ----------
+ request : BoundingBoxRequest
+ The bounding box request.
+
+ Returns
+ -------
+ requested_sdata : SpatialData
+ The SpatialData object containing the requested data.
+ Elements with no valid data are omitted.
+ """
+ requested_points = _bounding_box_query_points_dict(points_dict=self._sdata.points, request=request)
+
+ return SpatialData(points=requested_points)
+
+ def __call__(self, request: BaseSpatialRequest) -> SpatialData:
+ if isinstance(request, BoundingBoxRequest):
+ return self.bounding_box(request)
+ else:
+ raise TypeError("unknown request type")

spatialdata/_core/coordinate_system.py

@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 import json
 from typing import TYPE_CHECKING, Any, Optional, Union
 
@@ -43,7 +45,7 @@ def __repr__(self) -> str:
  return f"CoordinateSystem('{self.name}', {self._axes})"
 
  @staticmethod
- def from_dict(coord_sys: CoordSystem_t) -> "CoordinateSystem":
+ def from_dict(coord_sys: CoordSystem_t) -> CoordinateSystem:
  if "name" not in coord_sys.keys():
  raise ValueError("`coordinate_system` MUST have a name.")
  if "axes" not in coord_sys.keys():
@@ -81,7 +83,7 @@ def from_array(self, array: Any) -> None:
  raise NotImplementedError()
 
  @staticmethod
- def from_json(data: Union[str, bytes]) -> "CoordinateSystem":
+ def from_json(data: Union[str, bytes]) -> CoordinateSystem:
  coord_sys = json.loads(data)
  return CoordinateSystem.from_dict(coord_sys)
 
@@ -108,3 +110,21 @@ def axes_types(self) -> tuple[str, ...]:
 
  def __hash__(self) -> int:
  return hash(frozenset(self.to_dict()))
+
+
+def _get_spatial_axes(
+ coordinate_system: CoordinateSystem,
+) -> list[str]:
+ """Get the names of the spatial axes in a coordinate system.
+
+ Parameters
+ ----------
+ coordinate_system : CoordinateSystem
+ The coordinate system to get the spatial axes from.
+
+ Returns
+ -------
+ spatial_axis_names : List[str]
+ The names of the spatial axes.
+ """
+ return [axis.name for axis in coordinate_system._axes if axis.type == "space"]

tests/_core/test_spatial_query.py

@@ -0,0 +1,76 @@
+from dataclasses import FrozenInstanceError
+
+import numpy as np
+import pytest
+
+from spatialdata import PointsModel
+from spatialdata._core._spatial_query import (
+ BaseSpatialRequest,
+ BoundingBoxRequest,
+ _bounding_box_query_points,
+)
+from tests._core.conftest import c_cs, cyx_cs, czyx_cs, xy_cs
+
+
+def _make_points_element():
+ """Helper function to make a Points element."""
+ coordinates = np.array([[10, 10], [20, 20], [20, 30]], dtype=float)
+ return PointsModel.parse(coordinates)
+
+
+def test_bounding_box_request_immutable():
+ """Test that the bounding box request is immutable."""
+ request = BoundingBoxRequest(
+ coordinate_system=cyx_cs, min_coordinate=np.array([0, 0]), max_coordinate=np.array([10, 10])
+ )
+ isinstance(request, BaseSpatialRequest)
+
+ # fields should be immutable
+ with pytest.raises(FrozenInstanceError):
+ request.coordinate_system = czyx_cs
+ with pytest.raises(FrozenInstanceError):
+ request.min_coordinate = np.array([5, 5, 5])
+ with pytest.raises(FrozenInstanceError):
+ request.max_coordinate = np.array([5, 5, 5])
+
+
+def test_bounding_box_request_no_spatial_axes():
+ """Requests with no spatial axes should raise an error"""
+ with pytest.raises(ValueError):
+ _ = BoundingBoxRequest(coordinate_system=c_cs, min_coordinate=np.array([0]), max_coordinate=np.array([10]))
+
+
+def test_bounding_box_points():
+ """test the points bounding box_query"""
+ points_element = _make_points_element()
+ original_x = np.array(points_element["x"])
+ original_y = np.array(points_element["y"])
+
+ request = BoundingBoxRequest(
+ coordinate_system=xy_cs, min_coordinate=np.array([18, 25]), max_coordinate=np.array([22, 35])
+ )
+ points_result = _bounding_box_query_points(points_element, request)
+ np.testing.assert_allclose(points_result["x"], [20])
+ np.testing.assert_allclose(points_result["y"], [30])
+
+ # result should be valid points element
+ PointsModel.validate(points_result)
+
+ # original element should be unchanged
+ np.testing.assert_allclose(points_element["x"], original_x)
+ np.testing.assert_allclose(points_element["y"], original_y)
+
+
+def test_bounding_box_points_no_points():
+ """Points bounding box query with no points in range should
+ return a points element with length 0.
+ """
+ points_element = _make_points_element()
+ request = BoundingBoxRequest(
+ coordinate_system=xy_cs, min_coordinate=np.array([40, 50]), max_coordinate=np.array([45, 55])
+ )
+ points_result = _bounding_box_query_points(points_element, request)
+ assert len(points_result) == 0
+
+ # result should be valid points element
+ PointsModel.validate(points_result)

spatialdata/_core/_spatial_query.py

@@ -0,0 +1,87 @@
+from __future__ import annotations
+
+from dataclasses import dataclass
+from typing import TYPE_CHECKING
+
+import numpy as np
+import pyarrow as pa
+
+from spatialdata._core.coordinate_system import CoordinateSystem, _get_spatial_axes
+
+if TYPE_CHECKING:
+ pass
+
+
+@dataclass(frozen=True)
+class BaseSpatialRequest:
+ """Base class for spatial queries."""
+
+ coordinate_system: CoordinateSystem
+
+ def __post_init__(self) -> None:
+ # validate the coordinate system
+ spatial_axes = _get_spatial_axes(self.coordinate_system)
+ if len(spatial_axes) == 0:
+ raise ValueError("No spatial axes in the requested coordinate system")
+
+
+@dataclass(frozen=True)
+class BoundingBoxRequest(BaseSpatialRequest):
+ """Query with an axis-aligned bounding box.
+
+ Attributes
+ ----------
+ coordinate_system : CoordinateSystem
+ The coordinate system the coordinates are expressed in.
+ min_coordinate : np.ndarray
+ The coordinate of the lower left hand corner (i.e., minimum values)
+ of the bounding box.
+ max_coordiate : np.ndarray
+ The coordinate of the upper right hand corner (i.e., maximum values)
+ of the bounding box
+ """
+
+ min_coordinate: np.ndarray # type: ignore[type-arg]
+ max_coordinate: np.ndarray # type: ignore[type-arg]
+
+
+def _bounding_box_query_points(points: pa.Table, request: BoundingBoxRequest) -> pa.Table:
+ """Perform a spatial bounding box query on a points element.
+
+ Parameters
+ ----------
+ points : pa.Table
+ The points element to perform the query on.
+ request : BoundingBoxRequest
+ The request for the query.
+
+ Returns
+ -------
+ query_result : pa.Table
+ The points contained within the specified bounding box.
+ """
+ spatial_axes = _get_spatial_axes(request.coordinate_system)
+
+ for axis_index, axis_name in enumerate(spatial_axes):
+ # filter by lower bound
+ min_value = request.min_coordinate[axis_index]
+ points = points.filter(pa.compute.greater(points[axis_name], min_value))
+
+ # filter by upper bound
+ max_value = request.max_coordinate[axis_index]
+ points = points.filter(pa.compute.less(points[axis_name], max_value))
+
+ return points
+
+
+def _bounding_box_query_points_dict(
+ points_dict: dict[str, pa.Table], request: BoundingBoxRequest
+) -> dict[str, pa.Table]:
+ requested_points = {}
+ for points_name, points_data in points_dict.items():
+ points = _bounding_box_query_points(points_data, request)
+ if len(points) > 0:
+ # do not include elements with no data
+ requested_points[points_name] = points
+
+ return requested_points

spatialdata/_core/_spatialdata.py

@@ -14,6 +14,11 @@
 from ome_zarr.types import JSONDict
 from spatial_image import SpatialImage
 
+from spatialdata._core._spatial_query import (
+ BaseSpatialRequest,
+ BoundingBoxRequest,
+ _bounding_box_query_points_dict,
+)
 from spatialdata._core.coordinate_system import CoordinateSystem
 from spatialdata._core.core_utils import SpatialElement, get_dims, get_transform
 from spatialdata._core.models import (
@@ -142,6 +147,12 @@ def __init__(
  Table_s.validate(table)
  self._table = table
 
+ self._query = QueryManager(self)
+
+ @property
+ def query(self) -> QueryManager:
+ return self._query
+
  def _add_image_in_memory(self, name: str, image: Union[SpatialImage, MultiscaleSpatialImage]) -> None:
  if name in self._images:
  raise ValueError(f"Image {name} already exists in the dataset.")
@@ -618,3 +629,34 @@ def _gen_elements(self) -> Generator[SpatialElement, None, None]:
  for element_type in ["images", "labels", "points", "polygons", "shapes"]:
  d = getattr(SpatialData, element_type).fget(self)
  yield from d.values()
+
+
+class QueryManager:
+ """Perform queries on SpatialData objects"""
+
+ def __init__(self, sdata: SpatialData):
+ self._sdata = sdata
+
+ def bounding_box(self, request: BoundingBoxRequest) -> SpatialData:
+ """Perform a bounding box query on the SpatialData object.
+
+ Parameters
+ ----------
+ request : BoundingBoxRequest
+ The bounding box request.
+
+ Returns
+ -------
+ requested_sdata : SpatialData
+ The SpatialData object containing the requested data.
+ Elements with no valid data are omitted.
+ """
+ requested_points = _bounding_box_query_points_dict(points_dict=self._sdata.points, request=request)
+
+ return SpatialData(points=requested_points)
+
+ def __call__(self, request: BaseSpatialRequest) -> SpatialData:
+ if isinstance(request, BoundingBoxRequest):
+ return self.bounding_box(request)
+ else:
+ raise TypeError("unknown request type")

spatialdata/_core/coordinate_system.py

@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 import json
 from typing import TYPE_CHECKING, Any, Optional, Union
 
@@ -43,7 +45,7 @@ def __repr__(self) -> str:
  return f"CoordinateSystem('{self.name}', {self._axes})"
 
  @staticmethod
- def from_dict(coord_sys: CoordSystem_t) -> "CoordinateSystem":
+ def from_dict(coord_sys: CoordSystem_t) -> CoordinateSystem:
  if "name" not in coord_sys.keys():
  raise ValueError("`coordinate_system` MUST have a name.")
  if "axes" not in coord_sys.keys():
@@ -81,7 +83,7 @@ def from_array(self, array: Any) -> None:
  raise NotImplementedError()
 
  @staticmethod
- def from_json(data: Union[str, bytes]) -> "CoordinateSystem":
+ def from_json(data: Union[str, bytes]) -> CoordinateSystem:
  coord_sys = json.loads(data)
  return CoordinateSystem.from_dict(coord_sys)
 
@@ -108,3 +110,21 @@ def axes_types(self) -> tuple[str, ...]:
 
  def __hash__(self) -> int:
  return hash(frozenset(self.to_dict()))
+
+
+def _get_spatial_axes(
+ coordinate_system: CoordinateSystem,
+) -> list[str]:
+ """Get the names of the spatial axes in a coordinate system.
+
+ Parameters
+ ----------
+ coordinate_system : CoordinateSystem
+ The coordinate system to get the spatial axes from.
+
+ Returns
+ -------
+ spatial_axis_names : List[str]
+ The names of the spatial axes.
+ """
+ return [axis.name for axis in coordinate_system._axes if axis.type == "space"]

tests/_core/test_spatial_query.py

@@ -0,0 +1,76 @@
+from dataclasses import FrozenInstanceError
+
+import numpy as np
+import pytest
+
+from spatialdata import PointsModel
+from spatialdata._core._spatial_query import (
+ BaseSpatialRequest,
+ BoundingBoxRequest,
+ _bounding_box_query_points,
+)
+from tests._core.conftest import c_cs, cyx_cs, czyx_cs, xy_cs
+
+
+def _make_points_element():
+ """Helper function to make a Points element."""
+ coordinates = np.array([[10, 10], [20, 20], [20, 30]], dtype=float)
+ return PointsModel.parse(coordinates)
+
+
+def test_bounding_box_request_immutable():
+ """Test that the bounding box request is immutable."""
+ request = BoundingBoxRequest(
+ coordinate_system=cyx_cs, min_coordinate=np.array([0, 0]), max_coordinate=np.array([10, 10])
+ )
+ isinstance(request, BaseSpatialRequest)
+
+ # fields should be immutable
+ with pytest.raises(FrozenInstanceError):
+ request.coordinate_system = czyx_cs
+ with pytest.raises(FrozenInstanceError):
+ request.min_coordinate = np.array([5, 5, 5])
+ with pytest.raises(FrozenInstanceError):
+ request.max_coordinate = np.array([5, 5, 5])
+
+
+def test_bounding_box_request_no_spatial_axes():
+ """Requests with no spatial axes should raise an error"""
+ with pytest.raises(ValueError):
+ _ = BoundingBoxRequest(coordinate_system=c_cs, min_coordinate=np.array([0]), max_coordinate=np.array([10]))
+
+
+def test_bounding_box_points():
+ """test the points bounding box_query"""
+ points_element = _make_points_element()
+ original_x = np.array(points_element["x"])
+ original_y = np.array(points_element["y"])
+
+ request = BoundingBoxRequest(
+ coordinate_system=xy_cs, min_coordinate=np.array([18, 25]), max_coordinate=np.array([22, 35])
+ )
+ points_result = _bounding_box_query_points(points_element, request)
+ np.testing.assert_allclose(points_result["x"], [20])
+ np.testing.assert_allclose(points_result["y"], [30])
+
+ # result should be valid points element
+ PointsModel.validate(points_result)
+
+ # original element should be unchanged
+ np.testing.assert_allclose(points_element["x"], original_x)
+ np.testing.assert_allclose(points_element["y"], original_y)
+
+
+def test_bounding_box_points_no_points():
+ """Points bounding box query with no points in range should
+ return a points element with length 0.
+ """
+ points_element = _make_points_element()
+ request = BoundingBoxRequest(
+ coordinate_system=xy_cs, min_coordinate=np.array([40, 50]), max_coordinate=np.array([45, 55])
+ )
+ points_result = _bounding_box_query_points(points_element, request)
+ assert len(points_result) == 0
+
+ # result should be valid points element
+ PointsModel.validate(points_result)

spatialdata/_core/_spatial_query.py

@@ -0,0 +1,87 @@
+from __future__ import annotations
+
+from dataclasses import dataclass
+from typing import TYPE_CHECKING
+
+import numpy as np
+import pyarrow as pa
+
+from spatialdata._core.coordinate_system import CoordinateSystem, _get_spatial_axes
+
+if TYPE_CHECKING:
+ pass
+
+
+@dataclass(frozen=True)
+class BaseSpatialRequest:
+ """Base class for spatial queries."""
+
+ coordinate_system: CoordinateSystem
+
+ def __post_init__(self) -> None:
+ # validate the coordinate system
+ spatial_axes = _get_spatial_axes(self.coordinate_system)
+ if len(spatial_axes) == 0:
+ raise ValueError("No spatial axes in the requested coordinate system")
+
+
+@dataclass(frozen=True)
+class BoundingBoxRequest(BaseSpatialRequest):
+ """Query with an axis-aligned bounding box.
+
+ Attributes
+ ----------
+ coordinate_system : CoordinateSystem
+ The coordinate system the coordinates are expressed in.
+ min_coordinate : np.ndarray
+ The coordinate of the lower left hand corner (i.e., minimum values)
+ of the bounding box.
+ max_coordiate : np.ndarray
+ The coordinate of the upper right hand corner (i.e., maximum values)
+ of the bounding box
+ """
+
+ min_coordinate: np.ndarray # type: ignore[type-arg]
+ max_coordinate: np.ndarray # type: ignore[type-arg]
+
+
+def _bounding_box_query_points(points: pa.Table, request: BoundingBoxRequest) -> pa.Table:
+ """Perform a spatial bounding box query on a points element.
+
+ Parameters
+ ----------
+ points : pa.Table
+ The points element to perform the query on.
+ request : BoundingBoxRequest
+ The request for the query.
+
+ Returns
+ -------
+ query_result : pa.Table
+ The points contained within the specified bounding box.
+ """
+ spatial_axes = _get_spatial_axes(request.coordinate_system)
+
+ for axis_index, axis_name in enumerate(spatial_axes):
+ # filter by lower bound
+ min_value = request.min_coordinate[axis_index]
+ points = points.filter(pa.compute.greater(points[axis_name], min_value))
+
+ # filter by upper bound
+ max_value = request.max_coordinate[axis_index]
+ points = points.filter(pa.compute.less(points[axis_name], max_value))
+
+ return points
+
+
+def _bounding_box_query_points_dict(
+ points_dict: dict[str, pa.Table], request: BoundingBoxRequest
+) -> dict[str, pa.Table]:
+ requested_points = {}
+ for points_name, points_data in points_dict.items():
+ points = _bounding_box_query_points(points_data, request)
+ if len(points) > 0:
+ # do not include elements with no data
+ requested_points[points_name] = points
+
+ return requested_points

spatialdata/_core/_spatialdata.py

@@ -14,6 +14,11 @@
 from ome_zarr.types import JSONDict
 from spatial_image import SpatialImage
 
+from spatialdata._core._spatial_query import (
+ BaseSpatialRequest,
+ BoundingBoxRequest,
+ _bounding_box_query_points_dict,
+)
 from spatialdata._core.coordinate_system import CoordinateSystem
 from spatialdata._core.core_utils import SpatialElement, get_dims, get_transform
 from spatialdata._core.models import (
@@ -142,6 +147,12 @@ def __init__(
  Table_s.validate(table)
  self._table = table
 
+ self._query = QueryManager(self)
+
+ @property
+ def query(self) -> QueryManager:
+ return self._query
+
  def _add_image_in_memory(self, name: str, image: Union[SpatialImage, MultiscaleSpatialImage]) -> None:
  if name in self._images:
  raise ValueError(f"Image {name} already exists in the dataset.")
@@ -618,3 +629,34 @@ def _gen_elements(self) -> Generator[SpatialElement, None, None]:
  for element_type in ["images", "labels", "points", "polygons", "shapes"]:
  d = getattr(SpatialData, element_type).fget(self)
  yield from d.values()
+
+
+class QueryManager:
+ """Perform queries on SpatialData objects"""
+
+ def __init__(self, sdata: SpatialData):
+ self._sdata = sdata
+
+ def bounding_box(self, request: BoundingBoxRequest) -> SpatialData:
+ """Perform a bounding box query on the SpatialData object.
+
+ Parameters
+ ----------
+ request : BoundingBoxRequest
+ The bounding box request.
+
+ Returns
+ -------
+ requested_sdata : SpatialData
+ The SpatialData object containing the requested data.
+ Elements with no valid data are omitted.
+ """
+ requested_points = _bounding_box_query_points_dict(points_dict=self._sdata.points, request=request)
+
+ return SpatialData(points=requested_points)
+
+ def __call__(self, request: BaseSpatialRequest) -> SpatialData:
+ if isinstance(request, BoundingBoxRequest):
+ return self.bounding_box(request)
+ else:
+ raise TypeError("unknown request type")

spatialdata/_core/coordinate_system.py

@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 import json
 from typing import TYPE_CHECKING, Any, Optional, Union
 
@@ -43,7 +45,7 @@ def __repr__(self) -> str:
  return f"CoordinateSystem('{self.name}', {self._axes})"
 
  @staticmethod
- def from_dict(coord_sys: CoordSystem_t) -> "CoordinateSystem":
+ def from_dict(coord_sys: CoordSystem_t) -> CoordinateSystem:
  if "name" not in coord_sys.keys():
  raise ValueError("`coordinate_system` MUST have a name.")
  if "axes" not in coord_sys.keys():
@@ -81,7 +83,7 @@ def from_array(self, array: Any) -> None:
  raise NotImplementedError()
 
  @staticmethod
- def from_json(data: Union[str, bytes]) -> "CoordinateSystem":
+ def from_json(data: Union[str, bytes]) -> CoordinateSystem:
  coord_sys = json.loads(data)
  return CoordinateSystem.from_dict(coord_sys)
 
@@ -108,3 +110,21 @@ def axes_types(self) -> tuple[str, ...]:
 
  def __hash__(self) -> int:
  return hash(frozenset(self.to_dict()))
+
+
+def _get_spatial_axes(
+ coordinate_system: CoordinateSystem,
+) -> list[str]:
+ """Get the names of the spatial axes in a coordinate system.
+
+ Parameters
+ ----------
+ coordinate_system : CoordinateSystem
+ The coordinate system to get the spatial axes from.
+
+ Returns
+ -------
+ spatial_axis_names : List[str]
+ The names of the spatial axes.
+ """
+ return [axis.name for axis in coordinate_system._axes if axis.type == "space"]

tests/_core/test_spatial_query.py

@@ -0,0 +1,76 @@
+from dataclasses import FrozenInstanceError
+
+import numpy as np
+import pytest
+
+from spatialdata import PointsModel
+from spatialdata._core._spatial_query import (
+ BaseSpatialRequest,
+ BoundingBoxRequest,
+ _bounding_box_query_points,
+)
+from tests._core.conftest import c_cs, cyx_cs, czyx_cs, xy_cs
+
+
+def _make_points_element():
+ """Helper function to make a Points element."""
+ coordinates = np.array([[10, 10], [20, 20], [20, 30]], dtype=float)
+ return PointsModel.parse(coordinates)
+
+
+def test_bounding_box_request_immutable():
+ """Test that the bounding box request is immutable."""
+ request = BoundingBoxRequest(
+ coordinate_system=cyx_cs, min_coordinate=np.array([0, 0]), max_coordinate=np.array([10, 10])
+ )
+ isinstance(request, BaseSpatialRequest)
+
+ # fields should be immutable
+ with pytest.raises(FrozenInstanceError):
+ request.coordinate_system = czyx_cs
+ with pytest.raises(FrozenInstanceError):
+ request.min_coordinate = np.array([5, 5, 5])
+ with pytest.raises(FrozenInstanceError):
+ request.max_coordinate = np.array([5, 5, 5])
+
+
+def test_bounding_box_request_no_spatial_axes():
+ """Requests with no spatial axes should raise an error"""
+ with pytest.raises(ValueError):
+ _ = BoundingBoxRequest(coordinate_system=c_cs, min_coordinate=np.array([0]), max_coordinate=np.array([10]))
+
+
+def test_bounding_box_points():
+ """test the points bounding box_query"""
+ points_element = _make_points_element()
+ original_x = np.array(points_element["x"])
+ original_y = np.array(points_element["y"])
+
+ request = BoundingBoxRequest(
+ coordinate_system=xy_cs, min_coordinate=np.array([18, 25]), max_coordinate=np.array([22, 35])
+ )
+ points_result = _bounding_box_query_points(points_element, request)
+ np.testing.assert_allclose(points_result["x"], [20])
+ np.testing.assert_allclose(points_result["y"], [30])
+
+ # result should be valid points element
+ PointsModel.validate(points_result)
+
+ # original element should be unchanged
+ np.testing.assert_allclose(points_element["x"], original_x)
+ np.testing.assert_allclose(points_element["y"], original_y)
+
+
+def test_bounding_box_points_no_points():
+ """Points bounding box query with no points in range should
+ return a points element with length 0.
+ """
+ points_element = _make_points_element()
+ request = BoundingBoxRequest(
+ coordinate_system=xy_cs, min_coordinate=np.array([40, 50]), max_coordinate=np.array([45, 55])
+ )
+ points_result = _bounding_box_query_points(points_element, request)
+ assert len(points_result) == 0
+
+ # result should be valid points element
+ PointsModel.validate(points_result)

spatialdata/_core/_spatial_query.py

@@ -0,0 +1,87 @@
+from __future__ import annotations
+
+from dataclasses import dataclass
+from typing import TYPE_CHECKING
+
+import numpy as np
+import pyarrow as pa
+
+from spatialdata._core.coordinate_system import CoordinateSystem, _get_spatial_axes
+
+if TYPE_CHECKING:
+ pass
+
+
+@dataclass(frozen=True)
+class BaseSpatialRequest:
+ """Base class for spatial queries."""
+
+ coordinate_system: CoordinateSystem
+
+ def __post_init__(self) -> None:
+ # validate the coordinate system
+ spatial_axes = _get_spatial_axes(self.coordinate_system)
+ if len(spatial_axes) == 0:
+ raise ValueError("No spatial axes in the requested coordinate system")
+
+
+@dataclass(frozen=True)
+class BoundingBoxRequest(BaseSpatialRequest):
+ """Query with an axis-aligned bounding box.
+
+ Attributes
+ ----------
+ coordinate_system : CoordinateSystem
+ The coordinate system the coordinates are expressed in.
+ min_coordinate : np.ndarray
+ The coordinate of the lower left hand corner (i.e., minimum values)
+ of the bounding box.
+ max_coordiate : np.ndarray
+ The coordinate of the upper right hand corner (i.e., maximum values)
+ of the bounding box
+ """
+
+ min_coordinate: np.ndarray # type: ignore[type-arg]
+ max_coordinate: np.ndarray # type: ignore[type-arg]
+
+
+def _bounding_box_query_points(points: pa.Table, request: BoundingBoxRequest) -> pa.Table:
+ """Perform a spatial bounding box query on a points element.
+
+ Parameters
+ ----------
+ points : pa.Table
+ The points element to perform the query on.
+ request : BoundingBoxRequest
+ The request for the query.
+
+ Returns
+ -------
+ query_result : pa.Table
+ The points contained within the specified bounding box.
+ """
+ spatial_axes = _get_spatial_axes(request.coordinate_system)
+
+ for axis_index, axis_name in enumerate(spatial_axes):
+ # filter by lower bound
+ min_value = request.min_coordinate[axis_index]
+ points = points.filter(pa.compute.greater(points[axis_name], min_value))
+
+ # filter by upper bound
+ max_value = request.max_coordinate[axis_index]
+ points = points.filter(pa.compute.less(points[axis_name], max_value))
+
+ return points
+
+
+def _bounding_box_query_points_dict(
+ points_dict: dict[str, pa.Table], request: BoundingBoxRequest
+) -> dict[str, pa.Table]:
+ requested_points = {}
+ for points_name, points_data in points_dict.items():
+ points = _bounding_box_query_points(points_data, request)
+ if len(points) > 0:
+ # do not include elements with no data
+ requested_points[points_name] = points
+
+ return requested_points

spatialdata/_core/_spatialdata.py

@@ -14,6 +14,11 @@
 from ome_zarr.types import JSONDict
 from spatial_image import SpatialImage
 
+from spatialdata._core._spatial_query import (
+ BaseSpatialRequest,
+ BoundingBoxRequest,
+ _bounding_box_query_points_dict,
+)
 from spatialdata._core.coordinate_system import CoordinateSystem
 from spatialdata._core.core_utils import SpatialElement, get_dims, get_transform
 from spatialdata._core.models import (
@@ -142,6 +147,12 @@ def __init__(
  Table_s.validate(table)
  self._table = table
 
+ self._query = QueryManager(self)
+
+ @property
+ def query(self) -> QueryManager:
+ return self._query
+
  def _add_image_in_memory(self, name: str, image: Union[SpatialImage, MultiscaleSpatialImage]) -> None:
  if name in self._images:
  raise ValueError(f"Image {name} already exists in the dataset.")
@@ -618,3 +629,34 @@ def _gen_elements(self) -> Generator[SpatialElement, None, None]:
  for element_type in ["images", "labels", "points", "polygons", "shapes"]:
  d = getattr(SpatialData, element_type).fget(self)
  yield from d.values()
+
+
+class QueryManager:
+ """Perform queries on SpatialData objects"""
+
+ def __init__(self, sdata: SpatialData):
+ self._sdata = sdata
+
+ def bounding_box(self, request: BoundingBoxRequest) -> SpatialData:
+ """Perform a bounding box query on the SpatialData object.
+
+ Parameters
+ ----------
+ request : BoundingBoxRequest
+ The bounding box request.
+
+ Returns
+ -------
+ requested_sdata : SpatialData
+ The SpatialData object containing the requested data.
+ Elements with no valid data are omitted.
+ """
+ requested_points = _bounding_box_query_points_dict(points_dict=self._sdata.points, request=request)
+
+ return SpatialData(points=requested_points)
+
+ def __call__(self, request: BaseSpatialRequest) -> SpatialData:
+ if isinstance(request, BoundingBoxRequest):
+ return self.bounding_box(request)
+ else:
+ raise TypeError("unknown request type")

spatialdata/_core/coordinate_system.py

@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 import json
 from typing import TYPE_CHECKING, Any, Optional, Union
 
@@ -43,7 +45,7 @@ def __repr__(self) -> str:
  return f"CoordinateSystem('{self.name}', {self._axes})"
 
  @staticmethod
- def from_dict(coord_sys: CoordSystem_t) -> "CoordinateSystem":
+ def from_dict(coord_sys: CoordSystem_t) -> CoordinateSystem:
  if "name" not in coord_sys.keys():
  raise ValueError("`coordinate_system` MUST have a name.")
  if "axes" not in coord_sys.keys():
@@ -81,7 +83,7 @@ def from_array(self, array: Any) -> None:
  raise NotImplementedError()
 
  @staticmethod
- def from_json(data: Union[str, bytes]) -> "CoordinateSystem":
+ def from_json(data: Union[str, bytes]) -> CoordinateSystem:
  coord_sys = json.loads(data)
  return CoordinateSystem.from_dict(coord_sys)
 
@@ -108,3 +110,21 @@ def axes_types(self) -> tuple[str, ...]:
 
  def __hash__(self) -> int:
  return hash(frozenset(self.to_dict()))
+
+
+def _get_spatial_axes(
+ coordinate_system: CoordinateSystem,
+) -> list[str]:
+ """Get the names of the spatial axes in a coordinate system.
+
+ Parameters
+ ----------
+ coordinate_system : CoordinateSystem
+ The coordinate system to get the spatial axes from.
+
+ Returns
+ -------
+ spatial_axis_names : List[str]
+ The names of the spatial axes.
+ """
+ return [axis.name for axis in coordinate_system._axes if axis.type == "space"]

tests/_core/test_spatial_query.py

@@ -0,0 +1,76 @@
+from dataclasses import FrozenInstanceError
+
+import numpy as np
+import pytest
+
+from spatialdata import PointsModel
+from spatialdata._core._spatial_query import (
+ BaseSpatialRequest,
+ BoundingBoxRequest,
+ _bounding_box_query_points,
+)
+from tests._core.conftest import c_cs, cyx_cs, czyx_cs, xy_cs
+
+
+def _make_points_element():
+ """Helper function to make a Points element."""
+ coordinates = np.array([[10, 10], [20, 20], [20, 30]], dtype=float)
+ return PointsModel.parse(coordinates)
+
+
+def test_bounding_box_request_immutable():
+ """Test that the bounding box request is immutable."""
+ request = BoundingBoxRequest(
+ coordinate_system=cyx_cs, min_coordinate=np.array([0, 0]), max_coordinate=np.array([10, 10])
+ )
+ isinstance(request, BaseSpatialRequest)
+
+ # fields should be immutable
+ with pytest.raises(FrozenInstanceError):
+ request.coordinate_system = czyx_cs
+ with pytest.raises(FrozenInstanceError):
+ request.min_coordinate = np.array([5, 5, 5])
+ with pytest.raises(FrozenInstanceError):
+ request.max_coordinate = np.array([5, 5, 5])
+
+
+def test_bounding_box_request_no_spatial_axes():
+ """Requests with no spatial axes should raise an error"""
+ with pytest.raises(ValueError):
+ _ = BoundingBoxRequest(coordinate_system=c_cs, min_coordinate=np.array([0]), max_coordinate=np.array([10]))
+
+
+def test_bounding_box_points():
+ """test the points bounding box_query"""
+ points_element = _make_points_element()
+ original_x = np.array(points_element["x"])
+ original_y = np.array(points_element["y"])
+
+ request = BoundingBoxRequest(
+ coordinate_system=xy_cs, min_coordinate=np.array([18, 25]), max_coordinate=np.array([22, 35])
+ )
+ points_result = _bounding_box_query_points(points_element, request)
+ np.testing.assert_allclose(points_result["x"], [20])
+ np.testing.assert_allclose(points_result["y"], [30])
+
+ # result should be valid points element
+ PointsModel.validate(points_result)
+
+ # original element should be unchanged
+ np.testing.assert_allclose(points_element["x"], original_x)
+ np.testing.assert_allclose(points_element["y"], original_y)
+
+
+def test_bounding_box_points_no_points():
+ """Points bounding box query with no points in range should
+ return a points element with length 0.
+ """
+ points_element = _make_points_element()
+ request = BoundingBoxRequest(
+ coordinate_system=xy_cs, min_coordinate=np.array([40, 50]), max_coordinate=np.array([45, 55])
+ )
+ points_result = _bounding_box_query_points(points_element, request)
+ assert len(points_result) == 0
+
+ # result should be valid points element
+ PointsModel.validate(points_result)

spatialdata/_core/_spatial_query.py

@@ -0,0 +1,87 @@
+from __future__ import annotations
+
+from dataclasses import dataclass
+from typing import TYPE_CHECKING
+
+import numpy as np
+import pyarrow as pa
+
+from spatialdata._core.coordinate_system import CoordinateSystem, _get_spatial_axes
+
+if TYPE_CHECKING:
+ pass
+
+
+@dataclass(frozen=True)
+class BaseSpatialRequest:
+ """Base class for spatial queries."""
+
+ coordinate_system: CoordinateSystem
+
+ def __post_init__(self) -> None:
+ # validate the coordinate system
+ spatial_axes = _get_spatial_axes(self.coordinate_system)
+ if len(spatial_axes) == 0:
+ raise ValueError("No spatial axes in the requested coordinate system")
+
+
+@dataclass(frozen=True)
+class BoundingBoxRequest(BaseSpatialRequest):
+ """Query with an axis-aligned bounding box.
+
+ Attributes
+ ----------
+ coordinate_system : CoordinateSystem
+ The coordinate system the coordinates are expressed in.
+ min_coordinate : np.ndarray
+ The coordinate of the lower left hand corner (i.e., minimum values)
+ of the bounding box.
+ max_coordiate : np.ndarray
+ The coordinate of the upper right hand corner (i.e., maximum values)
+ of the bounding box
+ """
+
+ min_coordinate: np.ndarray # type: ignore[type-arg]
+ max_coordinate: np.ndarray # type: ignore[type-arg]
+
+
+def _bounding_box_query_points(points: pa.Table, request: BoundingBoxRequest) -> pa.Table:
+ """Perform a spatial bounding box query on a points element.
+
+ Parameters
+ ----------
+ points : pa.Table
+ The points element to perform the query on.
+ request : BoundingBoxRequest
+ The request for the query.
+
+ Returns
+ -------
+ query_result : pa.Table
+ The points contained within the specified bounding box.
+ """
+ spatial_axes = _get_spatial_axes(request.coordinate_system)
+
+ for axis_index, axis_name in enumerate(spatial_axes):
+ # filter by lower bound
+ min_value = request.min_coordinate[axis_index]
+ points = points.filter(pa.compute.greater(points[axis_name], min_value))
+
+ # filter by upper bound
+ max_value = request.max_coordinate[axis_index]
+ points = points.filter(pa.compute.less(points[axis_name], max_value))
+
+ return points
+
+
+def _bounding_box_query_points_dict(
+ points_dict: dict[str, pa.Table], request: BoundingBoxRequest
+) -> dict[str, pa.Table]:
+ requested_points = {}
+ for points_name, points_data in points_dict.items():
+ points = _bounding_box_query_points(points_data, request)
+ if len(points) > 0:
+ # do not include elements with no data
+ requested_points[points_name] = points
+
+ return requested_points

spatialdata/_core/_spatialdata.py

@@ -14,6 +14,11 @@
 from ome_zarr.types import JSONDict
 from spatial_image import SpatialImage
 
+from spatialdata._core._spatial_query import (
+ BaseSpatialRequest,
+ BoundingBoxRequest,
+ _bounding_box_query_points_dict,
+)
 from spatialdata._core.coordinate_system import CoordinateSystem
 from spatialdata._core.core_utils import SpatialElement, get_dims, get_transform
 from spatialdata._core.models import (
@@ -142,6 +147,12 @@ def __init__(
  Table_s.validate(table)
  self._table = table
 
+ self._query = QueryManager(self)
+
+ @property
+ def query(self) -> QueryManager:
+ return self._query
+
  def _add_image_in_memory(self, name: str, image: Union[SpatialImage, MultiscaleSpatialImage]) -> None:
  if name in self._images:
  raise ValueError(f"Image {name} already exists in the dataset.")
@@ -618,3 +629,34 @@ def _gen_elements(self) -> Generator[SpatialElement, None, None]:
  for element_type in ["images", "labels", "points", "polygons", "shapes"]:
  d = getattr(SpatialData, element_type).fget(self)
  yield from d.values()
+
+
+class QueryManager:
+ """Perform queries on SpatialData objects"""
+
+ def __init__(self, sdata: SpatialData):
+ self._sdata = sdata
+
+ def bounding_box(self, request: BoundingBoxRequest) -> SpatialData:
+ """Perform a bounding box query on the SpatialData object.
+
+ Parameters
+ ----------
+ request : BoundingBoxRequest
+ The bounding box request.
+
+ Returns
+ -------
+ requested_sdata : SpatialData
+ The SpatialData object containing the requested data.
+ Elements with no valid data are omitted.
+ """
+ requested_points = _bounding_box_query_points_dict(points_dict=self._sdata.points, request=request)
+
+ return SpatialData(points=requested_points)
+
+ def __call__(self, request: BaseSpatialRequest) -> SpatialData:
+ if isinstance(request, BoundingBoxRequest):
+ return self.bounding_box(request)
+ else:
+ raise TypeError("unknown request type")

spatialdata/_core/coordinate_system.py

@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 import json
 from typing import TYPE_CHECKING, Any, Optional, Union
 
@@ -43,7 +45,7 @@ def __repr__(self) -> str:
  return f"CoordinateSystem('{self.name}', {self._axes})"
 
  @staticmethod
- def from_dict(coord_sys: CoordSystem_t) -> "CoordinateSystem":
+ def from_dict(coord_sys: CoordSystem_t) -> CoordinateSystem:
  if "name" not in coord_sys.keys():
  raise ValueError("`coordinate_system` MUST have a name.")
  if "axes" not in coord_sys.keys():
@@ -81,7 +83,7 @@ def from_array(self, array: Any) -> None:
  raise NotImplementedError()
 
  @staticmethod
- def from_json(data: Union[str, bytes]) -> "CoordinateSystem":
+ def from_json(data: Union[str, bytes]) -> CoordinateSystem:
  coord_sys = json.loads(data)
  return CoordinateSystem.from_dict(coord_sys)
 
@@ -108,3 +110,21 @@ def axes_types(self) -> tuple[str, ...]:
 
  def __hash__(self) -> int:
  return hash(frozenset(self.to_dict()))
+
+
+def _get_spatial_axes(
+ coordinate_system: CoordinateSystem,
+) -> list[str]:
+ """Get the names of the spatial axes in a coordinate system.
+
+ Parameters
+ ----------
+ coordinate_system : CoordinateSystem
+ The coordinate system to get the spatial axes from.
+
+ Returns
+ -------
+ spatial_axis_names : List[str]
+ The names of the spatial axes.
+ """
+ return [axis.name for axis in coordinate_system._axes if axis.type == "space"]

tests/_core/test_spatial_query.py

@@ -0,0 +1,76 @@
+from dataclasses import FrozenInstanceError
+
+import numpy as np
+import pytest
+
+from spatialdata import PointsModel
+from spatialdata._core._spatial_query import (
+ BaseSpatialRequest,
+ BoundingBoxRequest,
+ _bounding_box_query_points,
+)
+from tests._core.conftest import c_cs, cyx_cs, czyx_cs, xy_cs
+
+
+def _make_points_element():
+ """Helper function to make a Points element."""
+ coordinates = np.array([[10, 10], [20, 20], [20, 30]], dtype=float)
+ return PointsModel.parse(coordinates)
+
+
+def test_bounding_box_request_immutable():
+ """Test that the bounding box request is immutable."""
+ request = BoundingBoxRequest(
+ coordinate_system=cyx_cs, min_coordinate=np.array([0, 0]), max_coordinate=np.array([10, 10])
+ )
+ isinstance(request, BaseSpatialRequest)
+
+ # fields should be immutable
+ with pytest.raises(FrozenInstanceError):
+ request.coordinate_system = czyx_cs
+ with pytest.raises(FrozenInstanceError):
+ request.min_coordinate = np.array([5, 5, 5])
+ with pytest.raises(FrozenInstanceError):
+ request.max_coordinate = np.array([5, 5, 5])
+
+
+def test_bounding_box_request_no_spatial_axes():
+ """Requests with no spatial axes should raise an error"""
+ with pytest.raises(ValueError):
+ _ = BoundingBoxRequest(coordinate_system=c_cs, min_coordinate=np.array([0]), max_coordinate=np.array([10]))
+
+
+def test_bounding_box_points():
+ """test the points bounding box_query"""
+ points_element = _make_points_element()
+ original_x = np.array(points_element["x"])
+ original_y = np.array(points_element["y"])
+
+ request = BoundingBoxRequest(
+ coordinate_system=xy_cs, min_coordinate=np.array([18, 25]), max_coordinate=np.array([22, 35])
+ )
+ points_result = _bounding_box_query_points(points_element, request)
+ np.testing.assert_allclose(points_result["x"], [20])
+ np.testing.assert_allclose(points_result["y"], [30])
+
+ # result should be valid points element
+ PointsModel.validate(points_result)
+
+ # original element should be unchanged
+ np.testing.assert_allclose(points_element["x"], original_x)
+ np.testing.assert_allclose(points_element["y"], original_y)
+
+
+def test_bounding_box_points_no_points():
+ """Points bounding box query with no points in range should
+ return a points element with length 0.
+ """
+ points_element = _make_points_element()
+ request = BoundingBoxRequest(
+ coordinate_system=xy_cs, min_coordinate=np.array([40, 50]), max_coordinate=np.array([45, 55])
+ )
+ points_result = _bounding_box_query_points(points_element, request)
+ assert len(points_result) == 0
+
+ # result should be valid points element
+ PointsModel.validate(points_result)

spatialdata/_core/_spatial_query.py

@@ -0,0 +1,87 @@
+from __future__ import annotations
+
+from dataclasses import dataclass
+from typing import TYPE_CHECKING
+
+import numpy as np
+import pyarrow as pa
+
+from spatialdata._core.coordinate_system import CoordinateSystem, _get_spatial_axes
+
+if TYPE_CHECKING:
+ pass
+
+
+@dataclass(frozen=True)
+class BaseSpatialRequest:
+ """Base class for spatial queries."""
+
+ coordinate_system: CoordinateSystem
+
+ def __post_init__(self) -> None:
+ # validate the coordinate system
+ spatial_axes = _get_spatial_axes(self.coordinate_system)
+ if len(spatial_axes) == 0:
+ raise ValueError("No spatial axes in the requested coordinate system")
+
+
+@dataclass(frozen=True)
+class BoundingBoxRequest(BaseSpatialRequest):
+ """Query with an axis-aligned bounding box.
+
+ Attributes
+ ----------
+ coordinate_system : CoordinateSystem
+ The coordinate system the coordinates are expressed in.
+ min_coordinate : np.ndarray
+ The coordinate of the lower left hand corner (i.e., minimum values)
+ of the bounding box.
+ max_coordiate : np.ndarray
+ The coordinate of the upper right hand corner (i.e., maximum values)
+ of the bounding box
+ """
+
+ min_coordinate: np.ndarray # type: ignore[type-arg]
+ max_coordinate: np.ndarray # type: ignore[type-arg]
+
+
+def _bounding_box_query_points(points: pa.Table, request: BoundingBoxRequest) -> pa.Table:
+ """Perform a spatial bounding box query on a points element.
+
+ Parameters
+ ----------
+ points : pa.Table
+ The points element to perform the query on.
+ request : BoundingBoxRequest
+ The request for the query.
+
+ Returns
+ -------
+ query_result : pa.Table
+ The points contained within the specified bounding box.
+ """
+ spatial_axes = _get_spatial_axes(request.coordinate_system)
+
+ for axis_index, axis_name in enumerate(spatial_axes):
+ # filter by lower bound
+ min_value = request.min_coordinate[axis_index]
+ points = points.filter(pa.compute.greater(points[axis_name], min_value))
+
+ # filter by upper bound
+ max_value = request.max_coordinate[axis_index]
+ points = points.filter(pa.compute.less(points[axis_name], max_value))
+
+ return points
+
+
+def _bounding_box_query_points_dict(
+ points_dict: dict[str, pa.Table], request: BoundingBoxRequest
+) -> dict[str, pa.Table]:
+ requested_points = {}
+ for points_name, points_data in points_dict.items():
+ points = _bounding_box_query_points(points_data, request)
+ if len(points) > 0:
+ # do not include elements with no data
+ requested_points[points_name] = points
+
+ return requested_points

spatialdata/_core/_spatialdata.py

@@ -14,6 +14,11 @@
 from ome_zarr.types import JSONDict
 from spatial_image import SpatialImage
 
+from spatialdata._core._spatial_query import (
+ BaseSpatialRequest,
+ BoundingBoxRequest,
+ _bounding_box_query_points_dict,
+)
 from spatialdata._core.coordinate_system import CoordinateSystem
 from spatialdata._core.core_utils import SpatialElement, get_dims, get_transform
 from spatialdata._core.models import (
@@ -142,6 +147,12 @@ def __init__(
  Table_s.validate(table)
  self._table = table
 
+ self._query = QueryManager(self)
+
+ @property
+ def query(self) -> QueryManager:
+ return self._query
+
  def _add_image_in_memory(self, name: str, image: Union[SpatialImage, MultiscaleSpatialImage]) -> None:
  if name in self._images:
  raise ValueError(f"Image {name} already exists in the dataset.")
@@ -618,3 +629,34 @@ def _gen_elements(self) -> Generator[SpatialElement, None, None]:
  for element_type in ["images", "labels", "points", "polygons", "shapes"]:
  d = getattr(SpatialData, element_type).fget(self)
  yield from d.values()
+
+
+class QueryManager:
+ """Perform queries on SpatialData objects"""
+
+ def __init__(self, sdata: SpatialData):
+ self._sdata = sdata
+
+ def bounding_box(self, request: BoundingBoxRequest) -> SpatialData:
+ """Perform a bounding box query on the SpatialData object.
+
+ Parameters
+ ----------
+ request : BoundingBoxRequest
+ The bounding box request.
+
+ Returns
+ -------
+ requested_sdata : SpatialData
+ The SpatialData object containing the requested data.
+ Elements with no valid data are omitted.
+ """
+ requested_points = _bounding_box_query_points_dict(points_dict=self._sdata.points, request=request)
+
+ return SpatialData(points=requested_points)
+
+ def __call__(self, request: BaseSpatialRequest) -> SpatialData:
+ if isinstance(request, BoundingBoxRequest):
+ return self.bounding_box(request)
+ else:
+ raise TypeError("unknown request type")

spatialdata/_core/coordinate_system.py

@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 import json
 from typing import TYPE_CHECKING, Any, Optional, Union
 
@@ -43,7 +45,7 @@ def __repr__(self) -> str:
  return f"CoordinateSystem('{self.name}', {self._axes})"
 
  @staticmethod
- def from_dict(coord_sys: CoordSystem_t) -> "CoordinateSystem":
+ def from_dict(coord_sys: CoordSystem_t) -> CoordinateSystem:
  if "name" not in coord_sys.keys():
  raise ValueError("`coordinate_system` MUST have a name.")
  if "axes" not in coord_sys.keys():
@@ -81,7 +83,7 @@ def from_array(self, array: Any) -> None:
  raise NotImplementedError()
 
  @staticmethod
- def from_json(data: Union[str, bytes]) -> "CoordinateSystem":
+ def from_json(data: Union[str, bytes]) -> CoordinateSystem:
  coord_sys = json.loads(data)
  return CoordinateSystem.from_dict(coord_sys)
 
@@ -108,3 +110,21 @@ def axes_types(self) -> tuple[str, ...]:
 
  def __hash__(self) -> int:
  return hash(frozenset(self.to_dict()))
+
+
+def _get_spatial_axes(
+ coordinate_system: CoordinateSystem,
+) -> list[str]:
+ """Get the names of the spatial axes in a coordinate system.
+
+ Parameters
+ ----------
+ coordinate_system : CoordinateSystem
+ The coordinate system to get the spatial axes from.
+
+ Returns
+ -------
+ spatial_axis_names : List[str]
+ The names of the spatial axes.
+ """
+ return [axis.name for axis in coordinate_system._axes if axis.type == "space"]

tests/_core/test_spatial_query.py

@@ -0,0 +1,76 @@
+from dataclasses import FrozenInstanceError
+
+import numpy as np
+import pytest
+
+from spatialdata import PointsModel
+from spatialdata._core._spatial_query import (
+ BaseSpatialRequest,
+ BoundingBoxRequest,
+ _bounding_box_query_points,
+)
+from tests._core.conftest import c_cs, cyx_cs, czyx_cs, xy_cs
+
+
+def _make_points_element():
+ """Helper function to make a Points element."""
+ coordinates = np.array([[10, 10], [20, 20], [20, 30]], dtype=float)
+ return PointsModel.parse(coordinates)
+
+
+def test_bounding_box_request_immutable():
+ """Test that the bounding box request is immutable."""
+ request = BoundingBoxRequest(
+ coordinate_system=cyx_cs, min_coordinate=np.array([0, 0]), max_coordinate=np.array([10, 10])
+ )
+ isinstance(request, BaseSpatialRequest)
+
+ # fields should be immutable
+ with pytest.raises(FrozenInstanceError):
+ request.coordinate_system = czyx_cs
+ with pytest.raises(FrozenInstanceError):
+ request.min_coordinate = np.array([5, 5, 5])
+ with pytest.raises(FrozenInstanceError):
+ request.max_coordinate = np.array([5, 5, 5])
+
+
+def test_bounding_box_request_no_spatial_axes():
+ """Requests with no spatial axes should raise an error"""
+ with pytest.raises(ValueError):
+ _ = BoundingBoxRequest(coordinate_system=c_cs, min_coordinate=np.array([0]), max_coordinate=np.array([10]))
+
+
+def test_bounding_box_points():
+ """test the points bounding box_query"""
+ points_element = _make_points_element()
+ original_x = np.array(points_element["x"])
+ original_y = np.array(points_element["y"])
+
+ request = BoundingBoxRequest(
+ coordinate_system=xy_cs, min_coordinate=np.array([18, 25]), max_coordinate=np.array([22, 35])
+ )
+ points_result = _bounding_box_query_points(points_element, request)
+ np.testing.assert_allclose(points_result["x"], [20])
+ np.testing.assert_allclose(points_result["y"], [30])
+
+ # result should be valid points element
+ PointsModel.validate(points_result)
+
+ # original element should be unchanged
+ np.testing.assert_allclose(points_element["x"], original_x)
+ np.testing.assert_allclose(points_element["y"], original_y)
+
+
+def test_bounding_box_points_no_points():
+ """Points bounding box query with no points in range should
+ return a points element with length 0.
+ """
+ points_element = _make_points_element()
+ request = BoundingBoxRequest(
+ coordinate_system=xy_cs, min_coordinate=np.array([40, 50]), max_coordinate=np.array([45, 55])
+ )
+ points_result = _bounding_box_query_points(points_element, request)
+ assert len(points_result) == 0
+
+ # result should be valid points element
+ PointsModel.validate(points_result)

spatialdata/_core/_spatial_query.py

@@ -0,0 +1,87 @@
+from __future__ import annotations
+
+from dataclasses import dataclass
+from typing import TYPE_CHECKING
+
+import numpy as np
+import pyarrow as pa
+
+from spatialdata._core.coordinate_system import CoordinateSystem, _get_spatial_axes
+
+if TYPE_CHECKING:
+ pass
+
+
+@dataclass(frozen=True)
+class BaseSpatialRequest:
+ """Base class for spatial queries."""
+
+ coordinate_system: CoordinateSystem
+
+ def __post_init__(self) -> None:
+ # validate the coordinate system
+ spatial_axes = _get_spatial_axes(self.coordinate_system)
+ if len(spatial_axes) == 0:
+ raise ValueError("No spatial axes in the requested coordinate system")
+
+
+@dataclass(frozen=True)
+class BoundingBoxRequest(BaseSpatialRequest):
+ """Query with an axis-aligned bounding box.
+
+ Attributes
+ ----------
+ coordinate_system : CoordinateSystem
+ The coordinate system the coordinates are expressed in.
+ min_coordinate : np.ndarray
+ The coordinate of the lower left hand corner (i.e., minimum values)
+ of the bounding box.
+ max_coordiate : np.ndarray
+ The coordinate of the upper right hand corner (i.e., maximum values)
+ of the bounding box
+ """
+
+ min_coordinate: np.ndarray # type: ignore[type-arg]
+ max_coordinate: np.ndarray # type: ignore[type-arg]
+
+
+def _bounding_box_query_points(points: pa.Table, request: BoundingBoxRequest) -> pa.Table:
+ """Perform a spatial bounding box query on a points element.
+
+ Parameters
+ ----------
+ points : pa.Table
+ The points element to perform the query on.
+ request : BoundingBoxRequest
+ The request for the query.
+
+ Returns
+ -------
+ query_result : pa.Table
+ The points contained within the specified bounding box.
+ """
+ spatial_axes = _get_spatial_axes(request.coordinate_system)
+
+ for axis_index, axis_name in enumerate(spatial_axes):
+ # filter by lower bound
+ min_value = request.min_coordinate[axis_index]
+ points = points.filter(pa.compute.greater(points[axis_name], min_value))
+
+ # filter by upper bound
+ max_value = request.max_coordinate[axis_index]
+ points = points.filter(pa.compute.less(points[axis_name], max_value))
+
+ return points
+
+
+def _bounding_box_query_points_dict(
+ points_dict: dict[str, pa.Table], request: BoundingBoxRequest
+) -> dict[str, pa.Table]:
+ requested_points = {}
+ for points_name, points_data in points_dict.items():
+ points = _bounding_box_query_points(points_data, request)
+ if len(points) > 0:
+ # do not include elements with no data
+ requested_points[points_name] = points
+
+ return requested_points

spatialdata/_core/_spatialdata.py

@@ -14,6 +14,11 @@
 from ome_zarr.types import JSONDict
 from spatial_image import SpatialImage
 
+from spatialdata._core._spatial_query import (
+ BaseSpatialRequest,
+ BoundingBoxRequest,
+ _bounding_box_query_points_dict,
+)
 from spatialdata._core.coordinate_system import CoordinateSystem
 from spatialdata._core.core_utils import SpatialElement, get_dims, get_transform
 from spatialdata._core.models import (
@@ -142,6 +147,12 @@ def __init__(
  Table_s.validate(table)
  self._table = table
 
+ self._query = QueryManager(self)
+
+ @property
+ def query(self) -> QueryManager:
+ return self._query
+
  def _add_image_in_memory(self, name: str, image: Union[SpatialImage, MultiscaleSpatialImage]) -> None:
  if name in self._images:
  raise ValueError(f"Image {name} already exists in the dataset.")
@@ -618,3 +629,34 @@ def _gen_elements(self) -> Generator[SpatialElement, None, None]:
  for element_type in ["images", "labels", "points", "polygons", "shapes"]:
  d = getattr(SpatialData, element_type).fget(self)
  yield from d.values()
+
+
+class QueryManager:
+ """Perform queries on SpatialData objects"""
+
+ def __init__(self, sdata: SpatialData):
+ self._sdata = sdata
+
+ def bounding_box(self, request: BoundingBoxRequest) -> SpatialData:
+ """Perform a bounding box query on the SpatialData object.
+
+ Parameters
+ ----------
+ request : BoundingBoxRequest
+ The bounding box request.
+
+ Returns
+ -------
+ requested_sdata : SpatialData
+ The SpatialData object containing the requested data.
+ Elements with no valid data are omitted.
+ """
+ requested_points = _bounding_box_query_points_dict(points_dict=self._sdata.points, request=request)
+
+ return SpatialData(points=requested_points)
+
+ def __call__(self, request: BaseSpatialRequest) -> SpatialData:
+ if isinstance(request, BoundingBoxRequest):
+ return self.bounding_box(request)
+ else:
+ raise TypeError("unknown request type")

spatialdata/_core/coordinate_system.py

@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 import json
 from typing import TYPE_CHECKING, Any, Optional, Union
 
@@ -43,7 +45,7 @@ def __repr__(self) -> str:
  return f"CoordinateSystem('{self.name}', {self._axes})"
 
  @staticmethod
- def from_dict(coord_sys: CoordSystem_t) -> "CoordinateSystem":
+ def from_dict(coord_sys: CoordSystem_t) -> CoordinateSystem:
  if "name" not in coord_sys.keys():
  raise ValueError("`coordinate_system` MUST have a name.")
  if "axes" not in coord_sys.keys():
@@ -81,7 +83,7 @@ def from_array(self, array: Any) -> None:
  raise NotImplementedError()
 
  @staticmethod
- def from_json(data: Union[str, bytes]) -> "CoordinateSystem":
+ def from_json(data: Union[str, bytes]) -> CoordinateSystem:
  coord_sys = json.loads(data)
  return CoordinateSystem.from_dict(coord_sys)
 
@@ -108,3 +110,21 @@ def axes_types(self) -> tuple[str, ...]:
 
  def __hash__(self) -> int:
  return hash(frozenset(self.to_dict()))
+
+
+def _get_spatial_axes(
+ coordinate_system: CoordinateSystem,
+) -> list[str]:
+ """Get the names of the spatial axes in a coordinate system.
+
+ Parameters
+ ----------
+ coordinate_system : CoordinateSystem
+ The coordinate system to get the spatial axes from.
+
+ Returns
+ -------
+ spatial_axis_names : List[str]
+ The names of the spatial axes.
+ """
+ return [axis.name for axis in coordinate_system._axes if axis.type == "space"]

tests/_core/test_spatial_query.py

@@ -0,0 +1,76 @@
+from dataclasses import FrozenInstanceError
+
+import numpy as np
+import pytest
+
+from spatialdata import PointsModel
+from spatialdata._core._spatial_query import (
+ BaseSpatialRequest,
+ BoundingBoxRequest,
+ _bounding_box_query_points,
+)
+from tests._core.conftest import c_cs, cyx_cs, czyx_cs, xy_cs
+
+
+def _make_points_element():
+ """Helper function to make a Points element."""
+ coordinates = np.array([[10, 10], [20, 20], [20, 30]], dtype=float)
+ return PointsModel.parse(coordinates)
+
+
+def test_bounding_box_request_immutable():
+ """Test that the bounding box request is immutable."""
+ request = BoundingBoxRequest(
+ coordinate_system=cyx_cs, min_coordinate=np.array([0, 0]), max_coordinate=np.array([10, 10])
+ )
+ isinstance(request, BaseSpatialRequest)
+
+ # fields should be immutable
+ with pytest.raises(FrozenInstanceError):
+ request.coordinate_system = czyx_cs
+ with pytest.raises(FrozenInstanceError):
+ request.min_coordinate = np.array([5, 5, 5])
+ with pytest.raises(FrozenInstanceError):
+ request.max_coordinate = np.array([5, 5, 5])
+
+
+def test_bounding_box_request_no_spatial_axes():
+ """Requests with no spatial axes should raise an error"""
+ with pytest.raises(ValueError):
+ _ = BoundingBoxRequest(coordinate_system=c_cs, min_coordinate=np.array([0]), max_coordinate=np.array([10]))
+
+
+def test_bounding_box_points():
+ """test the points bounding box_query"""
+ points_element = _make_points_element()
+ original_x = np.array(points_element["x"])
+ original_y = np.array(points_element["y"])
+
+ request = BoundingBoxRequest(
+ coordinate_system=xy_cs, min_coordinate=np.array([18, 25]), max_coordinate=np.array([22, 35])
+ )
+ points_result = _bounding_box_query_points(points_element, request)
+ np.testing.assert_allclose(points_result["x"], [20])
+ np.testing.assert_allclose(points_result["y"], [30])
+
+ # result should be valid points element
+ PointsModel.validate(points_result)
+
+ # original element should be unchanged
+ np.testing.assert_allclose(points_element["x"], original_x)
+ np.testing.assert_allclose(points_element["y"], original_y)
+
+
+def test_bounding_box_points_no_points():
+ """Points bounding box query with no points in range should
+ return a points element with length 0.
+ """
+ points_element = _make_points_element()
+ request = BoundingBoxRequest(
+ coordinate_system=xy_cs, min_coordinate=np.array([40, 50]), max_coordinate=np.array([45, 55])
+ )
+ points_result = _bounding_box_query_points(points_element, request)
+ assert len(points_result) == 0
+
+ # result should be valid points element
+ PointsModel.validate(points_result)

spatialdata/_core/_spatial_query.py

@@ -0,0 +1,87 @@
+from __future__ import annotations
+
+from dataclasses import dataclass
+from typing import TYPE_CHECKING
+
+import numpy as np
+import pyarrow as pa
+
+from spatialdata._core.coordinate_system import CoordinateSystem, _get_spatial_axes
+
+if TYPE_CHECKING:
+ pass
+
+
+@dataclass(frozen=True)
+class BaseSpatialRequest:
+ """Base class for spatial queries."""
+
+ coordinate_system: CoordinateSystem
+
+ def __post_init__(self) -> None:
+ # validate the coordinate system
+ spatial_axes = _get_spatial_axes(self.coordinate_system)
+ if len(spatial_axes) == 0:
+ raise ValueError("No spatial axes in the requested coordinate system")
+
+
+@dataclass(frozen=True)
+class BoundingBoxRequest(BaseSpatialRequest):
+ """Query with an axis-aligned bounding box.
+
+ Attributes
+ ----------
+ coordinate_system : CoordinateSystem
+ The coordinate system the coordinates are expressed in.
+ min_coordinate : np.ndarray
+ The coordinate of the lower left hand corner (i.e., minimum values)
+ of the bounding box.
+ max_coordiate : np.ndarray
+ The coordinate of the upper right hand corner (i.e., maximum values)
+ of the bounding box
+ """
+
+ min_coordinate: np.ndarray # type: ignore[type-arg]
+ max_coordinate: np.ndarray # type: ignore[type-arg]
+
+
+def _bounding_box_query_points(points: pa.Table, request: BoundingBoxRequest) -> pa.Table:
+ """Perform a spatial bounding box query on a points element.
+
+ Parameters
+ ----------
+ points : pa.Table
+ The points element to perform the query on.
+ request : BoundingBoxRequest
+ The request for the query.
+
+ Returns
+ -------
+ query_result : pa.Table
+ The points contained within the specified bounding box.
+ """
+ spatial_axes = _get_spatial_axes(request.coordinate_system)
+
+ for axis_index, axis_name in enumerate(spatial_axes):
+ # filter by lower bound
+ min_value = request.min_coordinate[axis_index]
+ points = points.filter(pa.compute.greater(points[axis_name], min_value))
+
+ # filter by upper bound
+ max_value = request.max_coordinate[axis_index]
+ points = points.filter(pa.compute.less(points[axis_name], max_value))
+
+ return points
+
+
+def _bounding_box_query_points_dict(
+ points_dict: dict[str, pa.Table], request: BoundingBoxRequest
+) -> dict[str, pa.Table]:
+ requested_points = {}
+ for points_name, points_data in points_dict.items():
+ points = _bounding_box_query_points(points_data, request)
+ if len(points) > 0:
+ # do not include elements with no data
+ requested_points[points_name] = points
+
+ return requested_points

spatialdata/_core/_spatialdata.py

@@ -14,6 +14,11 @@
 from ome_zarr.types import JSONDict
 from spatial_image import SpatialImage
 
+from spatialdata._core._spatial_query import (
+ BaseSpatialRequest,
+ BoundingBoxRequest,
+ _bounding_box_query_points_dict,
+)
 from spatialdata._core.coordinate_system import CoordinateSystem
 from spatialdata._core.core_utils import SpatialElement, get_dims, get_transform
 from spatialdata._core.models import (
@@ -142,6 +147,12 @@ def __init__(
  Table_s.validate(table)
  self._table = table
 
+ self._query = QueryManager(self)
+
+ @property
+ def query(self) -> QueryManager:
+ return self._query
+
  def _add_image_in_memory(self, name: str, image: Union[SpatialImage, MultiscaleSpatialImage]) -> None:
  if name in self._images:
  raise ValueError(f"Image {name} already exists in the dataset.")
@@ -618,3 +629,34 @@ def _gen_elements(self) -> Generator[SpatialElement, None, None]:
  for element_type in ["images", "labels", "points", "polygons", "shapes"]:
  d = getattr(SpatialData, element_type).fget(self)
  yield from d.values()
+
+
+class QueryManager:
+ """Perform queries on SpatialData objects"""
+
+ def __init__(self, sdata: SpatialData):
+ self._sdata = sdata
+
+ def bounding_box(self, request: BoundingBoxRequest) -> SpatialData:
+ """Perform a bounding box query on the SpatialData object.
+
+ Parameters
+ ----------
+ request : BoundingBoxRequest
+ The bounding box request.
+
+ Returns
+ -------
+ requested_sdata : SpatialData
+ The SpatialData object containing the requested data.
+ Elements with no valid data are omitted.
+ """
+ requested_points = _bounding_box_query_points_dict(points_dict=self._sdata.points, request=request)
+
+ return SpatialData(points=requested_points)
+
+ def __call__(self, request: BaseSpatialRequest) -> SpatialData:
+ if isinstance(request, BoundingBoxRequest):
+ return self.bounding_box(request)
+ else:
+ raise TypeError("unknown request type")

spatialdata/_core/coordinate_system.py

@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 import json
 from typing import TYPE_CHECKING, Any, Optional, Union
 
@@ -43,7 +45,7 @@ def __repr__(self) -> str:
  return f"CoordinateSystem('{self.name}', {self._axes})"
 
  @staticmethod
- def from_dict(coord_sys: CoordSystem_t) -> "CoordinateSystem":
+ def from_dict(coord_sys: CoordSystem_t) -> CoordinateSystem:
  if "name" not in coord_sys.keys():
  raise ValueError("`coordinate_system` MUST have a name.")
  if "axes" not in coord_sys.keys():
@@ -81,7 +83,7 @@ def from_array(self, array: Any) -> None:
  raise NotImplementedError()
 
  @staticmethod
- def from_json(data: Union[str, bytes]) -> "CoordinateSystem":
+ def from_json(data: Union[str, bytes]) -> CoordinateSystem:
  coord_sys = json.loads(data)
  return CoordinateSystem.from_dict(coord_sys)
 
@@ -108,3 +110,21 @@ def axes_types(self) -> tuple[str, ...]:
 
  def __hash__(self) -> int:
  return hash(frozenset(self.to_dict()))
+
+
+def _get_spatial_axes(
+ coordinate_system: CoordinateSystem,
+) -> list[str]:
+ """Get the names of the spatial axes in a coordinate system.
+
+ Parameters
+ ----------
+ coordinate_system : CoordinateSystem
+ The coordinate system to get the spatial axes from.
+
+ Returns
+ -------
+ spatial_axis_names : List[str]
+ The names of the spatial axes.
+ """
+ return [axis.name for axis in coordinate_system._axes if axis.type == "space"]

tests/_core/test_spatial_query.py

@@ -0,0 +1,76 @@
+from dataclasses import FrozenInstanceError
+
+import numpy as np
+import pytest
+
+from spatialdata import PointsModel
+from spatialdata._core._spatial_query import (
+ BaseSpatialRequest,
+ BoundingBoxRequest,
+ _bounding_box_query_points,
+)
+from tests._core.conftest import c_cs, cyx_cs, czyx_cs, xy_cs
+
+
+def _make_points_element():
+ """Helper function to make a Points element."""
+ coordinates = np.array([[10, 10], [20, 20], [20, 30]], dtype=float)
+ return PointsModel.parse(coordinates)
+
+
+def test_bounding_box_request_immutable():
+ """Test that the bounding box request is immutable."""
+ request = BoundingBoxRequest(
+ coordinate_system=cyx_cs, min_coordinate=np.array([0, 0]), max_coordinate=np.array([10, 10])
+ )
+ isinstance(request, BaseSpatialRequest)
+
+ # fields should be immutable
+ with pytest.raises(FrozenInstanceError):
+ request.coordinate_system = czyx_cs
+ with pytest.raises(FrozenInstanceError):
+ request.min_coordinate = np.array([5, 5, 5])
+ with pytest.raises(FrozenInstanceError):
+ request.max_coordinate = np.array([5, 5, 5])
+
+
+def test_bounding_box_request_no_spatial_axes():
+ """Requests with no spatial axes should raise an error"""
+ with pytest.raises(ValueError):
+ _ = BoundingBoxRequest(coordinate_system=c_cs, min_coordinate=np.array([0]), max_coordinate=np.array([10]))
+
+
+def test_bounding_box_points():
+ """test the points bounding box_query"""
+ points_element = _make_points_element()
+ original_x = np.array(points_element["x"])
+ original_y = np.array(points_element["y"])
+
+ request = BoundingBoxRequest(
+ coordinate_system=xy_cs, min_coordinate=np.array([18, 25]), max_coordinate=np.array([22, 35])
+ )
+ points_result = _bounding_box_query_points(points_element, request)
+ np.testing.assert_allclose(points_result["x"], [20])
+ np.testing.assert_allclose(points_result["y"], [30])
+
+ # result should be valid points element
+ PointsModel.validate(points_result)
+
+ # original element should be unchanged
+ np.testing.assert_allclose(points_element["x"], original_x)
+ np.testing.assert_allclose(points_element["y"], original_y)
+
+
+def test_bounding_box_points_no_points():
+ """Points bounding box query with no points in range should
+ return a points element with length 0.
+ """
+ points_element = _make_points_element()
+ request = BoundingBoxRequest(
+ coordinate_system=xy_cs, min_coordinate=np.array([40, 50]), max_coordinate=np.array([45, 55])
+ )
+ points_result = _bounding_box_query_points(points_element, request)
+ assert len(points_result) == 0
+
+ # result should be valid points element
+ PointsModel.validate(points_result)