xenium: label hex cell_ids without label_id; drop degenerate boundary polygons

src/spatialdata_io/readers/xenium.py

@@ -442,8 +442,12 @@ def xenium(
         # Map the table to the cell_labels element when a label_index mapping exists, so the instance
         # key resolves against the raster labels instead of the circles. If any table cell lacks a
         # boundary (so has no label), warn and leave the table on circles rather than crash.
+        label_index = None
         if cells_labels and cell_im is not None:
             label_index = _cell_id_to_label_index(cell_im)
+        elif cells_labels and cells_zarr_ctx is None:
+            label_index = _hex_cell_id_to_label_index(table.obs[str(XeniumKeys.CELL_ID)])
+        if label_index is not None:
             table_ids = table.obs[str(XeniumKeys.CELL_ID)].to_numpy()
             if pd.Index(table_ids).isin(label_index.index).all():
                 table.obs["cell_labels"] = label_index.loc[table_ids].to_numpy()
@@ -539,6 +543,28 @@ def _decode_cell_id_column(cell_id_column: pd.Series) -> pd.Series:
     return cell_id_column
 
 
+def _hex_cell_id_to_label_index(cell_ids: ArrayLike) -> pd.Series | None:
+    """cell_id (str) -> label_index (int) for hex cell_ids that have no ``label_id`` mapping.
+
+    XOA 1.3-1.x boundary files carry hex cell_ids but no ``label_id``; the mask values live only in
+    cells.zarr.zip, which stripped (e.g. GEO) deposits omit. The shifted-hex prefix decodes to the
+    cell's integer id, which is unique within a single dataset suffix, so it labels the reconstructed
+    raster and the table consistently. It need not equal the original XOA mask value. Returns None
+    when the ids are not decodable hex ids, span several dataset suffixes, decode to 0 (background)
+    or are not unique.
+    """
+    ids = np.asarray(cell_ids).astype(str)
+    if not len(ids):
+        return None
+    try:
+        prefix, suffix = prefix_suffix_uint32_from_cell_id_str(ids)
+    except (KeyError, ValueError):
+        return None
+    if len(np.unique(suffix)) != 1 or (prefix == 0).any() or len(np.unique(prefix)) != len(prefix):
+        return None
+    return pd.Series(prefix.astype(np.int64), index=ids)
+
+
 def _cell_id_to_label_index(indices_mapping: pd.DataFrame) -> pd.Series:
     """cell_id (str) -> label_index (int) lookup from an indices mapping."""
     return indices_mapping.set_index("cell_id")["label_index"]
@@ -637,6 +663,24 @@ def _get_polygons(
 
     group_ends = np.concatenate([group_starts[1:], [n]])
 
+    # A polygon ring needs >= 4 coordinates; XOA occasionally emits 2-3 vertex boundaries, which make
+    # shapely raise. Drop them, except when indices_mapping must align 1:1 with the polygons.
+    sizes = group_ends - group_starts
+    keep = sizes >= 4
+    if not keep.all() and (has_label_id or indices_mapping is None):
+        warnings.warn(
+            f"Dropping {int((~keep).sum())} degenerate polygons with fewer than 4 vertices from {file}.",
+            UserWarning,
+            stacklevel=3,
+        )
+        rows = np.repeat(keep, sizes)
+        coords = coords[rows]
+        id_col = id_col.filter(pa.array(rows))
+        if has_label_id:
+            id_arr = id_arr[rows]
+        group_ends = np.cumsum(sizes[keep])
+        group_starts = group_ends - sizes[keep]
+
     # offsets for ragged array:
     # offsets[0] (ring_offsets): describing to which rings the vertex positions belong to
     # offsets[1] (geom_offsets): describing to which polygons the rings belong to
@@ -690,7 +734,7 @@ def _get_polygons(
         #   fashion to the cells, so we could modify get_indices_mapping NOT to return None
         #   in that case, which would imply that we don't fall into this slow fallback. This
         #   could be work for a future version.
-        unique_ids = id_col.filter(np.concatenate([[True], change_mask])).to_pylist()
+        unique_ids = id_col.take(pa.array(group_starts)).to_pylist()
         index = _decode_cell_id_column(pd.Series(unique_ids))
         geo_df = GeoDataFrame({"geometry": geoms}, index=index.values)
         geo_df.index.name = "cell_id"
@@ -729,10 +773,11 @@ def _get_labels_from_boundaries(
         label = _cell_id_to_label_index(indices_mapping).loc[shapes.index].to_numpy()
     elif pd.api.types.is_integer_dtype(shapes.index):
         label = shapes.index.to_numpy()
+    elif (hex_label := _hex_cell_id_to_label_index(shapes.index)) is not None:
+        label = hex_label.to_numpy()
     else:
-        # Non-integer (hex) cell_ids with no label_index mapping: we have no integer to fill the
-        # raster with. For CSV bundles this cannot happen (hex ids always ship a ``label_id`` column,
-        # so ``indices_mapping`` is set), but guard it rather than crash cryptically on ``astype``.
+        # Non-integer cell_ids that are neither mapped by ``label_id`` nor decodable as 10x hex ids:
+        # we have no integer to fill the raster with, so fail clearly rather than on ``astype``.
         raise ValueError(
             "Cannot build raster labels: boundary polygons have non-integer cell_ids but no "
             "label_index mapping (missing 'label_id' column). Please report this at "

tests/test_xenium.py

@@ -340,6 +340,42 @@ def test_xenium_csv_only_hex_labels() -> None:
     assert len(_label_ids(sdata["nucleus_labels"])) >= 12  # >= cells (multinucleate)
 
 
+def test_xenium_degenerate_boundary_dropped(tmp_path: Path) -> None:
+    # XOA occasionally emits boundaries with < 4 coordinates, which shapely cannot turn into a ring
+    import shutil
+
+    for f in XENIUM_CSV_ONLY.iterdir():
+        shutil.copy(f, tmp_path / f.name)
+    b = pd.read_csv(XENIUM_CSV_ONLY / "cell_boundaries.csv.gz")
+    first = b["cell_id"].iloc[0]
+    b = pd.concat([b[b["cell_id"] == first].head(3), b[b["cell_id"] != first]])
+    b.to_csv(tmp_path / "cell_boundaries.csv.gz", index=False)
+    with pytest.warns(UserWarning, match="Dropping 1 degenerate polygons"):
+        sdata = xenium(tmp_path, morphology_mip=False, morphology_focus=False, aligned_images=False)
+    assert len(sdata["cell_boundaries"]) == 29
+    assert first not in sdata["cell_boundaries"].index
+    assert sdata["table"].n_obs == 30
+
+
+def test_xenium_csv_only_hex_without_label_id(tmp_path: Path) -> None:
+    # XOA 1.3-1.x deposits stripped of cells.zarr.zip: hex cell_ids but no label_id column. The labels
+    # fall back to the decoded hex prefix, consistently in the raster and in the table.
+    import shutil
+
+    for f in XENIUM_CSV_HEX.iterdir():
+        if f.name.endswith("_boundaries.csv.gz"):
+            pd.read_csv(f).drop(columns="label_id").to_csv(tmp_path / f.name, index=False)
+        else:
+            shutil.copy(f, tmp_path / f.name)
+    sdata = xenium(tmp_path, morphology_mip=False, morphology_focus=False, aligned_images=False, transcripts=False)
+    table = sdata["table"]
+    prefix, _ = prefix_suffix_uint32_from_cell_id_str(table.obs["cell_id"].to_numpy())
+    np.testing.assert_array_equal(table.obs["cell_labels"].to_numpy(), prefix)
+    assert _label_ids(sdata["cell_labels"]) == set(prefix.tolist())
+    assert table.uns["spatialdata_attrs"]["instance_key"] == "cell_labels"
+    assert _label_ids(sdata["nucleus_labels"]) <= set(prefix.tolist())
+
+
 @skip_if_below_python_version()
 @pytest.mark.parametrize(
     (