Fixes with transformations; added transform_to_data_extent()

docs/api.md

@@ -32,6 +32,7 @@ Operations on `SpatialData` objects.
  transform
  rasterize
  aggregate
+ filter_by_coordinate_system
 ```
 
 ### Operations Utilities

docs/design_doc.md

@@ -564,6 +564,7 @@ with coordinate systems:
  with axes: c, y, x
  with elements: /images/point8, /labels/point8
 """
+
 sdata0 = sdata.query.coordinate_system("point23", filter_rows=False)
 sdata1 = sdata.query.bounding_box((0, 20, 0, 300))
 sdata1 = sdata.query.polygon("/polygons/annotations")

src/spatialdata/_core/data_extent.py

@@ -20,9 +20,6 @@
 from spatialdata.models._utils import SpatialElement
 from spatialdata.models.models import PointsModel
 from spatialdata.transformations.operations import get_transformation
-from spatialdata.transformations.transformations import (
- BaseTransformation,
-)
 
 BoundingBoxDescription = dict[str, tuple[float, float]]
 
@@ -289,9 +286,7 @@ def _(e: GeoDataFrame, coordinate_system: str = "global", exact: bool = True) ->
  coordinate_system=coordinate_system,
  extent=extent,
  )
- t = get_transformation(e, to_coordinate_system=coordinate_system)
- assert isinstance(t, BaseTransformation)
- transformed = transform(e, t)
+ transformed = transform(e, to_coordinate_system=coordinate_system)
  return _get_extent_of_shapes(transformed)
 
 
@@ -305,9 +300,7 @@ def _(e: DaskDataFrame, coordinate_system: str = "global", exact: bool = True) -
  coordinate_system=coordinate_system,
  extent=extent,
  )
- t = get_transformation(e, to_coordinate_system=coordinate_system)
- assert isinstance(t, BaseTransformation)
- transformed = transform(e, t)
+ transformed = transform(e, to_coordinate_system=coordinate_system)
  return _get_extent_of_points(transformed)
 
 
@@ -353,8 +346,6 @@ def _compute_extent_in_coordinate_system(
  -------
  The bounding box description in the specified coordinate system.
  """
- transformation = get_transformation(element, to_coordinate_system=coordinate_system)
- assert isinstance(transformation, BaseTransformation)
  from spatialdata._core.query._utils import get_bounding_box_corners
 
  axes = get_axes_names(element)
@@ -368,8 +359,10 @@ def _compute_extent_in_coordinate_system(
  max_coordinate=max_coordinates,
  )
  df = pd.DataFrame(corners.data, columns=corners.axis.data.tolist())
- points = PointsModel.parse(df, coordinates={k: k for k in axes})
- transformed_corners = pd.DataFrame(transform(points, transformation).compute())
+ d = get_transformation(element, get_all=True)
+ assert isinstance(d, dict)
+ points = PointsModel.parse(df, coordinates={k: k for k in axes}, transformations=d)
+ transformed_corners = pd.DataFrame(transform(points, to_coordinate_system=coordinate_system).compute())
  # Make sure min and max values are in the same order as axes
  extent = {}
  for ax in axes:

src/spatialdata/_core/operations/_utils.py

@@ -0,0 +1,117 @@
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+from multiscale_spatial_image import MultiscaleSpatialImage
+from spatial_image import SpatialImage
+
+if TYPE_CHECKING:
+ from spatialdata._core.spatialdata import SpatialData
+
+
+def transform_to_data_extent(
+ sdata: SpatialData,
+ coordinate_system: str,
+ target_unit_to_pixels: float | None = None,
+ target_width: float | None = None,
+ target_height: float | None = None,
+ target_depth: float | None = None,
+) -> SpatialData:
+ """
+ Transform the spatial data to match the data extent, keeping the positioning, and making all transformations equal.
+
+ Parameters
+ ----------
+ sdata
+ The spatial data to transform.
+ coordinate_system
+ The coordinate system to use to compute the extent and to transform the data to.
+ target_unit_to_pixels
+ The required number of pixels per unit (units in the target coordinate system) of the data that will be
+ produced.
+ target_width
+ The width of the data extent, in pixels, for the data that will be produced.
+ target_height
+ The height of the data extent, in pixels, for the data that will be produced.
+ target_depth
+ The depth of the data extent, in pixels, for the data that will be produced.
+
+ Returns
+ -------
+ SpatialData
+ The transformed spatial data with downscaled and padded images and adjusted vector coordinates; all the
+ transformations will set to Identity and the coordinates of the vector data will be aligned to the pixel
+ coordinates.
+
+ Notes
+ -----
+ - The data extent is the smallest rectangle that contains all the images and geometries.
+ - MultiscaleSpatialImage objects will be converted to SpatialImage objects.
+ - This helper function will be deprecated when https://github.com/scverse/spatialdata/issues/308 is closed,
+ as this function will be easily recovered by `transform_to_coordinate_system()`
+ """
+ # TODO: change "all the elements have idenity" with "all the elements have the same transformation"
+ from spatialdata._core.data_extent import get_extent
+ from spatialdata._core.operations.rasterize import _compute_target_dimensions, rasterize
+ from spatialdata._core.spatialdata import SpatialData
+ from spatialdata.transformations.operations import get_transformation, set_transformation
+ from spatialdata.transformations.transformations import BaseTransformation, Scale, Sequence, Translation
+
+ sdata = sdata.filter_by_coordinate_system(coordinate_system=coordinate_system)
+ # calling transform_to_coordinate_system will likely decrease the resolution, let's use rasterize() instead
+ sdata_vector = SpatialData(shapes=dict(sdata.shapes), points=dict(sdata.points))
+ sdata_raster = SpatialData(images=dict(sdata.images), labels=dict(sdata.labels))
+ sdata_vector_transformed = sdata_vector.transform_to_coordinate_system(coordinate_system)
+
+ de = get_extent(sdata, coordinate_system=coordinate_system)
+ de_axes = tuple(de.keys())
+ translation_to_origin = Translation([-de[ax][0] for ax in de_axes], axes=de_axes)
+
+ sizes = [de[ax][1] - de[ax][0] for ax in de_axes]
+ target_width, target_height, target_depth = _compute_target_dimensions(
+ spatial_axes=de_axes,
+ min_coordinate=[0 for _ in de_axes],
+ max_coordinate=sizes,
+ target_unit_to_pixels=target_unit_to_pixels,
+ target_width=target_width,
+ target_height=target_height,
+ target_depth=target_depth,
+ )
+ scale_to_target_d = {"x": target_width / sizes[de_axes.index("x")], "y": target_height / sizes[de_axes.index("y")]}
+ if target_depth is not None:
+ scale_to_target_d["z"] = target_depth / sizes[de_axes.index("z")]
+ scale_to_target = Scale([scale_to_target_d[ax] for ax in de_axes], axes=de_axes)
+
+ for el in sdata_vector_transformed._gen_elements_values():
+ t = get_transformation(el, to_coordinate_system=coordinate_system)
+ assert isinstance(t, BaseTransformation)
+ sequence = Sequence([t, translation_to_origin, scale_to_target])
+ set_transformation(el, transformation=sequence, to_coordinate_system=coordinate_system)
+ sdata_vector_transformed_inplace = sdata_vector_transformed.transform_to_coordinate_system(
+ coordinate_system, maintain_positioning=True
+ )
+
+ sdata_to_return_elements = {
+ **sdata_vector_transformed_inplace.shapes,
+ **sdata_vector_transformed_inplace.points,
+ }
+
+ for _, element_name, element in sdata_raster._gen_elements():
+ if isinstance(element, (MultiscaleSpatialImage, SpatialImage)):
+ rasterized = rasterize(
+ element,
+ axes=de_axes,
+ min_coordinate=[de[ax][0] for ax in de_axes],
+ max_coordinate=[de[ax][1] for ax in de_axes],
+ target_coordinate_system=coordinate_system,
+ target_unit_to_pixels=None,
+ target_width=target_width,
+ target_height=None,
+ target_depth=None,
+ )
+ sdata_to_return_elements[element_name] = rasterized
+ else:
+ sdata_to_return_elements[element_name] = element
+ if sdata.table is not None:
+ sdata_to_return_elements["table"] = sdata.table
+ return SpatialData.from_elements_dict(sdata_to_return_elements)

src/spatialdata/_core/operations/aggregate.py

@@ -170,8 +170,8 @@ def aggregate(
  target_coordinate_system, # type: ignore[assignment]
  )
  if not (by_transform == values_transform and isinstance(values_transform, Identity)):
- by_ = transform(by_, by_transform)
- values_ = transform(values_, values_transform)
+ by_ = transform(by_, to_coordinate_system=target_coordinate_system)
+ values_ = transform(values_, to_coordinate_system=target_coordinate_system)
 
  # dispatch
  adata = None

src/spatialdata/_core/operations/rasterize.py

@@ -448,7 +448,8 @@ def _(
  if schema in (Labels2DModel, Labels3DModel):
  kwargs = {"prefilter": False, "order": 0}
  elif schema in (Image2DModel, Image3DModel):
- kwargs = {}
+ kwargs = {"order": 0}
+ # kwargs = {}
  else:
  raise ValueError(f"Unsupported schema {schema}")
 

docs/api.md

@@ -32,6 +32,7 @@ Operations on `SpatialData` objects.
  transform
  rasterize
  aggregate
+ filter_by_coordinate_system
 ```
 
 ### Operations Utilities

