exceeding rasterize's maximum label index

[LouisK92]

I am using sd.rasterize on cell polygons of >100k cells and run into the error.

---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[14], line 1
----> 1 sdata["cell_labels"] = sd.rasterize(
...
 9 )

File ~/miniconda3/envs/g3/lib/python3.11/site-packages/spatialdata/_core/operations/rasterize.py:351, in rasterize(data, axes, min_coordinate, max_coordinate, target_coordinate_system, target_unit_to_pixels, target_width, target_height, target_depth, sdata, value_key, table_name, return_regions_as_labels, agg_func, return_single_channel)
 349 return rasterized
 350 if model in (PointsModel, ShapesModel):
--> 351 return rasterize_shapes_points(
... 368 )
 369 raise ValueError(f"Unsupported model {model}.")

File ~/miniconda3/envs/g3/lib/python3.11/site-packages/spatialdata/_core/operations/rasterize.py:723, in rasterize_shapes_points(data, axes, min_coordinate, max_coordinate, target_coordinate_system, target_unit_to_pixels, target_width, target_height, target_depth, element_name, sdata, value_key, table_name, return_regions_as_labels, agg_func, return_single_channel)
 721 max_uint16 = np.iinfo(np.uint16).max
 722 if max_label > max_uint16:
--> 723 raise ValueError(f"Maximum label index is {max_label}. Values higher than {max_uint16} are not supported.")
 724 agg = agg.astype(np.uint16)
 725 return Labels2DModel.parse(agg, transformations=transformations)

ValueError: Maximum label index is 117340. Values higher than 65535 are not supported.

I was wondering if the threshold needs to be that low or could be increased?

To reproduce the error:

import numpy as np
import spatialdata as sd

def get_n_coords(n):
 coords = np.random.rand(n, 2)
 return coords
 
sdata1 = sd.SpatialData(points={"points": PointsModel.parse(get_n_coords(65534))})
sdata2 = sd.SpatialData(points={"points": PointsModel.parse(get_n_coords(65536))})

# sdata1 would work
img_extent = sd.get_extent(sdata2['points'])
sd.rasterize(
 sdata2["points"],
 ["x", "y"],
 min_coordinate=[int(img_extent["x"][0]), int(img_extent["y"][0])],
 max_coordinate=[int(img_extent["x"][1]), int(img_extent["y"][1])],
 target_coordinate_system="global",
 target_unit_to_pixels=1,
 return_regions_as_labels=True,
)

[LouisK92]

A quick workaround (assuming that cell labels don't overlap) would be the following. Note that I used circles instead of points for the example, which is closer to what I actually wanted. In case of points the .iloc fails since we'd deal with a dask dataframe in that case.

import numpy as np
import spatialdata as sd
from spatialdata.models import ShapesModel

def get_n_coords(n):
 coords = np.random.rand(n, 2)
 return coords
 
sdata = sd.SpatialData(shapes={"shapes": ShapesModel.parse(get_n_coords(65536), geometry=0, radius=np.ones(65536))})

img_extent = sd.get_extent(sdata['shapes']) # NOTE: Here I actually use sd.get_extent(sdata['image'])

N = 65535 # number of shapes that are rasterised per step
n_cells = len(sdata['shapes'])
n_iter = n_cells // N + bool(n_cells % N)

rasterize_args = {
 "min_coordinate": [int(img_extent["x"][0]), int(img_extent["y"][0])],
 "max_coordinate": [int(img_extent["x"][1]), int(img_extent["y"][1])],
 "target_coordinate_system": "global",
 "target_unit_to_pixels": 1,
 "return_regions_as_labels": True,
}

for i in range(n_iter):
 labels_image_ = sd.rasterize(sdata['shapes'].iloc[i*N:min((i+1)*N, n_cells)], ["x", "y"], **rasterize_args)
 if i == 0:
 labels_image = labels_image_
 else:
 labels_image.values[labels_image_.values > 0] = labels_image_.values[labels_image_.values > 0]

sdata['labels'] = labels_image

[LucaMarconato]

Related discussion: #187