feat: Add lazy table loading via anndata.experimental.read_lazy

src/spatialdata/_core/query/relational_query.py

@@ -246,6 +246,18 @@ def _region_as_str_if_list_of_len_one(region: list[str]) -> str | list[str]:
     return region if len(region) > 1 else region[0]
 
 
+def _obs_as_dataframe(table: AnnData) -> pd.DataFrame:
+    """Return ``table.obs`` as a pandas DataFrame.
+
+    Tables read with ``anndata.experimental.read_lazy`` expose ``obs`` as an xarray
+    ``Dataset2D``, which does not implement the pandas methods (``reset_index``,
+    ``groupby``) the join helpers rely on. Materialize it in that case; ``obs`` is the
+    small axis of the table, so this does not pull in ``X``, which stays lazy.
+    """
+    obs = table.obs
+    return obs if isinstance(obs, pd.DataFrame) else obs.to_memory()
+
+
 def _right_exclusive_join_spatialelement_table(
     element_dict: dict[str, dict[str, Any]],
     table: AnnData,
@@ -256,7 +268,7 @@ def _right_exclusive_join_spatialelement_table(
     if isinstance(regions, str):
         regions = [regions]
     # reset_index so group_df.index gives integer positions — safe with duplicate obs names
-    obs = table.obs.reset_index()
+    obs = _obs_as_dataframe(table).reset_index()
     groups_df = obs.groupby(by=region_column_name, observed=False)
     keep = np.zeros(len(table), dtype=bool)
     has_match = False
@@ -301,7 +313,7 @@ def _right_join_spatialelement_table(
     regions, region_column_name, instance_key = get_table_keys(table)
     if isinstance(regions, str):
         regions = [regions]
-    groups_df = table.obs.groupby(by=region_column_name, observed=False)
+    groups_df = _obs_as_dataframe(table).groupby(by=region_column_name, observed=False)
     for element_type, name_element in element_dict.items():
         for name, element in name_element.items():
             if name in regions:
@@ -343,7 +355,7 @@ def _inner_join_spatialelement_table(
     regions, region_column_name, instance_key = get_table_keys(table)
     if isinstance(regions, str):
         regions = [regions]
-    obs = table.obs.reset_index()
+    obs = _obs_as_dataframe(table).reset_index()
     groups_df = obs.groupby(by=region_column_name, observed=False)
     joined_indices = None
     for element_type, name_element in element_dict.items():
@@ -404,7 +416,7 @@ def _left_exclusive_join_spatialelement_table(
     regions, region_column_name, instance_key = get_table_keys(table)
     if isinstance(regions, str):
         regions = [regions]
-    groups_df = table.obs.groupby(by=region_column_name, observed=False)
+    groups_df = _obs_as_dataframe(table).groupby(by=region_column_name, observed=False)
     for element_type, name_element in element_dict.items():
         for name, element in name_element.items():
             if name in regions:
@@ -442,7 +454,7 @@ def _left_join_spatialelement_table(
     regions, region_column_name, instance_key = get_table_keys(table)
     if isinstance(regions, str):
         regions = [regions]
-    obs = table.obs.reset_index()
+    obs = _obs_as_dataframe(table).reset_index()
     groups_df = obs.groupby(by=region_column_name, observed=False)
     joined_indices = None
     for element_type, name_element in element_dict.items():
@@ -1076,7 +1088,12 @@ def get_values(
         if origin == "obs":
             df = obs[value_key_values].copy()
         if origin == "var":
-            matched_table.obs = pd.DataFrame(obs)
+            # When the table came from anndata.experimental.read_lazy, obs is a Dataset2D, not a
+            # DataFrame, and pd.DataFrame(obs) returns a malformed frame. Materialize via to_memory().
+            if isinstance(obs, pd.DataFrame):
+                matched_table.obs = pd.DataFrame(obs)
+            else:
+                matched_table.obs = obs.to_memory()
             if table_layer is None:
                 x = matched_table[:, value_key_values].X
             else:
@@ -1085,6 +1102,11 @@ def get_values(
                 x = matched_table[:, value_key_values].layers[table_layer]
             import scipy
 
+            if isinstance(x, da.Array):
+                # A lazy table backs X with dask, and pd.DataFrame() cannot consume that.
+                # get_values returns in-memory values by contract, and this is a selection of
+                # the requested columns only, so the materialization is bounded by those.
+                x = x.compute()
             if isinstance(x, scipy.sparse.csr_matrix | scipy.sparse.csc_matrix | scipy.sparse.coo_matrix):
                 x = x.todense()
             df = pd.DataFrame(x, columns=value_key_values)

src/spatialdata/_core/spatialdata.py

@@ -1868,6 +1868,7 @@ def read(
         file_path: str | Path | UPath | zarr.Group,
         selection: tuple[str] | None = None,
         reconsolidate_metadata: bool = False,
+        lazy: bool = False,
     ) -> SpatialData:
         """
         Read a SpatialData object from a Zarr storage (on-disk or remote).
@@ -1880,6 +1881,28 @@ def read(
             The elements to read (images, labels, points, shapes, table). If None, all elements are read.
         reconsolidate_metadata
             If the consolidated metadata store got corrupted this can lead to errors when trying to read the data.
+        lazy
+            If True, read tables lazily using anndata.experimental.read_lazy.
+            This keeps large tables out of memory until needed. Requires anndata >= 0.12.
+            Note: Images, labels, and points are always read lazily (using Dask).
+            This parameter only affects tables, which are normally loaded into memory.
+
+            When the stored ``X`` is sparse, the lazy table's ``X`` is a Dask array whose
+            blocks are ``scipy.sparse`` matrices, and Dask's array reductions
+            (``X.sum()``, ``X.mean()``, ``X.max()``, ``X.std()``, ...) are **not
+            supported**. They raise a ``TypeError`` or ``IndexError`` while the graph is
+            being built -- before ``.compute()`` is ever reached -- because Dask derives
+            the result metadata by calling the corresponding NumPy reduction on a
+            ``scipy.sparse`` block, and ``scipy.sparse`` does not accept the
+            ``keepdims``/``ndmin`` arguments NumPy passes down. This is a
+            Dask/``scipy.sparse`` interoperability limitation, not something this reader
+            introduces. Slicing and ``.compute()`` work normally, so reduce a materialized
+            block instead::
+
+                table.X[:1000].compute().sum()      # works
+                table.X.sum()                        # raises
+
+            or use ``dask.array.map_blocks`` with a function that handles sparse blocks.
 
         Returns
         -------
@@ -1892,7 +1915,7 @@ def read(
 
             _write_consolidated_metadata(file_path)
 
-        return read_zarr(file_path, selection=selection)
+        return read_zarr(file_path, selection=selection, lazy=lazy)
 
     @property
     def images(self) -> Images:

src/spatialdata/_io/io_table.py

@@ -24,8 +24,34 @@
 from spatialdata.models import TableModel, get_table_keys
 
 
-def _read_table(store: str | Path) -> AnnData:
-    table = read_anndata_zarr(str(store))
+def _read_table(store: str | Path, lazy: bool = False) -> AnnData:
+    """
+    Read a table from a zarr store.
+
+    Parameters
+    ----------
+    store
+        Path to the zarr store containing the table.
+    lazy
+        If True, read the table lazily using ``anndata.experimental.read_lazy``.
+        This keeps large matrices (X, layers) as dask arrays backed by zarr,
+        so they are only loaded into memory on demand. Requires anndata >= 0.12.
+
+    Returns
+    -------
+    The AnnData table, either lazily loaded or in-memory.
+
+    Raises
+    ------
+    ImportError
+        If ``lazy=True`` but anndata >= 0.12 is not installed.
+    """
+    if lazy:
+        from anndata.experimental import read_lazy
+
+        table = read_lazy(str(store))
+    else:
+        table = read_anndata_zarr(str(store))
 
     f = zarr.open(Path(store), mode="r")  # Path avoids zarr v3 URL-parsing special chars (e.g. #) in names
     version = _parse_version(f, expect_attrs_key=False)

src/spatialdata/_io/io_zarr.py

@@ -3,6 +3,7 @@
 import os
 import warnings
 from collections.abc import Callable
+from functools import partial
 from json import JSONDecodeError
 from pathlib import Path
 from typing import Any, Literal
@@ -126,6 +127,7 @@ def read_zarr(
     store: str | Path | UPath | zarr.Group,
     selection: None | tuple[str] = None,
     on_bad_files: Literal[BadFileHandleMethod.ERROR, BadFileHandleMethod.WARN] = BadFileHandleMethod.ERROR,
+    lazy: bool = False,
 ) -> SpatialData:
     """
     Read a SpatialData dataset from a zarr store (on-disk or remote).
@@ -149,6 +151,28 @@ def read_zarr(
           object is returned containing only elements that could be read. Failures can only be
           determined from the warnings.
 
+    lazy
+        If True, read tables lazily using anndata.experimental.read_lazy.
+        This keeps large tables out of memory until needed. Requires anndata >= 0.12.
+        Note: Images, labels, and points are always read lazily (using Dask).
+        This parameter only affects tables, which are normally loaded into memory.
+
+        When the stored ``X`` is sparse, the lazy table's ``X`` is a Dask array whose
+        blocks are ``scipy.sparse`` matrices, and Dask's array reductions
+        (``X.sum()``, ``X.mean()``, ``X.max()``, ``X.std()``, ...) are **not supported**.
+        They raise a ``TypeError`` or ``IndexError`` while the graph is being built --
+        before ``.compute()`` is ever reached -- because Dask derives the result metadata
+        by calling the corresponding NumPy reduction on a ``scipy.sparse`` block, and
+        ``scipy.sparse`` does not accept the ``keepdims``/``ndmin`` arguments NumPy passes
+        down. This is a Dask/``scipy.sparse`` interoperability limitation, not something
+        this reader introduces. Slicing and ``.compute()`` work normally, so reduce a
+        materialized block instead::
+
+            table.X[:1000].compute().sum()      # works
+            table.X.sum()                        # raises
+
+        or use ``dask.array.map_blocks`` with a function that handles sparse blocks.
+
     Returns
     -------
     A SpatialData object.
@@ -195,7 +219,7 @@ def read_zarr(
         "labels": (_read_multiscale, "labels", labels),
         "points": (_read_points, "points", points),
         "shapes": (_read_shapes, "shapes", shapes),
-        "tables": (_read_table, "tables", tables),
+        "tables": (partial(_read_table, lazy=lazy), "tables", tables),
     }
     for group_name, (
         read_func,

src/spatialdata/_utils.py

@@ -225,11 +225,17 @@ def _inplace_fix_subset_categorical_obs(subset_adata: AnnData, original_adata: A
     """
     if not hasattr(subset_adata, "obs") or not hasattr(original_adata, "obs"):
         return
-    obs = pd.DataFrame(subset_adata.obs)
+    # Tables read via anndata.experimental.read_lazy have a Dataset2D obs instead of a DataFrame;
+    # pd.DataFrame() would silently malform it, so materialize with to_memory() in that case.
+    obs = pd.DataFrame(subset_adata.obs) if isinstance(subset_adata.obs, pd.DataFrame) else subset_adata.obs.to_memory()
+    original_obs = (
+        original_adata.obs if isinstance(original_adata.obs, pd.DataFrame) else original_adata.obs.to_memory()
+    )
+
     for column in obs.columns:
         is_categorical = isinstance(obs[column].dtype, pd.CategoricalDtype)
         if is_categorical:
-            c = obs[column].cat.set_categories(original_adata.obs[column].cat.categories)
+            c = obs[column].cat.set_categories(original_obs[column].cat.categories)
             obs[column] = c
     subset_adata.obs = obs
 

src/spatialdata/models/models.py

@@ -55,6 +55,25 @@
 ATTRS_KEY = "spatialdata_attrs"
 
 
+def _is_lazy_anndata(adata: AnnData) -> bool:
+    """Check if an AnnData object is lazily loaded.
+
+    Lazy AnnData objects (from anndata.experimental.read_lazy) have obs/var
+    stored as xarray Dataset2D instead of pandas DataFrame.
+
+    Parameters
+    ----------
+    adata
+        The AnnData object to check.
+
+    Returns
+    -------
+    True if the AnnData is lazily loaded, False otherwise.
+    """
+    # Check if obs is not a pandas DataFrame (lazy AnnData uses xarray Dataset2D)
+    return not isinstance(adata.obs, pd.DataFrame)
+
+
 def _parse_transformations(element: SpatialElement, transformations: MappingToCoordinateSystem_t | None = None) -> None:
     _validate_mapping_to_coordinate_system_type(transformations)
     transformations_in_element = _get_transformations(element)
@@ -1085,6 +1104,13 @@ def _validate_table_annotation_metadata(cls, data: AnnData) -> None:
             raise ValueError(f"`{attr[cls.REGION_KEY_KEY]}` not found in `adata.obs`. Please create the column.")
         if attr[cls.INSTANCE_KEY] not in data.obs:
             raise ValueError(f"`{attr[cls.INSTANCE_KEY]}` not found in `adata.obs`. Please create the column.")
+
+        # Skip detailed dtype/value validation for lazy-loaded AnnData
+        # These checks would trigger data loading, defeating the purpose of lazy loading
+        # Validation will occur when data is actually computed/accessed
+        if _is_lazy_anndata(data):
+            return
+
         instance_col = data.obs[attr[cls.INSTANCE_KEY]]
         dtype = instance_col.dtype
 
@@ -1154,14 +1180,19 @@ def validate(
         if ATTRS_KEY not in data.uns:
             return data
 
+        # Check if this is a lazy-loaded AnnData (from anndata.experimental.read_lazy)
+        # Lazy AnnData has xarray-based obs/var, which requires different validation
+        is_lazy = _is_lazy_anndata(data)
+
         _, region_key, instance_key = get_table_keys(data)
         if region_key is not None:
             if region_key not in data.obs:
                 raise ValueError(
                     f"Region key `{region_key}` not in `adata.obs`. Please create the column and parse "
                     f"using TableModel.parse(adata)."
                 )
-            if not isinstance(data.obs[region_key].dtype, CategoricalDtype):
+            # Skip dtype validation for lazy tables (would require loading data)
+            if not is_lazy and not isinstance(data.obs[region_key].dtype, CategoricalDtype):
                 raise ValueError(
                     f"`table.obs[{region_key}]` must be of type `categorical`, not `{type(data.obs[region_key])}`."
                 )
@@ -1171,7 +1202,8 @@ def validate(
                     f"Instance key `{instance_key}` not in `adata.obs`. Please create the column and parse"
                     f" using TableModel.parse(adata)."
                 )
-            if data.obs[instance_key].isnull().values.any():
+            # Skip null check for lazy tables (would require loading data)
+            if not is_lazy and data.obs[instance_key].isnull().values.any():
                 raise ValueError("`table.obs[instance_key]` must not contain null values, but it does.")
 
         cls._validate_table_annotation_metadata(data)