docs/design_doc.md

@@ -564,6 +564,7 @@ with coordinate systems:
  with axes: c, y, x
  with elements: /images/point8, /labels/point8
 """
+
 sdata0 = sdata.query.coordinate_system("point23", filter_rows=False)
 sdata1 = sdata.query.bounding_box((0, 20, 0, 300))
 sdata1 = sdata.query.polygon("/polygons/annotations")

src/spatialdata/_core/data_extent.py

@@ -20,9 +20,6 @@
 from spatialdata.models._utils import SpatialElement
 from spatialdata.models.models import PointsModel
 from spatialdata.transformations.operations import get_transformation
-from spatialdata.transformations.transformations import (
- BaseTransformation,
-)
 
 BoundingBoxDescription = dict[str, tuple[float, float]]
 
@@ -289,9 +286,7 @@ def _(e: GeoDataFrame, coordinate_system: str = "global", exact: bool = True) ->
  coordinate_system=coordinate_system,
  extent=extent,
  )
- t = get_transformation(e, to_coordinate_system=coordinate_system)
- assert isinstance(t, BaseTransformation)
- transformed = transform(e, t)
+ transformed = transform(e, to_coordinate_system=coordinate_system)
  return _get_extent_of_shapes(transformed)
 
 
@@ -305,9 +300,7 @@ def _(e: DaskDataFrame, coordinate_system: str = "global", exact: bool = True) -
  coordinate_system=coordinate_system,
  extent=extent,
  )
- t = get_transformation(e, to_coordinate_system=coordinate_system)
- assert isinstance(t, BaseTransformation)
- transformed = transform(e, t)
+ transformed = transform(e, to_coordinate_system=coordinate_system)
  return _get_extent_of_points(transformed)
 
 
@@ -353,8 +346,6 @@ def _compute_extent_in_coordinate_system(
  -------
  The bounding box description in the specified coordinate system.
  """
- transformation = get_transformation(element, to_coordinate_system=coordinate_system)
- assert isinstance(transformation, BaseTransformation)
  from spatialdata._core.query._utils import get_bounding_box_corners
 
  axes = get_axes_names(element)
@@ -368,8 +359,10 @@ def _compute_extent_in_coordinate_system(
  max_coordinate=max_coordinates,
  )
  df = pd.DataFrame(corners.data, columns=corners.axis.data.tolist())
- points = PointsModel.parse(df, coordinates={k: k for k in axes})
- transformed_corners = pd.DataFrame(transform(points, transformation).compute())
+ d = get_transformation(element, get_all=True)
+ assert isinstance(d, dict)
+ points = PointsModel.parse(df, coordinates={k: k for k in axes}, transformations=d)
+ transformed_corners = pd.DataFrame(transform(points, to_coordinate_system=coordinate_system).compute())
  # Make sure min and max values are in the same order as axes
  extent = {}
  for ax in axes:

src/spatialdata/_core/operations/_utils.py

@@ -0,0 +1,117 @@
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+from multiscale_spatial_image import MultiscaleSpatialImage
+from spatial_image import SpatialImage
+
+if TYPE_CHECKING:
+ from spatialdata._core.spatialdata import SpatialData
+
+
+def transform_to_data_extent(
+ sdata: SpatialData,
+ coordinate_system: str,
+ target_unit_to_pixels: float | None = None,
+ target_width: float | None = None,
+ target_height: float | None = None,
+ target_depth: float | None = None,
+) -> SpatialData:
+ """
+ Transform the spatial data to match the data extent, keeping the positioning, and making all transformations equal.
+
+ Parameters
+ ----------
+ sdata
+ The spatial data to transform.
+ coordinate_system
+ The coordinate system to use to compute the extent and to transform the data to.
+ target_unit_to_pixels
+ The required number of pixels per unit (units in the target coordinate system) of the data that will be
+ produced.
+ target_width
+ The width of the data extent, in pixels, for the data that will be produced.
+ target_height
+ The height of the data extent, in pixels, for the data that will be produced.
+ target_depth
+ The depth of the data extent, in pixels, for the data that will be produced.
+
+ Returns
+ -------
+ SpatialData
+ The transformed spatial data with downscaled and padded images and adjusted vector coordinates; all the
+ transformations will set to Identity and the coordinates of the vector data will be aligned to the pixel
+ coordinates.
+
+ Notes
+ -----
+ - The data extent is the smallest rectangle that contains all the images and geometries.
+ - MultiscaleSpatialImage objects will be converted to SpatialImage objects.
+ - This helper function will be deprecated when https://github.com/scverse/spatialdata/issues/308 is closed,
+ as this function will be easily recovered by `transform_to_coordinate_system()`
+ """
+ # TODO: change "all the elements have idenity" with "all the elements have the same transformation"
+ from spatialdata._core.data_extent import get_extent
+ from spatialdata._core.operations.rasterize import _compute_target_dimensions, rasterize
+ from spatialdata._core.spatialdata import SpatialData
+ from spatialdata.transformations.operations import get_transformation, set_transformation
+ from spatialdata.transformations.transformations import BaseTransformation, Scale, Sequence, Translation
+
+ sdata = sdata.filter_by_coordinate_system(coordinate_system=coordinate_system)
+ # calling transform_to_coordinate_system will likely decrease the resolution, let's use rasterize() instead
+ sdata_vector = SpatialData(shapes=dict(sdata.shapes), points=dict(sdata.points))
+ sdata_raster = SpatialData(images=dict(sdata.images), labels=dict(sdata.labels))
+ sdata_vector_transformed = sdata_vector.transform_to_coordinate_system(coordinate_system)
+
+ de = get_extent(sdata, coordinate_system=coordinate_system)
+ de_axes = tuple(de.keys())
+ translation_to_origin = Translation([-de[ax][0] for ax in de_axes], axes=de_axes)
+
+ sizes = [de[ax][1] - de[ax][0] for ax in de_axes]
+ target_width, target_height, target_depth = _compute_target_dimensions(
+ spatial_axes=de_axes,
+ min_coordinate=[0 for _ in de_axes],
+ max_coordinate=sizes,
+ target_unit_to_pixels=target_unit_to_pixels,
+ target_width=target_width,
+ target_height=target_height,
+ target_depth=target_depth,
+ )
+ scale_to_target_d = {"x": target_width / sizes[de_axes.index("x")], "y": target_height / sizes[de_axes.index("y")]}
+ if target_depth is not None:
+ scale_to_target_d["z"] = target_depth / sizes[de_axes.index("z")]
+ scale_to_target = Scale([scale_to_target_d[ax] for ax in de_axes], axes=de_axes)
+
+ for el in sdata_vector_transformed._gen_elements_values():
+ t = get_transformation(el, to_coordinate_system=coordinate_system)
+ assert isinstance(t, BaseTransformation)
+ sequence = Sequence([t, translation_to_origin, scale_to_target])
+ set_transformation(el, transformation=sequence, to_coordinate_system=coordinate_system)
+ sdata_vector_transformed_inplace = sdata_vector_transformed.transform_to_coordinate_system(
+ coordinate_system, maintain_positioning=True
+ )
+
+ sdata_to_return_elements = {
+ **sdata_vector_transformed_inplace.shapes,
+ **sdata_vector_transformed_inplace.points,
+ }
+
+ for _, element_name, element in sdata_raster._gen_elements():
+ if isinstance(element, (MultiscaleSpatialImage, SpatialImage)):
+ rasterized = rasterize(
+ element,
+ axes=de_axes,
+ min_coordinate=[de[ax][0] for ax in de_axes],
+ max_coordinate=[de[ax][1] for ax in de_axes],
+ target_coordinate_system=coordinate_system,
+ target_unit_to_pixels=None,
+ target_width=target_width,
+ target_height=None,
+ target_depth=None,
+ )
+ sdata_to_return_elements[element_name] = rasterized
+ else:
+ sdata_to_return_elements[element_name] = element
+ if sdata.table is not None:
+ sdata_to_return_elements["table"] = sdata.table
+ return SpatialData.from_elements_dict(sdata_to_return_elements)

src/spatialdata/_core/operations/aggregate.py

@@ -170,8 +170,8 @@ def aggregate(
  target_coordinate_system, # type: ignore[assignment]
  )
  if not (by_transform == values_transform and isinstance(values_transform, Identity)):
- by_ = transform(by_, by_transform)
- values_ = transform(values_, values_transform)
+ by_ = transform(by_, to_coordinate_system=target_coordinate_system)
+ values_ = transform(values_, to_coordinate_system=target_coordinate_system)
 
  # dispatch
  adata = None

src/spatialdata/_core/operations/rasterize.py

