rasterize converts labels to images

[aeisenbarth]

Describe the bug

When applying transform_to_data_extent on a SpatialData with labels, they are merged into images.

To Reproduce
Steps to reproduce the behavior:

from spatialdata.datasets import blobs
sdata = blobs()
sdata_rasterized = transform_to_data_extent(sdata, coordinate_system="global", target_unit_to_pixels=1)
assert len(sdata_rasterized.labels) == len(sdata.labels)

In [61]: sdata_rasterized
Out[61]: 
SpatialData object
├── Images
│     ├── 'blobs_image': DataArray[cyx] (3, 512, 512)
│     ├── 'blobs_labels': DataArray[cyx] (1, 512, 512)
│     ├── 'blobs_multiscale_image': DataArray[cyx] (3, 512, 512)
│     └── 'blobs_multiscale_labels': DataArray[cyx] (1, 512, 512)
├── Points
│     └── 'blobs_points': DataFrame with shape: (<Delayed>, 4) (2D points)
├── Shapes
│     ├── 'blobs_circles': GeoDataFrame shape: (5, 2) (2D shapes)
│     ├── 'blobs_multipolygons': GeoDataFrame shape: (2, 1) (2D shapes)
│     └── 'blobs_polygons': GeoDataFrame shape: (5, 1) (2D shapes)
└── Tables
      └── 'table': AnnData (26, 3)
with coordinate systems:
    ▸ 'global', with elements:
        blobs_image (Images), blobs_labels (Images), blobs_multiscale_image (Images), blobs_multiscale_labels (Images), blobs_points (Points), blobs_circles (Shapes), blobs_multipolygons (Shapes), blobs_polygons (Shapes)

Expected behavior
All element types should be preserved. Labels should stay in labels.

Additional context
This is caused by the following line:

spatialdata/src/spatialdata/_core/operations/rasterize.py

Lines 330 to 331 in 6694c91

	if model in (Labels2DModel, Labels3DModel) and not return_regions_as_labels:
	model = Image2DModel if model == Labels2DModel else Image3DModel

The behavior is undocumented and not intuitive. As it was intentionally introduced, there seems to exist a use case, but in other use cases it is a potential source of bugs. Some control over it is already provided through the flag return_regions_as_labels=False.

• We should document it.
• We should rethink whether it is safer to change the default of return_regions_as_labels, especially if it is for internal use.
• We must either make sure all functions calling rasterize pass the correct value for return_regions_as_labels that is appropriate in each case, or further expose the flag to their callers.
  A possible fix would be to make transform_to_data_extent call rasterize(…, return_regions_as_labels=True)

  spatialdata/src/spatialdata/_core/operations/_utils.py

  Lines 118 to 128 in 6694c91

  	rasterized = rasterize(
  	element,
  	axes=data_extent_axes,
  	min_coordinate=[data_extent[ax][0] for ax in data_extent_axes],
  	max_coordinate=[data_extent[ax][1] for ax in data_extent_axes],
  	target_coordinate_system=coordinate_system,
  	target_unit_to_pixels=None,
  	target_width=target_width,
  	target_height=None,
  	target_depth=None,
  	)