@@ -448,7 +448,8 @@ def _(
  if schema in (Labels2DModel, Labels3DModel):
  kwargs = {"prefilter": False, "order": 0}
  elif schema in (Image2DModel, Image3DModel):
- kwargs = {}
+ kwargs = {"order": 0}
+ # kwargs = {}
  else:
  raise ValueError(f"Unsupported schema {schema}")
 

docs/api.md

@@ -32,6 +32,7 @@ Operations on `SpatialData` objects.
  transform
  rasterize
  aggregate
+ filter_by_coordinate_system
 ```
 
 ### Operations Utilities

docs/design_doc.md

@@ -564,6 +564,7 @@ with coordinate systems:
  with axes: c, y, x
  with elements: /images/point8, /labels/point8
 """
+
 sdata0 = sdata.query.coordinate_system("point23", filter_rows=False)
 sdata1 = sdata.query.bounding_box((0, 20, 0, 300))
 sdata1 = sdata.query.polygon("/polygons/annotations")

src/spatialdata/_core/data_extent.py

@@ -20,9 +20,6 @@
 from spatialdata.models._utils import SpatialElement
 from spatialdata.models.models import PointsModel
 from spatialdata.transformations.operations import get_transformation
-from spatialdata.transformations.transformations import (
- BaseTransformation,
-)
 
 BoundingBoxDescription = dict[str, tuple[float, float]]
 
@@ -289,9 +286,7 @@ def _(e: GeoDataFrame, coordinate_system: str = "global", exact: bool = True) ->
  coordinate_system=coordinate_system,
  extent=extent,
  )
- t = get_transformation(e, to_coordinate_system=coordinate_system)
- assert isinstance(t, BaseTransformation)
- transformed = transform(e, t)
+ transformed = transform(e, to_coordinate_system=coordinate_system)
  return _get_extent_of_shapes(transformed)
 
 
@@ -305,9 +300,7 @@ def _(e: DaskDataFrame, coordinate_system: str = "global", exact: bool = True) -
  coordinate_system=coordinate_system,
  extent=extent,
  )
- t = get_transformation(e, to_coordinate_system=coordinate_system)
- assert isinstance(t, BaseTransformation)
- transformed = transform(e, t)
+ transformed = transform(e, to_coordinate_system=coordinate_system)
  return _get_extent_of_points(transformed)
 
 
@@ -353,8 +346,6 @@ def _compute_extent_in_coordinate_system(
  -------
  The bounding box description in the specified coordinate system.
  """
- transformation = get_transformation(element, to_coordinate_system=coordinate_system)
- assert isinstance(transformation, BaseTransformation)
  from spatialdata._core.query._utils import get_bounding_box_corners
 
  axes = get_axes_names(element)
@@ -368,8 +359,10 @@ def _compute_extent_in_coordinate_system(
  max_coordinate=max_coordinates,
  )
  df = pd.DataFrame(corners.data, columns=corners.axis.data.tolist())
- points = PointsModel.parse(df, coordinates={k: k for k in axes})
- transformed_corners = pd.DataFrame(transform(points, transformation).compute())
+ d = get_transformation(element, get_all=True)
+ assert isinstance(d, dict)
+ points = PointsModel.parse(df, coordinates={k: k for k in axes}, transformations=d)
+ transformed_corners = pd.DataFrame(transform(points, to_coordinate_system=coordinate_system).compute())
  # Make sure min and max values are in the same order as axes
  extent = {}
  for ax in axes:

src/spatialdata/_core/operations/_utils.py

@@ -0,0 +1,117 @@
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+from multiscale_spatial_image import MultiscaleSpatialImage
+from spatial_image import SpatialImage
+
+if TYPE_CHECKING:
+ from spatialdata._core.spatialdata import SpatialData
+
+
+def transform_to_data_extent(
+ sdata: SpatialData,
+ coordinate_system: str,
+ target_unit_to_pixels: float | None = None,
+ target_width: float | None = None,
+ target_height: float | None = None,
+ target_depth: float | None = None,
+) -> SpatialData:
+ """
+ Transform the spatial data to match the data extent, keeping the positioning, and making all transformations equal.
+
+ Parameters
+ ----------
+ sdata
+ The spatial data to transform.
+ coordinate_system
+ The coordinate system to use to compute the extent and to transform the data to.
+ target_unit_to_pixels
+ The required number of pixels per unit (units in the target coordinate system) of the data that will be
+ produced.
+ target_width
+ The width of the data extent, in pixels, for the data that will be produced.
+ target_height
+ The height of the data extent, in pixels, for the data that will be produced.
+ target_depth
+ The depth of the data extent, in pixels, for the data that will be produced.
+
+ Returns
+ -------
+ SpatialData
+ The transformed spatial data with downscaled and padded images and adjusted vector coordinates; all the
+ transformations will set to Identity and the coordinates of the vector data will be aligned to the pixel
+ coordinates.
+
+ Notes
+ -----
+ - The data extent is the smallest rectangle that contains all the images and geometries.
+ - MultiscaleSpatialImage objects will be converted to SpatialImage objects.
+ - This helper function will be deprecated when https://github.com/scverse/spatialdata/issues/308 is closed,
+ as this function will be easily recovered by `transform_to_coordinate_system()`
+ """
+ # TODO: change "all the elements have idenity" with "all the elements have the same transformation"
+ from spatialdata._core.data_extent import get_extent
+ from spatialdata._core.operations.rasterize import _compute_target_dimensions, rasterize
+ from spatialdata._core.spatialdata import SpatialData
+ from spatialdata.transformations.operations import get_transformation, set_transformation
+ from spatialdata.transformations.transformations import BaseTransformation, Scale, Sequence, Translation
+
+ sdata = sdata.filter_by_coordinate_system(coordinate_system=coordinate_system)
+ # calling transform_to_coordinate_system will likely decrease the resolution, let's use rasterize() instead
+ sdata_vector = SpatialData(shapes=dict(sdata.shapes), points=dict(sdata.points))
+ sdata_raster = SpatialData(images=dict(sdata.images), labels=dict(sdata.labels))
+ sdata_vector_transformed = sdata_vector.transform_to_coordinate_system(coordinate_system)
+
+ de = get_extent(sdata, coordinate_system=coordinate_system)
+ de_axes = tuple(de.keys())
+ translation_to_origin = Translation([-de[ax][0] for ax in de_axes], axes=de_axes)
+
+ sizes = [de[ax][1] - de[ax][0] for ax in de_axes]
+ target_width, target_height, target_depth = _compute_target_dimensions(
+ spatial_axes=de_axes,
+ min_coordinate=[0 for _ in de_axes],
+ max_coordinate=sizes,
+ target_unit_to_pixels=target_unit_to_pixels,
+ target_width=target_width,
+ target_height=target_height,
+ target_depth=target_depth,
+ )
+ scale_to_target_d = {"x": target_width / sizes[de_axes.index("x")], "y": target_height / sizes[de_axes.index("y")]}
+ if target_depth is not None:
+ scale_to_target_d["z"] = target_depth / sizes[de_axes.index("z")]
+ scale_to_target = Scale([scale_to_target_d[ax] for ax in de_axes], axes=de_axes)
+
+ for el in sdata_vector_transformed._gen_elements_values():
+ t = get_transformation(el, to_coordinate_system=coordinate_system)
+ assert isinstance(t, BaseTransformation)
+ sequence = Sequence([t, translation_to_origin, scale_to_target])
+ set_transformation(el, transformation=sequence, to_coordinate_system=coordinate_system)
+ sdata_vector_transformed_inplace = sdata_vector_transformed.transform_to_coordinate_system(
+ coordinate_system, maintain_positioning=True
+ )
+
+ sdata_to_return_elements = {
+ **sdata_vector_transformed_inplace.shapes,
+ **sdata_vector_transformed_inplace.points,
+ }
+
+ for _, element_name, element in sdata_raster._gen_elements():
+ if isinstance(element, (MultiscaleSpatialImage, SpatialImage)):
+ rasterized = rasterize(
+ element,
+ axes=de_axes,
+ min_coordinate=[de[ax][0] for ax in de_axes],
+ max_coordinate=[de[ax][1] for ax in de_axes],
+ target_coordinate_system=coordinate_system,
+ target_unit_to_pixels=None,
+ target_width=target_width,
+ target_height=None,
+ target_depth=None,
+ )
+ sdata_to_return_elements[element_name] = rasterized
+ else:
+ sdata_to_return_elements[element_name] = element
+ if sdata.table is not None:
+ sdata_to_return_elements["table"] = sdata.table
+ return SpatialData.from_elements_dict(sdata_to_return_elements)

src/spatialdata/_core/operations/aggregate.py

@@ -170,8 +170,8 @@ def aggregate(
  target_coordinate_system, # type: ignore[assignment]
  )
  if not (by_transform == values_transform and isinstance(values_transform, Identity)):
- by_ = transform(by_, by_transform)
- values_ = transform(values_, values_transform)
+ by_ = transform(by_, to_coordinate_system=target_coordinate_system)
+ values_ = transform(values_, to_coordinate_system=target_coordinate_system)
 
  # dispatch
  adata = None

src/spatialdata/_core/operations/rasterize.py

@@ -448,7 +448,8 @@ def _(
  if schema in (Labels2DModel, Labels3DModel):
  kwargs = {"prefilter": False, "order": 0}
  elif schema in (Image2DModel, Image3DModel):
- kwargs = {}
+ kwargs = {"order": 0}
+ # kwargs = {}
  else:
  raise ValueError(f"Unsupported schema {schema}")
 

docs/api.md

@@ -32,6 +32,7 @@ Operations on `SpatialData` objects.
  transform
  rasterize
  aggregate
+ filter_by_coordinate_system
 ```
 
 ### Operations Utilities

docs/design_doc.md

@@ -564,6 +564,7 @@ with coordinate systems:
  with axes: c, y, x
  with elements: /images/point8, /labels/point8
 """
+
 sdata0 = sdata.query.coordinate_system("point23", filter_rows=False)
 sdata1 = sdata.query.bounding_box((0, 20, 0, 300))
 sdata1 = sdata.query.polygon("/polygons/annotations")

src/spatialdata/_core/data_extent.py

@@ -20,9 +20,6 @@
 from spatialdata.models._utils import SpatialElement
 from spatialdata.models.models import PointsModel
 from spatialdata.transformations.operations import get_transformation
-from spatialdata.transformations.transformations import (
- BaseTransformation,
-)
 
 BoundingBoxDescription = dict[str, tuple[float, float]]
 
@@ -289,9 +286,7 @@ def _(e: GeoDataFrame, coordinate_system: str = "global", exact: bool = True) ->
  coordinate_system=coordinate_system,
  extent=extent,
  )
- t = get_transformation(e, to_coordinate_system=coordinate_system)
- assert isinstance(t, BaseTransformation)
- transformed = transform(e, t)
+ transformed = transform(e, to_coordinate_system=coordinate_system)
  return _get_extent_of_shapes(transformed)
 
 
@@ -305,9 +300,7 @@ def _(e: DaskDataFrame, coordinate_system: str = "global", exact: bool = True) -
  coordinate_system=coordinate_system,
  extent=extent,
  )
- t = get_transformation(e, to_coordinate_system=coordinate_system)
- assert isinstance(t, BaseTransformation)
- transformed = transform(e, t)
+ transformed = transform(e, to_coordinate_system=coordinate_system)
  return _get_extent_of_points(transformed)
 
 
@@ -353,8 +346,6 @@ def _compute_extent_in_coordinate_system(
  -------
  The bounding box description in the specified coordinate system.
  """
- transformation = get_transformation(element, to_coordinate_system=coordinate_system)
- assert isinstance(transformation, BaseTransformation)
  from spatialdata._core.query._utils import get_bounding_box_corners
 
  axes = get_axes_names(element)
@@ -368,8 +359,10 @@ def _compute_extent_in_coordinate_system(
  max_coordinate=max_coordinates,
  )
  df = pd.DataFrame(corners.data, columns=corners.axis.data.tolist())
- points = PointsModel.parse(df, coordinates={k: k for k in axes})
- transformed_corners = pd.DataFrame(transform(points, transformation).compute())
+ d = get_transformation(element, get_all=True)
+ assert isinstance(d, dict)
+ points = PointsModel.parse(df, coordinates={k: k for k in axes}, transformations=d)
+ transformed_corners = pd.DataFrame(transform(points, to_coordinate_system=coordinate_system).compute())
  # Make sure min and max values are in the same order as axes
  extent = {}
  for ax in axes:

src/spatialdata/_core/operations/_utils.py

@@ -0,0 +1,117 @@
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+from multiscale_spatial_image import MultiscaleSpatialImage
+from spatial_image import SpatialImage
+
+if TYPE_CHECKING:
+ from spatialdata._core.spatialdata import SpatialData
+
+
+def transform_to_data_extent(
+ sdata: SpatialData,
+ coordinate_system: str,
+ target_unit_to_pixels: float | None = None,
+ target_width: float | None = None,
+ target_height: float | None = None,
+ target_depth: float | None = None,
+) -> SpatialData:
+ """
+ Transform the spatial data to match the data extent, keeping the positioning, and making all transformations equal.
+
+ Parameters
+ ----------
+ sdata
+ The spatial data to transform.
+ coordinate_system
+ The coordinate system to use to compute the extent and to transform the data to.
+ target_unit_to_pixels
+ The required number of pixels per unit (units in the target coordinate system) of the data that will be
+ produced.
+ target_width
+ The width of the data extent, in pixels, for the data that will be produced.
+ target_height
+ The height of the data extent, in pixels, for the data that will be produced.
+ target_depth
+ The depth of the data extent, in pixels, for the data that will be produced.
+
+ Returns
+ -------
+ SpatialData
+ The transformed spatial data with downscaled and padded images and adjusted vector coordinates; all the
+ transformations will set to Identity and the coordinates of the vector data will be aligned to the pixel
+ coordinates.
+
+ Notes
+ -----
+ - The data extent is the smallest rectangle that contains all the images and geometries.
+ - MultiscaleSpatialImage objects will be converted to SpatialImage objects.
+ - This helper function will be deprecated when https://github.com/scverse/spatialdata/issues/308 is closed,
+ as this function will be easily recovered by `transform_to_coordinate_system()`
+ """
+ # TODO: change "all the elements have idenity" with "all the elements have the same transformation"
+ from spatialdata._core.data_extent import get_extent
+ from spatialdata._core.operations.rasterize import _compute_target_dimensions, rasterize
+ from spatialdata._core.spatialdata import SpatialData
+ from spatialdata.transformations.operations import get_transformation, set_transformation
+ from spatialdata.transformations.transformations import BaseTransformation, Scale, Sequence, Translation
+
+ sdata = sdata.filter_by_coordinate_system(coordinate_system=coordinate_system)
+ # calling transform_to_coordinate_system will likely decrease the resolution, let's use rasterize() instead
+ sdata_vector = SpatialData(shapes=dict(sdata.shapes), points=dict(sdata.points))
+ sdata_raster = SpatialData(images=dict(sdata.images), labels=dict(sdata.labels))
+ sdata_vector_transformed = sdata_vector.transform_to_coordinate_system(coordinate_system)
+
+ de = get_extent(sdata, coordinate_system=coordinate_system)
+ de_axes = tuple(de.keys())
+ translation_to_origin = Translation([-de[ax][0] for ax in de_axes], axes=de_axes)
+
+ sizes = [de[ax][1] - de[ax][0] for ax in de_axes]
+ target_width, target_height, target_depth = _compute_target_dimensions(
+ spatial_axes=de_axes,
+ min_coordinate=[0 for _ in de_axes],
+ max_coordinate=sizes,
+ target_unit_to_pixels=target_unit_to_pixels,
+ target_width=target_width,
+ target_height=target_height,
+ target_depth=target_depth,
+ )
+ scale_to_target_d = {"x": target_width / sizes[de_axes.index("x")], "y": target_height / sizes[de_axes.index("y")]}
+ if target_depth is not None:
+ scale_to_target_d["z"] = target_depth / sizes[de_axes.index("z")]
+ scale_to_target = Scale([scale_to_target_d[ax] for ax in de_axes], axes=de_axes)
+
+ for el in sdata_vector_transformed._gen_elements_values():
+ t = get_transformation(el, to_coordinate_system=coordinate_system)
+ assert isinstance(t, BaseTransformation)
+ sequence = Sequence([t, translation_to_origin, scale_to_target])
+ set_transformation(el, transformation=sequence, to_coordinate_system=coordinate_system)
+ sdata_vector_transformed_inplace = sdata_vector_transformed.transform_to_coordinate_system(
+ coordinate_system, maintain_positioning=True
+ )
+
+ sdata_to_return_elements = {
+ **sdata_vector_transformed_inplace.shapes,
+ **sdata_vector_transformed_inplace.points,
+ }
+
+ for _, element_name, element in sdata_raster._gen_elements():
+ if isinstance(element, (MultiscaleSpatialImage, SpatialImage)):
+ rasterized = rasterize(
+ element,
+ axes=de_axes,
+ min_coordinate=[de[ax][0] for ax in de_axes],
+ max_coordinate=[de[ax][1] for ax in de_axes],
+ target_coordinate_system=coordinate_system,
+ target_unit_to_pixels=None,
+ target_width=target_width,
+ target_height=None,
+ target_depth=None,
+ )
+ sdata_to_return_elements[element_name] = rasterized
+ else:
+ sdata_to_return_elements[element_name] = element
+ if sdata.table is not None:
+ sdata_to_return_elements["table"] = sdata.table
+ return SpatialData.from_elements_dict(sdata_to_return_elements)

src/spatialdata/_core/operations/aggregate.py

@@ -170,8 +170,8 @@ def aggregate(
  target_coordinate_system, # type: ignore[assignment]
  )
  if not (by_transform == values_transform and isinstance(values_transform, Identity)):
- by_ = transform(by_, by_transform)
- values_ = transform(values_, values_transform)
+ by_ = transform(by_, to_coordinate_system=target_coordinate_system)
+ values_ = transform(values_, to_coordinate_system=target_coordinate_system)
 
  # dispatch
  adata = None

src/spatialdata/_core/operations/rasterize.py

@@ -448,7 +448,8 @@ def _(
  if schema in (Labels2DModel, Labels3DModel):
  kwargs = {"prefilter": False, "order": 0}
  elif schema in (Image2DModel, Image3DModel):
- kwargs = {}
+ kwargs = {"order": 0}
+ # kwargs = {}
  else:
  raise ValueError(f"Unsupported schema {schema}")
 

docs/api.md

@@ -32,6 +32,7 @@ Operations on `SpatialData` objects.
  transform
  rasterize
  aggregate
+ filter_by_coordinate_system
 ```
 
 ### Operations Utilities

docs/design_doc.md

@@ -564,6 +564,7 @@ with coordinate systems:
  with axes: c, y, x
  with elements: /images/point8, /labels/point8
 """
+
 sdata0 = sdata.query.coordinate_system("point23", filter_rows=False)
 sdata1 = sdata.query.bounding_box((0, 20, 0, 300))
 sdata1 = sdata.query.polygon("/polygons/annotations")

src/spatialdata/_core/data_extent.py

@@ -20,9 +20,6 @@
 from spatialdata.models._utils import SpatialElement
 from spatialdata.models.models import PointsModel
 from spatialdata.transformations.operations import get_transformation
-from spatialdata.transformations.transformations import (
- BaseTransformation,
-)
 
 BoundingBoxDescription = dict[str, tuple[float, float]]
 
@@ -289,9 +286,7 @@ def _(e: GeoDataFrame, coordinate_system: str = "global", exact: bool = True) ->
  coordinate_system=coordinate_system,
  extent=extent,
  )
- t = get_transformation(e, to_coordinate_system=coordinate_system)
- assert isinstance(t, BaseTransformation)
- transformed = transform(e, t)
+ transformed = transform(e, to_coordinate_system=coordinate_system)
  return _get_extent_of_shapes(transformed)
 
 
@@ -305,9 +300,7 @@ def _(e: DaskDataFrame, coordinate_system: str = "global", exact: bool = True) -
  coordinate_system=coordinate_system,
  extent=extent,
  )
- t = get_transformation(e, to_coordinate_system=coordinate_system)
- assert isinstance(t, BaseTransformation)
- transformed = transform(e, t)
+ transformed = transform(e, to_coordinate_system=coordinate_system)
  return _get_extent_of_points(transformed)
 
 
@@ -353,8 +346,6 @@ def _compute_extent_in_coordinate_system(
  -------
  The bounding box description in the specified coordinate system.
  """
- transformation = get_transformation(element, to_coordinate_system=coordinate_system)
- assert isinstance(transformation, BaseTransformation)
  from spatialdata._core.query._utils import get_bounding_box_corners
 
  axes = get_axes_names(element)
@@ -368,8 +359,10 @@ def _compute_extent_in_coordinate_system(
  max_coordinate=max_coordinates,
  )
  df = pd.DataFrame(corners.data, columns=corners.axis.data.tolist())
- points = PointsModel.parse(df, coordinates={k: k for k in axes})
- transformed_corners = pd.DataFrame(transform(points, transformation).compute())
+ d = get_transformation(element, get_all=True)
+ assert isinstance(d, dict)
+ points = PointsModel.parse(df, coordinates={k: k for k in axes}, transformations=d)
+ transformed_corners = pd.DataFrame(transform(points, to_coordinate_system=coordinate_system).compute())
  # Make sure min and max values are in the same order as axes
  extent = {}
  for ax in axes:

src/spatialdata/_core/operations/_utils.py

@@ -0,0 +1,117 @@
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+from multiscale_spatial_image import MultiscaleSpatialImage
+from spatial_image import SpatialImage
+
+if TYPE_CHECKING:
+ from spatialdata._core.spatialdata import SpatialData
+
+
+def transform_to_data_extent(
+ sdata: SpatialData,
+ coordinate_system: str,
+ target_unit_to_pixels: float | None = None,
+ target_width: float | None = None,
+ target_height: float | None = None,
+ target_depth: float | None = None,
+) -> SpatialData:
+ """
+ Transform the spatial data to match the data extent, keeping the positioning, and making all transformations equal.
+
+ Parameters
+ ----------
+ sdata
+ The spatial data to transform.
+ coordinate_system
+ The coordinate system to use to compute the extent and to transform the data to.
+ target_unit_to_pixels
+ The required number of pixels per unit (units in the target coordinate system) of the data that will be
+ produced.
+ target_width
+ The width of the data extent, in pixels, for the data that will be produced.
+ target_height
+ The height of the data extent, in pixels, for the data that will be produced.
+ target_depth
+ The depth of the data extent, in pixels, for the data that will be produced.
+
+ Returns
+ -------
+ SpatialData
+ The transformed spatial data with downscaled and padded images and adjusted vector coordinates; all the
+ transformations will set to Identity and the coordinates of the vector data will be aligned to the pixel
+ coordinates.
+
+ Notes
+ -----
+ - The data extent is the smallest rectangle that contains all the images and geometries.
+ - MultiscaleSpatialImage objects will be converted to SpatialImage objects.
+ - This helper function will be deprecated when https://github.com/scverse/spatialdata/issues/308 is closed,
+ as this function will be easily recovered by `transform_to_coordinate_system()`
+ """
+ # TODO: change "all the elements have idenity" with "all the elements have the same transformation"
+ from spatialdata._core.data_extent import get_extent
+ from spatialdata._core.operations.rasterize import _compute_target_dimensions, rasterize
+ from spatialdata._core.spatialdata import SpatialData
+ from spatialdata.transformations.operations import get_transformation, set_transformation
+ from spatialdata.transformations.transformations import BaseTransformation, Scale, Sequence, Translation
+
+ sdata = sdata.filter_by_coordinate_system(coordinate_system=coordinate_system)
+ # calling transform_to_coordinate_system will likely decrease the resolution, let's use rasterize() instead
+ sdata_vector = SpatialData(shapes=dict(sdata.shapes), points=dict(sdata.points))
+ sdata_raster = SpatialData(images=dict(sdata.images), labels=dict(sdata.labels))
+ sdata_vector_transformed = sdata_vector.transform_to_coordinate_system(coordinate_system)
+
+ de = get_extent(sdata, coordinate_system=coordinate_system)
+ de_axes = tuple(de.keys())
+ translation_to_origin = Translation([-de[ax][0] for ax in de_axes], axes=de_axes)
+
+ sizes = [de[ax][1] - de[ax][0] for ax in de_axes]
+ target_width, target_height, target_depth = _compute_target_dimensions(
+ spatial_axes=de_axes,
+ min_coordinate=[0 for _ in de_axes],
+ max_coordinate=sizes,
+ target_unit_to_pixels=target_unit_to_pixels,
+ target_width=target_width,
+ target_height=target_height,
+ target_depth=target_depth,
+ )
+ scale_to_target_d = {"x": target_width / sizes[de_axes.index("x")], "y": target_height / sizes[de_axes.index("y")]}
+ if target_depth is not None:
+ scale_to_target_d["z"] = target_depth / sizes[de_axes.index("z")]
+ scale_to_target = Scale([scale_to_target_d[ax] for ax in de_axes], axes=de_axes)
+
+ for el in sdata_vector_transformed._gen_elements_values():
+ t = get_transformation(el, to_coordinate_system=coordinate_system)
+ assert isinstance(t, BaseTransformation)
+ sequence = Sequence([t, translation_to_origin, scale_to_target])
+ set_transformation(el, transformation=sequence, to_coordinate_system=coordinate_system)
+ sdata_vector_transformed_inplace = sdata_vector_transformed.transform_to_coordinate_system(
+ coordinate_system, maintain_positioning=True
+ )
+
+ sdata_to_return_elements = {
+ **sdata_vector_transformed_inplace.shapes,
+ **sdata_vector_transformed_inplace.points,
+ }
+
+ for _, element_name, element in sdata_raster._gen_elements():
+ if isinstance(element, (MultiscaleSpatialImage, SpatialImage)):
+ rasterized = rasterize(
+ element,
+ axes=de_axes,
+ min_coordinate=[de[ax][0] for ax in de_axes],
+ max_coordinate=[de[ax][1] for ax in de_axes],
+ target_coordinate_system=coordinate_system,
+ target_unit_to_pixels=None,
+ target_width=target_width,
+ target_height=None,
+ target_depth=None,
+ )
+ sdata_to_return_elements[element_name] = rasterized
+ else:
+ sdata_to_return_elements[element_name] = element
+ if sdata.table is not None:
+ sdata_to_return_elements["table"] = sdata.table
+ return SpatialData.from_elements_dict(sdata_to_return_elements)

src/spatialdata/_core/operations/aggregate.py

@@ -170,8 +170,8 @@ def aggregate(
  target_coordinate_system, # type: ignore[assignment]
  )
  if not (by_transform == values_transform and isinstance(values_transform, Identity)):
- by_ = transform(by_, by_transform)
- values_ = transform(values_, values_transform)
+ by_ = transform(by_, to_coordinate_system=target_coordinate_system)
+ values_ = transform(values_, to_coordinate_system=target_coordinate_system)
 
  # dispatch
  adata = None

src/spatialdata/_core/operations/rasterize.py

@@ -448,7 +448,8 @@ def _(
  if schema in (Labels2DModel, Labels3DModel):
  kwargs = {"prefilter": False, "order": 0}
  elif schema in (Image2DModel, Image3DModel):
- kwargs = {}
+ kwargs = {"order": 0}
+ # kwargs = {}
  else:
  raise ValueError(f"Unsupported schema {schema}")
 

docs/api.md

@@ -32,6 +32,7 @@ Operations on `SpatialData` objects.
  transform
  rasterize
  aggregate
+ filter_by_coordinate_system
 ```
 
 ### Operations Utilities

docs/design_doc.md

@@ -564,6 +564,7 @@ with coordinate systems:
  with axes: c, y, x
  with elements: /images/point8, /labels/point8
 """
+
 sdata0 = sdata.query.coordinate_system("point23", filter_rows=False)
 sdata1 = sdata.query.bounding_box((0, 20, 0, 300))
 sdata1 = sdata.query.polygon("/polygons/annotations")

src/spatialdata/_core/data_extent.py

@@ -20,9 +20,6 @@
 from spatialdata.models._utils import SpatialElement
 from spatialdata.models.models import PointsModel
 from spatialdata.transformations.operations import get_transformation
-from spatialdata.transformations.transformations import (
- BaseTransformation,
-)
 
 BoundingBoxDescription = dict[str, tuple[float, float]]
 
@@ -289,9 +286,7 @@ def _(e: GeoDataFrame, coordinate_system: str = "global", exact: bool = True) ->
  coordinate_system=coordinate_system,
  extent=extent,
  )
- t = get_transformation(e, to_coordinate_system=coordinate_system)
- assert isinstance(t, BaseTransformation)
- transformed = transform(e, t)
+ transformed = transform(e, to_coordinate_system=coordinate_system)
  return _get_extent_of_shapes(transformed)
 
 
@@ -305,9 +300,7 @@ def _(e: DaskDataFrame, coordinate_system: str = "global", exact: bool = True) -
  coordinate_system=coordinate_system,
  extent=extent,
  )
- t = get_transformation(e, to_coordinate_system=coordinate_system)
- assert isinstance(t, BaseTransformation)
- transformed = transform(e, t)
+ transformed = transform(e, to_coordinate_system=coordinate_system)
  return _get_extent_of_points(transformed)
 
 
@@ -353,8 +346,6 @@ def _compute_extent_in_coordinate_system(
  -------
  The bounding box description in the specified coordinate system.
  """
- transformation = get_transformation(element, to_coordinate_system=coordinate_system)
- assert isinstance(transformation, BaseTransformation)
  from spatialdata._core.query._utils import get_bounding_box_corners
 
  axes = get_axes_names(element)
@@ -368,8 +359,10 @@ def _compute_extent_in_coordinate_system(
  max_coordinate=max_coordinates,
  )
  df = pd.DataFrame(corners.data, columns=corners.axis.data.tolist())
- points = PointsModel.parse(df, coordinates={k: k for k in axes})
- transformed_corners = pd.DataFrame(transform(points, transformation).compute())
+ d = get_transformation(element, get_all=True)
+ assert isinstance(d, dict)
+ points = PointsModel.parse(df, coordinates={k: k for k in axes}, transformations=d)
+ transformed_corners = pd.DataFrame(transform(points, to_coordinate_system=coordinate_system).compute())
  # Make sure min and max values are in the same order as axes
  extent = {}
  for ax in axes:

src/spatialdata/_core/operations/_utils.py

@@ -0,0 +1,117 @@
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+from multiscale_spatial_image import MultiscaleSpatialImage
+from spatial_image import SpatialImage
+
+if TYPE_CHECKING:
+ from spatialdata._core.spatialdata import SpatialData
+
+
+def transform_to_data_extent(
+ sdata: SpatialData,
+ coordinate_system: str,
+ target_unit_to_pixels: float | None = None,
+ target_width: float | None = None,
+ target_height: float | None = None,
+ target_depth: float | None = None,
+) -> SpatialData:
+ """
+ Transform the spatial data to match the data extent, keeping the positioning, and making all transformations equal.
+
+ Parameters
+ ----------
+ sdata
+ The spatial data to transform.
+ coordinate_system
+ The coordinate system to use to compute the extent and to transform the data to.
+ target_unit_to_pixels
+ The required number of pixels per unit (units in the target coordinate system) of the data that will be
+ produced.
+ target_width
+ The width of the data extent, in pixels, for the data that will be produced.
+ target_height
+ The height of the data extent, in pixels, for the data that will be produced.
+ target_depth
+ The depth of the data extent, in pixels, for the data that will be produced.
+
+ Returns
+ -------
+ SpatialData
+ The transformed spatial data with downscaled and padded images and adjusted vector coordinates; all the
+ transformations will set to Identity and the coordinates of the vector data will be aligned to the pixel
+ coordinates.
+
+ Notes
+ -----
+ - The data extent is the smallest rectangle that contains all the images and geometries.
+ - MultiscaleSpatialImage objects will be converted to SpatialImage objects.
+ - This helper function will be deprecated when https://github.com/scverse/spatialdata/issues/308 is closed,
+ as this function will be easily recovered by `transform_to_coordinate_system()`
+ """
+ # TODO: change "all the elements have idenity" with "all the elements have the same transformation"
+ from spatialdata._core.data_extent import get_extent
+ from spatialdata._core.operations.rasterize import _compute_target_dimensions, rasterize
+ from spatialdata._core.spatialdata import SpatialData
+ from spatialdata.transformations.operations import get_transformation, set_transformation
+ from spatialdata.transformations.transformations import BaseTransformation, Scale, Sequence, Translation
+
+ sdata = sdata.filter_by_coordinate_system(coordinate_system=coordinate_system)
+ # calling transform_to_coordinate_system will likely decrease the resolution, let's use rasterize() instead
+ sdata_vector = SpatialData(shapes=dict(sdata.shapes), points=dict(sdata.points))
+ sdata_raster = SpatialData(images=dict(sdata.images), labels=dict(sdata.labels))
+ sdata_vector_transformed = sdata_vector.transform_to_coordinate_system(coordinate_system)
+
+ de = get_extent(sdata, coordinate_system=coordinate_system)
+ de_axes = tuple(de.keys())
+ translation_to_origin = Translation([-de[ax][0] for ax in de_axes], axes=de_axes)
+
+ sizes = [de[ax][1] - de[ax][0] for ax in de_axes]
+ target_width, target_height, target_depth = _compute_target_dimensions(
+ spatial_axes=de_axes,
+ min_coordinate=[0 for _ in de_axes],
+ max_coordinate=sizes,
+ target_unit_to_pixels=target_unit_to_pixels,
+ target_width=target_width,
+ target_height=target_height,
+ target_depth=target_depth,
+ )
+ scale_to_target_d = {"x": target_width / sizes[de_axes.index("x")], "y": target_height / sizes[de_axes.index("y")]}
+ if target_depth is not None:
+ scale_to_target_d["z"] = target_depth / sizes[de_axes.index("z")]
+ scale_to_target = Scale([scale_to_target_d[ax] for ax in de_axes], axes=de_axes)
+
+ for el in sdata_vector_transformed._gen_elements_values():
+ t = get_transformation(el, to_coordinate_system=coordinate_system)
+ assert isinstance(t, BaseTransformation)
+ sequence = Sequence([t, translation_to_origin, scale_to_target])
+ set_transformation(el, transformation=sequence, to_coordinate_system=coordinate_system)
+ sdata_vector_transformed_inplace = sdata_vector_transformed.transform_to_coordinate_system(
+ coordinate_system, maintain_positioning=True
+ )
+
+ sdata_to_return_elements = {
+ **sdata_vector_transformed_inplace.shapes,
+ **sdata_vector_transformed_inplace.points,
+ }
+
+ for _, element_name, element in sdata_raster._gen_elements():
+ if isinstance(element, (MultiscaleSpatialImage, SpatialImage)):
+ rasterized = rasterize(
+ element,
+ axes=de_axes,
+ min_coordinate=[de[ax][0] for ax in de_axes],
+ max_coordinate=[de[ax][1] for ax in de_axes],
+ target_coordinate_system=coordinate_system,
+ target_unit_to_pixels=None,
+ target_width=target_width,
+ target_height=None,
+ target_depth=None,
+ )
+ sdata_to_return_elements[element_name] = rasterized
+ else:
+ sdata_to_return_elements[element_name] = element
+ if sdata.table is not None:
+ sdata_to_return_elements["table"] = sdata.table
+ return SpatialData.from_elements_dict(sdata_to_return_elements)

src/spatialdata/_core/operations/aggregate.py

@@ -170,8 +170,8 @@ def aggregate(
  target_coordinate_system, # type: ignore[assignment]
  )
  if not (by_transform == values_transform and isinstance(values_transform, Identity)):
- by_ = transform(by_, by_transform)
- values_ = transform(values_, values_transform)
+ by_ = transform(by_, to_coordinate_system=target_coordinate_system)
+ values_ = transform(values_, to_coordinate_system=target_coordinate_system)
 
  # dispatch
  adata = None

src/spatialdata/_core/operations/rasterize.py

@@ -448,7 +448,8 @@ def _(
  if schema in (Labels2DModel, Labels3DModel):
  kwargs = {"prefilter": False, "order": 0}
  elif schema in (Image2DModel, Image3DModel):
- kwargs = {}
+ kwargs = {"order": 0}
+ # kwargs = {}
  else:
  raise ValueError(f"Unsupported schema {schema}")
 

docs/api.md

@@ -32,6 +32,7 @@ Operations on `SpatialData` objects.
  transform
  rasterize
  aggregate
+ filter_by_coordinate_system
 ```
 
 ### Operations Utilities

docs/design_doc.md

@@ -564,6 +564,7 @@ with coordinate systems:
  with axes: c, y, x
  with elements: /images/point8, /labels/point8
 """
+
 sdata0 = sdata.query.coordinate_system("point23", filter_rows=False)
 sdata1 = sdata.query.bounding_box((0, 20, 0, 300))
 sdata1 = sdata.query.polygon("/polygons/annotations")

src/spatialdata/_core/data_extent.py

@@ -20,9 +20,6 @@
 from spatialdata.models._utils import SpatialElement
 from spatialdata.models.models import PointsModel
 from spatialdata.transformations.operations import get_transformation
-from spatialdata.transformations.transformations import (
- BaseTransformation,
-)
 
 BoundingBoxDescription = dict[str, tuple[float, float]]
 
@@ -289,9 +286,7 @@ def _(e: GeoDataFrame, coordinate_system: str = "global", exact: bool = True) ->
  coordinate_system=coordinate_system,
  extent=extent,
  )
- t = get_transformation(e, to_coordinate_system=coordinate_system)
- assert isinstance(t, BaseTransformation)
- transformed = transform(e, t)
+ transformed = transform(e, to_coordinate_system=coordinate_system)
  return _get_extent_of_shapes(transformed)
 
 
@@ -305,9 +300,7 @@ def _(e: DaskDataFrame, coordinate_system: str = "global", exact: bool = True) -
  coordinate_system=coordinate_system,
  extent=extent,
  )
- t = get_transformation(e, to_coordinate_system=coordinate_system)
- assert isinstance(t, BaseTransformation)
- transformed = transform(e, t)
+ transformed = transform(e, to_coordinate_system=coordinate_system)
  return _get_extent_of_points(transformed)
 
 
@@ -353,8 +346,6 @@ def _compute_extent_in_coordinate_system(
  -------
  The bounding box description in the specified coordinate system.
  """
- transformation = get_transformation(element, to_coordinate_system=coordinate_system)
- assert isinstance(transformation, BaseTransformation)
  from spatialdata._core.query._utils import get_bounding_box_corners
 
  axes = get_axes_names(element)
@@ -368,8 +359,10 @@ def _compute_extent_in_coordinate_system(
  max_coordinate=max_coordinates,
  )
  df = pd.DataFrame(corners.data, columns=corners.axis.data.tolist())
- points = PointsModel.parse(df, coordinates={k: k for k in axes})
- transformed_corners = pd.DataFrame(transform(points, transformation).compute())
+ d = get_transformation(element, get_all=True)
+ assert isinstance(d, dict)
+ points = PointsModel.parse(df, coordinates={k: k for k in axes}, transformations=d)
+ transformed_corners = pd.DataFrame(transform(points, to_coordinate_system=coordinate_system).compute())
  # Make sure min and max values are in the same order as axes
  extent = {}
  for ax in axes:

src/spatialdata/_core/operations/_utils.py

@@ -0,0 +1,117 @@
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+from multiscale_spatial_image import MultiscaleSpatialImage
+from spatial_image import SpatialImage
+
+if TYPE_CHECKING:
+ from spatialdata._core.spatialdata import SpatialData
+
+
+def transform_to_data_extent(
+ sdata: SpatialData,
+ coordinate_system: str,
+ target_unit_to_pixels: float | None = None,
+ target_width: float | None = None,
+ target_height: float | None = None,
+ target_depth: float | None = None,
+) -> SpatialData:
+ """
+ Transform the spatial data to match the data extent, keeping the positioning, and making all transformations equal.
+
+ Parameters
+ ----------
+ sdata
+ The spatial data to transform.
+ coordinate_system
+ The coordinate system to use to compute the extent and to transform the data to.
+ target_unit_to_pixels
+ The required number of pixels per unit (units in the target coordinate system) of the data that will be
+ produced.
+ target_width
+ The width of the data extent, in pixels, for the data that will be produced.
+ target_height
+ The height of the data extent, in pixels, for the data that will be produced.
+ target_depth
+ The depth of the data extent, in pixels, for the data that will be produced.
+
+ Returns
+ -------
+ SpatialData
+ The transformed spatial data with downscaled and padded images and adjusted vector coordinates; all the
+ transformations will set to Identity and the coordinates of the vector data will be aligned to the pixel
+ coordinates.
+
+ Notes
+ -----
+ - The data extent is the smallest rectangle that contains all the images and geometries.
+ - MultiscaleSpatialImage objects will be converted to SpatialImage objects.
+ - This helper function will be deprecated when https://github.com/scverse/spatialdata/issues/308 is closed,
+ as this function will be easily recovered by `transform_to_coordinate_system()`
+ """
+ # TODO: change "all the elements have idenity" with "all the elements have the same transformation"
+ from spatialdata._core.data_extent import get_extent
+ from spatialdata._core.operations.rasterize import _compute_target_dimensions, rasterize
+ from spatialdata._core.spatialdata import SpatialData
+ from spatialdata.transformations.operations import get_transformation, set_transformation
+ from spatialdata.transformations.transformations import BaseTransformation, Scale, Sequence, Translation
+
+ sdata = sdata.filter_by_coordinate_system(coordinate_system=coordinate_system)
+ # calling transform_to_coordinate_system will likely decrease the resolution, let's use rasterize() instead
+ sdata_vector = SpatialData(shapes=dict(sdata.shapes), points=dict(sdata.points))
+ sdata_raster = SpatialData(images=dict(sdata.images), labels=dict(sdata.labels))
+ sdata_vector_transformed = sdata_vector.transform_to_coordinate_system(coordinate_system)
+
+ de = get_extent(sdata, coordinate_system=coordinate_system)
+ de_axes = tuple(de.keys())
+ translation_to_origin = Translation([-de[ax][0] for ax in de_axes], axes=de_axes)
+
+ sizes = [de[ax][1] - de[ax][0] for ax in de_axes]
+ target_width, target_height, target_depth = _compute_target_dimensions(
+ spatial_axes=de_axes,
+ min_coordinate=[0 for _ in de_axes],
+ max_coordinate=sizes,
+ target_unit_to_pixels=target_unit_to_pixels,
+ target_width=target_width,
+ target_height=target_height,
+ target_depth=target_depth,
+ )
+ scale_to_target_d = {"x": target_width / sizes[de_axes.index("x")], "y": target_height / sizes[de_axes.index("y")]}
+ if target_depth is not None:
+ scale_to_target_d["z"] = target_depth / sizes[de_axes.index("z")]
+ scale_to_target = Scale([scale_to_target_d[ax] for ax in de_axes], axes=de_axes)
+
+ for el in sdata_vector_transformed._gen_elements_values():
+ t = get_transformation(el, to_coordinate_system=coordinate_system)
+ assert isinstance(t, BaseTransformation)
+ sequence = Sequence([t, translation_to_origin, scale_to_target])
+ set_transformation(el, transformation=sequence, to_coordinate_system=coordinate_system)
+ sdata_vector_transformed_inplace = sdata_vector_transformed.transform_to_coordinate_system(
+ coordinate_system, maintain_positioning=True
+ )
+
+ sdata_to_return_elements = {
+ **sdata_vector_transformed_inplace.shapes,
+ **sdata_vector_transformed_inplace.points,
+ }
+
+ for _, element_name, element in sdata_raster._gen_elements():
+ if isinstance(element, (MultiscaleSpatialImage, SpatialImage)):
+ rasterized = rasterize(
+ element,
+ axes=de_axes,
+ min_coordinate=[de[ax][0] for ax in de_axes],
+ max_coordinate=[de[ax][1] for ax in de_axes],
+ target_coordinate_system=coordinate_system,
+ target_unit_to_pixels=None,
+ target_width=target_width,
+ target_height=None,
+ target_depth=None,
+ )
+ sdata_to_return_elements[element_name] = rasterized
+ else:
+ sdata_to_return_elements[element_name] = element
+ if sdata.table is not None:
+ sdata_to_return_elements["table"] = sdata.table
+ return SpatialData.from_elements_dict(sdata_to_return_elements)

src/spatialdata/_core/operations/aggregate.py

@@ -170,8 +170,8 @@ def aggregate(
  target_coordinate_system, # type: ignore[assignment]
  )
  if not (by_transform == values_transform and isinstance(values_transform, Identity)):
- by_ = transform(by_, by_transform)
- values_ = transform(values_, values_transform)
+ by_ = transform(by_, to_coordinate_system=target_coordinate_system)
+ values_ = transform(values_, to_coordinate_system=target_coordinate_system)
 
  # dispatch
  adata = None

src/spatialdata/_core/operations/rasterize.py

@@ -448,7 +448,8 @@ def _(
  if schema in (Labels2DModel, Labels3DModel):
  kwargs = {"prefilter": False, "order": 0}
  elif schema in (Image2DModel, Image3DModel):
- kwargs = {}
+ kwargs = {"order": 0}
+ # kwargs = {}
  else:
  raise ValueError(f"Unsupported schema {schema}")
 

docs/api.md

@@ -32,6 +32,7 @@ Operations on `SpatialData` objects.
  transform
  rasterize
  aggregate
+ filter_by_coordinate_system
 ```
 
 ### Operations Utilities

docs/design_doc.md

@@ -564,6 +564,7 @@ with coordinate systems:
  with axes: c, y, x
  with elements: /images/point8, /labels/point8
 """
+
 sdata0 = sdata.query.coordinate_system("point23", filter_rows=False)
 sdata1 = sdata.query.bounding_box((0, 20, 0, 300))
 sdata1 = sdata.query.polygon("/polygons/annotations")

src/spatialdata/_core/data_extent.py

@@ -20,9 +20,6 @@
 from spatialdata.models._utils import SpatialElement
 from spatialdata.models.models import PointsModel
 from spatialdata.transformations.operations import get_transformation
-from spatialdata.transformations.transformations import (
- BaseTransformation,
-)
 
 BoundingBoxDescription = dict[str, tuple[float, float]]
 
@@ -289,9 +286,7 @@ def _(e: GeoDataFrame, coordinate_system: str = "global", exact: bool = True) ->
  coordinate_system=coordinate_system,
  extent=extent,
  )
- t = get_transformation(e, to_coordinate_system=coordinate_system)
- assert isinstance(t, BaseTransformation)
- transformed = transform(e, t)
+ transformed = transform(e, to_coordinate_system=coordinate_system)
  return _get_extent_of_shapes(transformed)
 
 
@@ -305,9 +300,7 @@ def _(e: DaskDataFrame, coordinate_system: str = "global", exact: bool = True) -
  coordinate_system=coordinate_system,
  extent=extent,
  )
- t = get_transformation(e, to_coordinate_system=coordinate_system)
- assert isinstance(t, BaseTransformation)
- transformed = transform(e, t)
+ transformed = transform(e, to_coordinate_system=coordinate_system)
  return _get_extent_of_points(transformed)
 
 
@@ -353,8 +346,6 @@ def _compute_extent_in_coordinate_system(
  -------
  The bounding box description in the specified coordinate system.
  """
- transformation = get_transformation(element, to_coordinate_system=coordinate_system)
- assert isinstance(transformation, BaseTransformation)
  from spatialdata._core.query._utils import get_bounding_box_corners
 
  axes = get_axes_names(element)
@@ -368,8 +359,10 @@ def _compute_extent_in_coordinate_system(
  max_coordinate=max_coordinates,
  )
  df = pd.DataFrame(corners.data, columns=corners.axis.data.tolist())
- points = PointsModel.parse(df, coordinates={k: k for k in axes})
- transformed_corners = pd.DataFrame(transform(points, transformation).compute())
+ d = get_transformation(element, get_all=True)
+ assert isinstance(d, dict)
+ points = PointsModel.parse(df, coordinates={k: k for k in axes}, transformations=d)
+ transformed_corners = pd.DataFrame(transform(points, to_coordinate_system=coordinate_system).compute())
  # Make sure min and max values are in the same order as axes
  extent = {}
  for ax in axes:

src/spatialdata/_core/operations/_utils.py

@@ -0,0 +1,117 @@
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+from multiscale_spatial_image import MultiscaleSpatialImage
+from spatial_image import SpatialImage
+
+if TYPE_CHECKING:
+ from spatialdata._core.spatialdata import SpatialData
+
+
+def transform_to_data_extent(
+ sdata: SpatialData,
+ coordinate_system: str,
+ target_unit_to_pixels: float | None = None,
+ target_width: float | None = None,
+ target_height: float | None = None,
+ target_depth: float | None = None,
+) -> SpatialData:
+ """
+ Transform the spatial data to match the data extent, keeping the positioning, and making all transformations equal.
+
+ Parameters
+ ----------
+ sdata
+ The spatial data to transform.
+ coordinate_system
+ The coordinate system to use to compute the extent and to transform the data to.
+ target_unit_to_pixels
+ The required number of pixels per unit (units in the target coordinate system) of the data that will be
+ produced.
+ target_width
+ The width of the data extent, in pixels, for the data that will be produced.
+ target_height
+ The height of the data extent, in pixels, for the data that will be produced.
+ target_depth
+ The depth of the data extent, in pixels, for the data that will be produced.
+
+ Returns
+ -------
+ SpatialData
+ The transformed spatial data with downscaled and padded images and adjusted vector coordinates; all the
+ transformations will set to Identity and the coordinates of the vector data will be aligned to the pixel
+ coordinates.
+
+ Notes
+ -----
+ - The data extent is the smallest rectangle that contains all the images and geometries.
+ - MultiscaleSpatialImage objects will be converted to SpatialImage objects.
+ - This helper function will be deprecated when https://github.com/scverse/spatialdata/issues/308 is closed,
+ as this function will be easily recovered by `transform_to_coordinate_system()`
+ """
+ # TODO: change "all the elements have idenity" with "all the elements have the same transformation"
+ from spatialdata._core.data_extent import get_extent
+ from spatialdata._core.operations.rasterize import _compute_target_dimensions, rasterize
+ from spatialdata._core.spatialdata import SpatialData
+ from spatialdata.transformations.operations import get_transformation, set_transformation
+ from spatialdata.transformations.transformations import BaseTransformation, Scale, Sequence, Translation
+
+ sdata = sdata.filter_by_coordinate_system(coordinate_system=coordinate_system)
+ # calling transform_to_coordinate_system will likely decrease the resolution, let's use rasterize() instead
+ sdata_vector = SpatialData(shapes=dict(sdata.shapes), points=dict(sdata.points))
+ sdata_raster = SpatialData(images=dict(sdata.images), labels=dict(sdata.labels))
+ sdata_vector_transformed = sdata_vector.transform_to_coordinate_system(coordinate_system)
+
+ de = get_extent(sdata, coordinate_system=coordinate_system)
+ de_axes = tuple(de.keys())
+ translation_to_origin = Translation([-de[ax][0] for ax in de_axes], axes=de_axes)
+
+ sizes = [de[ax][1] - de[ax][0] for ax in de_axes]
+ target_width, target_height, target_depth = _compute_target_dimensions(
+ spatial_axes=de_axes,
+ min_coordinate=[0 for _ in de_axes],
+ max_coordinate=sizes,
+ target_unit_to_pixels=target_unit_to_pixels,
+ target_width=target_width,
+ target_height=target_height,
+ target_depth=target_depth,
+ )
+ scale_to_target_d = {"x": target_width / sizes[de_axes.index("x")], "y": target_height / sizes[de_axes.index("y")]}
+ if target_depth is not None:
+ scale_to_target_d["z"] = target_depth / sizes[de_axes.index("z")]
+ scale_to_target = Scale([scale_to_target_d[ax] for ax in de_axes], axes=de_axes)
+
+ for el in sdata_vector_transformed._gen_elements_values():
+ t = get_transformation(el, to_coordinate_system=coordinate_system)
+ assert isinstance(t, BaseTransformation)
+ sequence = Sequence([t, translation_to_origin, scale_to_target])
+ set_transformation(el, transformation=sequence, to_coordinate_system=coordinate_system)
+ sdata_vector_transformed_inplace = sdata_vector_transformed.transform_to_coordinate_system(
+ coordinate_system, maintain_positioning=True
+ )
+
+ sdata_to_return_elements = {
+ **sdata_vector_transformed_inplace.shapes,
+ **sdata_vector_transformed_inplace.points,
+ }
+
+ for _, element_name, element in sdata_raster._gen_elements():
+ if isinstance(element, (MultiscaleSpatialImage, SpatialImage)):
+ rasterized = rasterize(
+ element,
+ axes=de_axes,
+ min_coordinate=[de[ax][0] for ax in de_axes],
+ max_coordinate=[de[ax][1] for ax in de_axes],
+ target_coordinate_system=coordinate_system,
+ target_unit_to_pixels=None,
+ target_width=target_width,
+ target_height=None,
+ target_depth=None,
+ )
+ sdata_to_return_elements[element_name] = rasterized
+ else:
+ sdata_to_return_elements[element_name] = element
+ if sdata.table is not None:
+ sdata_to_return_elements["table"] = sdata.table
+ return SpatialData.from_elements_dict(sdata_to_return_elements)

src/spatialdata/_core/operations/aggregate.py

@@ -170,8 +170,8 @@ def aggregate(
  target_coordinate_system, # type: ignore[assignment]
  )
  if not (by_transform == values_transform and isinstance(values_transform, Identity)):
- by_ = transform(by_, by_transform)
- values_ = transform(values_, values_transform)
+ by_ = transform(by_, to_coordinate_system=target_coordinate_system)
+ values_ = transform(values_, to_coordinate_system=target_coordinate_system)
 
  # dispatch
  adata = None

src/spatialdata/_core/operations/rasterize.py

@@ -448,7 +448,8 @@ def _(
  if schema in (Labels2DModel, Labels3DModel):
  kwargs = {"prefilter": False, "order": 0}
  elif schema in (Image2DModel, Image3DModel):
- kwargs = {}
+ kwargs = {"order": 0}
+ # kwargs = {}
  else:
  raise ValueError(f"Unsupported schema {schema}")