diff --git a/src/spatialdata/_core/centroids.py b/src/spatialdata/_core/centroids.py
index cde583e1..060b6feb 100644
--- a/src/spatialdata/_core/centroids.py
+++ b/src/spatialdata/_core/centroids.py
@@ -1,65 +1,150 @@
 from __future__ import annotations
 
 from functools import singledispatch
+from typing import Literal
 
 import numpy as np
 import pandas as pd
 import xarray as xr
+from anndata import AnnData
 from dask.dataframe import DataFrame as DaskDataFrame
 from geopandas import GeoDataFrame
 from shapely import MultiPolygon, Point, Polygon
 from xarray import DataArray, DataTree
 
 from spatialdata._core.operations.transform import transform
-from spatialdata.models import get_axes_names
+from spatialdata._core.spatialdata import SpatialData
+from spatialdata.models import get_axes_names, get_table_keys
 from spatialdata.models._utils import SpatialElement
-from spatialdata.models.models import Labels2DModel, Labels3DModel, PointsModel, get_model
+from spatialdata.models.models import Labels2DModel, Labels3DModel, PointsModel, ShapesModel, TableModel, get_model
 from spatialdata.transformations.operations import get_transformation
 from spatialdata.transformations.transformations import BaseTransformation
 
 BoundingBoxDescription = dict[str, tuple[float, float]]
 
+PersistAs = Literal["Points", "adata"]
 
-def _validate_coordinate_system(e: SpatialElement, coordinate_system: str) -> None:
+
+def _validate_coordinate_system(e: SpatialElement, coordinate_system: str) -> dict[str, BaseTransformation]:
+    """Return the element's ``{coordinate_system: transformation}`` map, raising if ``coordinate_system`` is absent."""
     d = get_transformation(e, get_all=True)
     assert isinstance(d, dict)
-    assert coordinate_system in d, (
-        f"No transformation to coordinate system {coordinate_system} is available for the given element.\n"
-        f"Available coordinate systems: {list(d.keys())}"
-    )
+    if coordinate_system not in d:
+        raise ValueError(
+            f"No transformation to coordinate system {coordinate_system!r} is available for the given element. "
+            f"Available coordinate systems: {list(d.keys())}"
+        )
+    return d
+
+
+def _validate_persist_args(persist_as: str, coordinate_system: str | None, *, allow_adata: bool) -> None:
+    if persist_as not in ("Points", "adata"):
+        raise ValueError(f"`persist_as` must be 'Points' or 'adata', got {persist_as!r}.")
+    if persist_as == "adata" and not allow_adata:
+        raise ValueError(
+            "persist_as='adata' writes centroids into the element's annotating table, which needs the "
+            "`SpatialData` object: call `get_centroids(sdata, element_name, ..., persist_as='adata')`. "
+            "To get the centroids as a standalone element instead, use persist_as='Points'."
+        )
+    # ``coordinate_system=None`` means "intrinsic coordinates, do not transform". An intrinsic Points
+    # element is ill-defined (Points always carry a coordinate system), so intrinsic coords are only
+    # meaningful when writing into a table (persist_as='adata').
+    if coordinate_system is None and persist_as != "adata":
+        raise ValueError("`coordinate_system=None` (intrinsic coordinates) is only supported with persist_as='adata'.")
+
+
+def _is_identity(transformation: BaseTransformation, e: SpatialElement) -> bool:
+    """Return True if ``transformation`` is a mathematical no-op (its affine matrix equals the identity).
+
+    Compares the affine matrix rather than checking ``isinstance(..., Identity)``, so a numerically
+    identity ``Scale([1, 1])`` / ``Translation([0, 0])`` / ``Affine(np.eye(...))`` / ``Sequence([Identity()])``
+    also counts as intrinsic.
+    """
+    axes = get_axes_names(e)
+    matrix = transformation.to_affine_matrix(input_axes=axes, output_axes=axes)
+    return bool(np.allclose(matrix, np.eye(len(axes) + 1)))
+
+
+def _validate_return_area(
+    e: SpatialElement, coordinate_system: str, transformation: BaseTransformation, return_area: bool
+) -> None:
+    """Validate ``return_area`` for the ``persist_as="Points"`` path.
+
+    Points have no area (raises). For labels/shapes the area is computed in intrinsic (untransformed)
+    units, but the Points centroids *are* transformed into ``coordinate_system``; returning both under a
+    non-identity transform would leave area and centroids in different coordinate systems, so intrinsic
+    coordinates (a system whose transformation is the identity) are required.
+    """
+    if not return_area:
+        return
+    if get_model(e) is PointsModel:
+        raise ValueError("`return_area` is not supported for points elements (points have no area).")
+    if not _is_identity(transformation, e):
+        raise ValueError(
+            f"`return_area=True` requires intrinsic coordinates: pass `coordinate_system=None` or a "
+            f"coordinate system whose transformation is the identity. The area is not transformed, so "
+            f"pairing it with the non-identity transform to {coordinate_system!r} would put area and "
+            f"centroids in different coordinate systems."
+        )
 
 
 @singledispatch
 def get_centroids(
-    e: SpatialElement,
-    coordinate_system: str = "global",
+    e: SpatialElement | SpatialData,
+    coordinate_system: str | None = "global",
     return_background: bool = False,
-) -> DaskDataFrame:
+    return_area: bool = False,
+    persist_as: PersistAs = "Points",
+) -> DaskDataFrame | AnnData | None:
     """
-    Get the centroids of the geometries contained in a SpatialElement, as a new Points element.
+    Get the centroids of the geometries contained in a SpatialElement.
 
     Parameters
     ----------
     e
-        The SpatialElement. Only points, shapes (circles, polygons and multipolygons) and labels are supported.
+        The SpatialElement (points, shapes — circles, polygons and multipolygons — or labels), or a
+        :class:`~spatialdata.SpatialData` object. When a ``SpatialData`` is passed, the second
+        positional argument is the name of the element to measure (see the ``SpatialData`` overload).
     coordinate_system
-        The coordinate system in which the centroids are computed.
+        The coordinate system in which the centroids are computed. ``None`` returns the intrinsic
+        coordinates without applying any transformation (only supported with ``persist_as="adata"``).
     return_background
-        If True, the centroid of the background label (0) is included in the output.
+        If True, the centroid of the background label (0) is included in the output (labels only).
+    return_area
+        If True, also return the per-instance area: the pixel/voxel count for labels and the geometric
+        area for shapes (``pi * r**2`` for circles). Not supported for points (raises). With
+        ``persist_as="Points"`` the area is added as a feature column of the returned Points element.
+        The area is always computed in intrinsic (untransformed) units, so it is only allowed together
+        with intrinsic coordinates: ``return_area=True`` requires ``coordinate_system=None`` or a
+        coordinate system whose transformation is the identity (otherwise area and centroids would live
+        in different coordinate systems).
+    persist_as
+        ``"Points"`` (default) returns the centroids as a new Points element, transformed into
+        ``coordinate_system``. ``"adata"`` writes the centroids (and area) into the element's
+        annotating table and is only available through the :class:`~spatialdata.SpatialData` overload,
+        which can resolve that table.
+
+    Returns
+    -------
+    A Points element (``persist_as="Points"``). With ``persist_as="adata"`` (``SpatialData`` overload),
+    ``None`` when written in place, or the new ``AnnData`` table when ``inplace=False``.
 
     Notes
     -----
-    For :class:`~shapely.Multipolygon`s, the centroids are the average of the centroids of the polygons that constitute
-    each :class:`~shapely.Multipolygon`.
+    For :class:`~shapely.MultiPolygon`s, the centroid is the geometric (area-weighted) centroid computed by
+    ``geopandas`` (``GeoDataFrame.centroid``), i.e. sub-polygons contribute in proportion to their area. For
+    multiscale labels the centroids are computed on the full-resolution ``scale0`` level.
     """
     raise ValueError(f"The object type {type(e)} is not supported.")
 
 
-def _get_centroids_for_labels(xdata: xr.DataArray) -> pd.DataFrame:
+def _get_centroids_for_labels(xdata: xr.DataArray, return_area: bool = False) -> pd.DataFrame:
     """
     Compute centroids for all labels in a DataArray in a single O(n_voxels) pass.
 
-    Works for any number of spatial dimensions (2D and 3D labels).
+    Works for any number of spatial dimensions (2D and 3D labels). When ``return_area`` is True, an
+    ``area`` column (the per-label pixel/voxel count) is added; it is already computed for the
+    centroids, so this is free.
     """
     arr = xdata.data.compute()
     axes = list(xdata.dims)
@@ -74,69 +159,253 @@ def _get_centroids_for_labels(xdata: xr.DataArray) -> pd.DataFrame:
     coord_grids = np.meshgrid(*[xdata[ax].values for ax in axes], indexing="ij")
     data: dict[str, np.ndarray] = {}
     for ax, grid in zip(axes, coord_grids, strict=True):
-        coord_sums = np.bincount(flat_inverse, weights=grid.ravel().astype(float))
+        coord_sums = np.bincount(flat_inverse, weights=grid.ravel().astype(float, copy=False))
         data[ax] = coord_sums / counts  # counts > 0 by construction (unique guarantees this)
 
-    return pd.DataFrame(data, index=label_ids)
+    df = pd.DataFrame(data, index=label_ids)
+    if return_area:
+        df["area"] = counts.astype(float)
+    return df
 
 
-@get_centroids.register(DataArray)
-@get_centroids.register(DataTree)
-def _(
-    e: DataArray | DataTree,
-    coordinate_system: str = "global",
-    return_background: bool = False,
-) -> DaskDataFrame:
-    """Get the centroids of a Labels element (2D or 3D)."""
-    model = get_model(e)
-    if model not in [Labels2DModel, Labels3DModel]:
-        raise ValueError("Expected a `Labels` element. Found an `Image` instead.")
-    _validate_coordinate_system(e, coordinate_system)
-
-    if isinstance(e, DataTree):
-        assert len(e["scale0"]) == 1
-        e = next(iter(e["scale0"].values()))
-
-    df = _get_centroids_for_labels(e)
-    if not return_background and 0 in df.index:
-        df = df.drop(index=0)  # drop the background label
-    t = get_transformation(e, coordinate_system)
-    centroids = PointsModel.parse(df, transformations={coordinate_system: t})
-    return transform(centroids, to_coordinate_system=coordinate_system)
-
-
-@get_centroids.register(GeoDataFrame)
-def _(e: GeoDataFrame, coordinate_system: str = "global") -> DaskDataFrame:
-    """Get the centroids of a Shapes element (circles or polygons/multipolygons)."""
-    _validate_coordinate_system(e, coordinate_system)
-    t = get_transformation(e, coordinate_system)
-    assert isinstance(t, BaseTransformation)
-    # separate points from (multi-)polygons
+def _get_centroids_for_shapes(e: GeoDataFrame, return_area: bool) -> tuple[pd.DataFrame, np.ndarray | None]:
+    """Intrinsic per-shape centroids (``x, y`` columns indexed by the element's index) and optional area."""
     first_geometry = e["geometry"].iloc[0]
     if isinstance(first_geometry, Point):
         xy = e.geometry.get_coordinates().values
+        # shapely .area is 0 for circles (Point geometry); the radius column carries the size.
+        area = np.pi * np.asarray(e["radius"], dtype=float) ** 2 if return_area else None
     else:
         assert isinstance(first_geometry, Polygon | MultiPolygon), (
             f"Expected a GeoDataFrame either composed entirely of circles (Points with the `radius` column) or"
             f" Polygons/MultiPolygons. Found {type(first_geometry)} instead."
         )
         xy = e.centroid.get_coordinates().values
+        area = e.geometry.area.to_numpy() if return_area else None
     xy_df = pd.DataFrame(xy, columns=["x", "y"], index=e.index.copy())
-    points = PointsModel.parse(xy_df, transformations={coordinate_system: t})
-    return transform(points, to_coordinate_system=coordinate_system)
+    return xy_df, area
 
 
-@get_centroids.register(DaskDataFrame)
-def _(e: DaskDataFrame, coordinate_system: str = "global") -> DaskDataFrame:
-    """Get the centroids of a Points element."""
-    _validate_coordinate_system(e, coordinate_system)
-    axes = get_axes_names(e)
-    assert axes in [("x", "y"), ("x", "y", "z")]
-    coords = e[list(axes)].compute()
+def _intrinsic_centroid_frame(
+    element: SpatialElement, return_background: bool, return_area: bool
+) -> tuple[pd.DataFrame, np.ndarray | None, SpatialElement]:
+    """Per-instance intrinsic centroids (coordinate columns, indexed by instance id), optional area.
+
+    Also returns the element the centroids live on (for labels, the ``scale0`` level of a multiscale
+    raster), which carries the transformation to apply downstream.
+    """
+    model = get_model(element)
+    if model in (Labels2DModel, Labels3DModel):
+        raster = next(iter(element["scale0"].values())) if isinstance(element, DataTree) else element
+        df = _get_centroids_for_labels(raster, return_area=return_area)
+        if not return_background and 0 in df.index:
+            df = df.drop(index=0)  # drop the background label (its area, if any, goes with it)
+        area = df.pop("area").to_numpy() if return_area else None
+        return df, area, raster
+    if model is ShapesModel:
+        xy_df, area = _get_centroids_for_shapes(element, return_area)
+        return xy_df, area, element
+    if model is PointsModel:
+        if return_area:
+            raise ValueError("`return_area` is not supported for points elements (points have no area).")
+        axes = get_axes_names(element)
+        if axes not in [("x", "y"), ("x", "y", "z")]:
+            raise ValueError(f"Expected points axes to be ('x', 'y') or ('x', 'y', 'z'), got {axes}.")
+        return element[list(axes)].compute(), None, element
+    raise ValueError(
+        f"Centroids are not supported for elements modeled by {model.__name__}; "
+        f"expected a Labels, Shapes or Points element."
+    )
+
+
+def _points_from_centroids(
+    df: pd.DataFrame, area: np.ndarray | None, e: SpatialElement, coordinate_system: str
+) -> DaskDataFrame:
+    """Build a Points element from intrinsic centroids, transformed into ``coordinate_system``."""
+    out = df.assign(area=np.asarray(area, dtype=float)) if area is not None else df
     t = get_transformation(e, coordinate_system)
     assert isinstance(t, BaseTransformation)
-    centroids = PointsModel.parse(coords, transformations={coordinate_system: t})
-    return transform(centroids, to_coordinate_system=coordinate_system)
+    points = PointsModel.parse(out, transformations={coordinate_system: t})
+    return transform(points, to_coordinate_system=coordinate_system)
+
+
+@get_centroids.register(DataArray)
+@get_centroids.register(DataTree)
+@get_centroids.register(GeoDataFrame)
+@get_centroids.register(DaskDataFrame)
+def _(
+    e: SpatialElement,
+    coordinate_system: str | None = "global",
+    return_background: bool = False,
+    return_area: bool = False,
+    persist_as: PersistAs = "Points",
+) -> DaskDataFrame:
+    """Get the centroids of a Labels, Shapes or Points element."""
+    _validate_persist_args(persist_as, coordinate_system, allow_adata=False)
+    assert coordinate_system is not None  # guaranteed by _validate_persist_args (allow_adata=False)
+    transformations = _validate_coordinate_system(e, coordinate_system)
+    _validate_return_area(e, coordinate_system, transformations[coordinate_system], return_area)
+    df, area, raster = _intrinsic_centroid_frame(e, return_background, return_area)
+    return _points_from_centroids(df, area, raster, coordinate_system)
+
+
+def _resolve_annotating_table(sdata: SpatialData, element_name: str, table_name: str | None) -> str:
+    """Resolve the single table that annotates ``element_name`` (where centroids are written)."""
+    from spatialdata._core.query.relational_query import get_element_annotators
+
+    if table_name is not None:
+        if table_name not in sdata.tables:
+            raise KeyError(f"Table {table_name!r} not found in `sdata.tables`.")
+        return table_name
+    annotators = sorted(get_element_annotators(sdata, element_name))
+    if not annotators:
+        raise ValueError(
+            f"Element {element_name!r} has no annotating table to write centroids into. Use "
+            f"persist_as='Points' to get the centroids as a Points element instead, or annotate the "
+            f"element with a table first."
+        )
+    if len(annotators) > 1:
+        raise ValueError(
+            f"Element {element_name!r} is annotated by multiple tables ({', '.join(annotators)}); "
+            f"pass `table_name=` to choose one."
+        )
+    return annotators[0]
+
+
+def _write_centroids_into_table(
+    table: AnnData,
+    element_name: str,
+    centroids: pd.DataFrame,
+    area: np.ndarray | None,
+    overwrite: bool,
+) -> None:
+    """Write centroids into ``obsm["spatial"]`` and area into ``obs["area"]`` at the element's rows.
+
+    Only the table rows annotating ``element_name`` are touched (a table may annotate several
+    elements); instances annotated but absent from the element are written as NaN. Rows annotating
+    *other* elements are preserved. If the target already holds values for the element's own rows,
+    the write is refused unless ``overwrite=True`` (guards against clobbering data — e.g. a
+    pre-existing, unrelated ``obs["area"]`` column).
+    """
+    if not centroids.index.is_unique:
+        raise ValueError(f"Cannot persist centroids for {element_name!r}: its instance index has duplicate values.")
+    _, region_key, instance_key = get_table_keys(table)
+    mask = (table.obs[region_key].astype(str) == str(element_name)).to_numpy()
+    if not mask.any():
+        raise ValueError(f"The resolved table does not annotate element {element_name!r} (no matching rows).")
+
+    # Map each annotated instance to its centroid row (-1 where absent -> NaN). A *total* miss means the
+    # instance ids never align with the element index (mismatched dtype, or no shared instances).
+    keys = table.obs[instance_key].to_numpy()[mask]
+    idx = centroids.index.get_indexer(keys)
+    if (idx == -1).all():
+        raise ValueError(
+            f"No instance id annotating {element_name!r} is present in the element; check the table's "
+            f"`{instance_key}` values and dtype."
+        )
+    hit = idx != -1
+
+    def _scatter(values: np.ndarray) -> np.ndarray:
+        """Gather ``values`` (ordered like ``centroids``) onto the masked rows, NaN where absent."""
+        out = np.full((len(idx), *values.shape[1:]), np.nan)
+        out[hit] = values[idx[hit]]
+        return out
+
+    ndim = centroids.shape[1]
+    spatial = np.full((table.n_obs, ndim), np.nan)
+    existing = table.obsm.get(TableModel.SPATIAL_KEY)
+    if existing is not None:
+        existing = np.asarray(existing)
+        if existing.shape == (table.n_obs, ndim):
+            spatial = existing.astype(float, copy=True)  # preserve other regions' coordinates
+            if not overwrite and np.isfinite(spatial[mask]).any():
+                raise ValueError(
+                    f"obsm['{TableModel.SPATIAL_KEY}'] already holds coordinates for rows annotating "
+                    f"{element_name!r}; pass overwrite=True to replace them."
+                )
+        elif existing.shape[0] == table.n_obs:
+            raise ValueError(
+                f"Existing obsm['{TableModel.SPATIAL_KEY}'] {existing.shape} is incompatible with {ndim}-D "
+                f"centroids for "
+                f"{element_name!r}; refusing to overwrite other regions. Persist with persist_as='Points' instead."
+            )
+    spatial[mask] = _scatter(centroids.to_numpy(dtype=float))
+    table.obsm[TableModel.SPATIAL_KEY] = spatial
+
+    if area is not None:
+        col = np.full(table.n_obs, np.nan)
+        if TableModel.AREA_KEY in table.obs:
+            existing_area = table.obs[TableModel.AREA_KEY]
+            if not overwrite and (existing_area.notna().to_numpy() & mask).any():
+                raise ValueError(
+                    f"obs['{TableModel.AREA_KEY}'] already holds values for rows annotating {element_name!r} "
+                    f"(this may be an unrelated column); pass overwrite=True to replace them."
+                )
+            # Preserve values outside the element's rows; non-numeric entries coerce to NaN (only reachable
+            # with overwrite=True, where replacing that column is the intent).
+            col = pd.to_numeric(existing_area, errors="coerce").to_numpy(dtype=float, copy=True)
+        col[mask] = _scatter(np.asarray(area, dtype=float))
+        table.obs[TableModel.AREA_KEY] = col
+
+
+@get_centroids.register(SpatialData)
+def _get_centroids_sdata(
+    e: SpatialData,
+    element_name: str,
+    coordinate_system: str | None = "global",
+    return_background: bool = False,
+    return_area: bool = False,
+    persist_as: PersistAs = "Points",
+    table_name: str | None = None,
+    inplace: bool = True,
+    overwrite: bool = False,
+) -> DaskDataFrame | AnnData | None:
+    """Get the centroids of ``element_name``, or (``persist_as="adata"``) write them into its annotating table.
+
+    With ``persist_as="adata"`` the centroids go into ``obsm["spatial"]`` (and area into ``obs["area"]``) of the
+    resolved annotating table (``table_name=`` disambiguates). ``inplace=True`` (default) mutates that table and
+    returns ``None``; ``inplace=False`` writes into a copy of *only that table* and returns the new ``AnnData``,
+    leaving ``e`` untouched. ``overwrite`` (default ``False``) guards the write: if ``obsm["spatial"]`` or
+    ``obs["area"]`` already holds values for the element's own rows (e.g. a pre-existing, unrelated ``obs["area"]``
+    column), the write is refused unless ``overwrite=True``. ``persist_as="Points"`` behaves like calling
+    :func:`get_centroids` on the element.
+
+    ``obsm["spatial"]`` stores **intrinsic** (untransformed) coordinates: we do not persist the coordinate
+    system alongside the centroids, but each row already records which element it annotates (via the table's
+    region/instance keys), so the intrinsic system of every row is recoverable and can be transformed
+    downstream on demand. Accordingly ``persist_as="adata"`` requires intrinsic coordinates —
+    ``coordinate_system=None`` or a system whose transformation is the identity. The columns are always
+    ordered ``x, y[, z]`` (never ``y, x``), regardless of the element's own axis order.
+    """
+    _validate_persist_args(persist_as, coordinate_system, allow_adata=True)
+    element = e[element_name]
+
+    if persist_as == "Points":
+        return get_centroids(
+            element,
+            coordinate_system=coordinate_system,
+            return_background=return_background,
+            return_area=return_area,
+        )
+
+    # persist_as == "adata": resolve the annotating table and write the intrinsic centroids into it.
+    if coordinate_system is not None:
+        transformations = _validate_coordinate_system(element, coordinate_system)
+        if not _is_identity(transformations[coordinate_system], element):
+            raise ValueError(
+                f"persist_as='adata' stores intrinsic coordinates in obsm['{TableModel.SPATIAL_KEY}'] and cannot "
+                f"apply the non-identity transform to {coordinate_system!r}. Pass coordinate_system=None (or a "
+                f"coordinate system whose transformation is the identity); transform per element downstream if "
+                f"you need another coordinate system."
+            )
+    table_name = _resolve_annotating_table(e, element_name, table_name)
+    df, area, _ = _intrinsic_centroid_frame(element, return_background, return_area)
+    coord_cols = sorted(df.columns)  # canonical x, y[, z] (squidpy obsm["spatial"] order)
+    centroids = df[coord_cols]
+
+    table = e.tables[table_name] if inplace else e.tables[table_name].copy()
+    _write_centroids_into_table(table, element_name, centroids, area, overwrite)
+    return None if inplace else table
 
 
 ##
diff --git a/src/spatialdata/_core/operations/vectorize.py b/src/spatialdata/_core/operations/vectorize.py
index 41458458..d12a3f1a 100644
--- a/src/spatialdata/_core/operations/vectorize.py
+++ b/src/spatialdata/_core/operations/vectorize.py
@@ -1,7 +1,7 @@
 from __future__ import annotations
 
 from functools import singledispatch
-from typing import TYPE_CHECKING, Any
+from typing import TYPE_CHECKING, Any, cast
 
 import dask
 import numpy as np
@@ -137,7 +137,8 @@ def _get_centroids(element: SpatialElement) -> pd.DataFrame:
     if INTRINSIC_COORDINATE_SYSTEM in d:
         raise RuntimeError(f"The name {INTRINSIC_COORDINATE_SYSTEM} is reserved.")
     d[INTRINSIC_COORDINATE_SYSTEM] = Identity()
-    centroids = get_centroids(element, coordinate_system=INTRINSIC_COORDINATE_SYSTEM).compute()
+    # get_centroids returns DaskDataFrame for elements (only the SpatialData overload returns SpatialData).
+    centroids = cast(DaskDataFrame, get_centroids(element, coordinate_system=INTRINSIC_COORDINATE_SYSTEM)).compute()
     del d[INTRINSIC_COORDINATE_SYSTEM]
     return centroids
 
diff --git a/src/spatialdata/_core/spatialdata.py b/src/spatialdata/_core/spatialdata.py
index fb55ab08..6a168bf4 100644
--- a/src/spatialdata/_core/spatialdata.py
+++ b/src/spatialdata/_core/spatialdata.py
@@ -543,6 +543,39 @@ def aggregate(
             **kwargs,
         )
 
+    def get_centroids(
+        self,
+        element_name: str,
+        coordinate_system: str | None = "global",
+        return_background: bool = False,
+        return_area: bool = False,
+        persist_as: Literal["Points", "adata"] = "Points",
+        table_name: str | None = None,
+        inplace: bool = True,
+        overwrite: bool = False,
+    ) -> DaskDataFrame | AnnData | None:
+        """Get the centroids of ``element_name``, or persist them into its annotating table.
+
+        Convenience method for :func:`spatialdata.get_centroids` called on ``self``; see that function for the
+        complete docstring. With ``persist_as="adata"`` the centroids are written into ``obsm["spatial"]`` (and area
+        into ``obs["area"]``) of the resolved annotating table; ``inplace=True`` mutates it and returns ``None``,
+        ``inplace=False`` returns a modified copy of that table. ``overwrite=True`` allows replacing pre-existing
+        values in those keys for the element's rows (refused by default).
+        """
+        from spatialdata._core.centroids import _get_centroids_sdata
+
+        return _get_centroids_sdata(
+            self,
+            element_name,
+            coordinate_system=coordinate_system,
+            return_background=return_background,
+            return_area=return_area,
+            persist_as=persist_as,
+            table_name=table_name,
+            inplace=inplace,
+            overwrite=overwrite,
+        )
+
     def is_backed(self) -> bool:
         """Check if the data is backed by a Zarr storage or if it is in-memory."""
         return self.path is not None
diff --git a/src/spatialdata/dataloader/datasets.py b/src/spatialdata/dataloader/datasets.py
index 03879abc..07abacf6 100644
--- a/src/spatialdata/dataloader/datasets.py
+++ b/src/spatialdata/dataloader/datasets.py
@@ -5,12 +5,13 @@
 from functools import partial
 from itertools import chain
 from types import MappingProxyType
-from typing import Any
+from typing import Any, cast
 
 import anndata as ad
 import numpy as np
 import pandas as pd
 from anndata import AnnData
+from dask.dataframe import DataFrame as DaskDataFrame
 from geopandas import GeoDataFrame
 from pandas import CategoricalDtype
 from scipy.sparse import issparse
@@ -494,7 +495,7 @@ def _get_tile_coords(
 
     # extent, aka the tile size
     extent = (circles.radius * 2).values.reshape(-1, 1)
-    centroids_points = get_centroids(circles, coordinate_system=cs)
+    centroids_points = cast(DaskDataFrame, get_centroids(circles, coordinate_system=cs))
     axes = get_axes_names(centroids_points)
     centroids_numpy = centroids_points.compute().values
 
diff --git a/src/spatialdata/models/models.py b/src/spatialdata/models/models.py
index a818bdad..d6302154 100644
--- a/src/spatialdata/models/models.py
+++ b/src/spatialdata/models/models.py
@@ -950,6 +950,9 @@ class TableModel:
     REGION_KEY_KEY = "region_key"
     INSTANCE_KEY = "instance_key"
     ATTRS_KEY = ATTRS_KEY
+    # squidpy-style storage keys for centroids/area persisted into a table (see get_centroids).
+    SPATIAL_KEY = "spatial"
+    AREA_KEY = "area"
 
     @classmethod
     def _validate_set_region_key(cls, data: AnnData, region_key: str | None = None) -> None:
diff --git a/tests/core/test_centroids.py b/tests/core/test_centroids.py
index f8c8be1d..3ad08f92 100644
--- a/tests/core/test_centroids.py
+++ b/tests/core/test_centroids.py
@@ -11,11 +11,20 @@
 from spatialdata._core.centroids import get_centroids
 from spatialdata._core.query.relational_query import get_element_instances
 from spatialdata.models import Labels2DModel, Labels3DModel, PointsModel, TableModel, get_axes_names
-from spatialdata.transformations import Affine, Identity, get_transformation, set_transformation
+from spatialdata.transformations import Affine, Identity, Scale, get_transformation, set_transformation
 
 RNG = default_rng(42)
 
 
+def _assert_obsm_matches_points(table: AnnData, pts: pd.DataFrame) -> None:
+    # written obsm["spatial"] must match the element-level Points centroids on shared (non-background) ids.
+    inst = table.obs["instance_id"].to_numpy()
+    written = pd.DataFrame(table.obsm["spatial"], index=inst, columns=["x", "y"])
+    common = pts.index.intersection(written.index[inst != 0])
+    assert len(common) > 0
+    assert np.allclose(written.loc[common].to_numpy(), pts.loc[common][["x", "y"]].to_numpy())
+
+
 def _get_affine() -> Affine:
     theta: float = math.pi / 18
     k = 10.0
@@ -172,9 +181,150 @@ def test_get_centroids_labels(
         assert np.allclose(centroids.compute().values, centroids_transformed)
 
 
+def test_get_centroids_labels_area(labels):
+    # area for labels is the per-label pixel count; it rides along as a feature column of the Points.
+    element = labels["labels2d"]
+    centroids = get_centroids(element, return_area=True)
+    assert "area" in centroids.columns
+    ids, counts = np.unique(np.asarray(element.data), return_counts=True)
+    expected = dict(zip(ids, counts, strict=True))
+    got = centroids[["area"]].compute()["area"]
+    assert not (got.index == 0).any()  # background dropped
+    for label_id, area in got.items():
+        assert area == expected[label_id]
+
+
+def test_get_centroids_shapes_area_circles(shapes):
+    element = shapes["circles"]
+    centroids = get_centroids(element, return_area=True)
+    expected = np.pi * np.asarray(element["radius"], dtype=float) ** 2
+    assert np.allclose(centroids["area"].compute().to_numpy(), expected)
+
+
+@pytest.mark.parametrize("shapes_name", ["poly", "multipoly"])
+def test_get_centroids_shapes_area_polygons(shapes, shapes_name: str):
+    element = shapes[shapes_name]
+    centroids = get_centroids(element, return_area=True)
+    assert np.allclose(centroids["area"].compute().to_numpy(), element.geometry.area.to_numpy())
+
+
+def test_get_centroids_points_area_raises(points):
+    with pytest.raises(ValueError, match="not supported for points"):
+        get_centroids(points["points_0"], return_area=True)
+
+
+def test_get_centroids_area_non_identity_transform_raises(full_sdata):
+    # persist_as='Points' transforms the centroids but the area is intrinsic (untransformed); requesting
+    # both under a non-identity transform must raise rather than mix coordinate systems.
+    set_transformation(full_sdata["labels2d"], affine, "aligned")
+    with pytest.raises(ValueError, match="requires intrinsic coordinates"):
+        get_centroids(full_sdata["labels2d"], coordinate_system="aligned", return_area=True)
+
+
+def test_get_centroids_element_persist_adata_raises(labels):
+    # an element on its own has no annotating table; persist_as='adata' needs the SpatialData.
+    with pytest.raises(ValueError, match="persist_as='adata'"):
+        get_centroids(labels["labels2d"], persist_as="adata")
+
+
+def test_get_centroids_sdata_persist_into_table(full_sdata):
+    # `table` has 100 rows annotating `labels2d`; the `instance_id` column holds the label ids each row
+    # refers to (the value at a given row position is the label, not the row's position). The background
+    # label 0 has no centroid, so its row stays NaN.
+    table = full_sdata["table"]
+    assert "spatial" not in table.obsm
+    out = get_centroids(full_sdata, "labels2d", coordinate_system="global", return_area=True, persist_as="adata")
+    assert out is None  # inplace=True (default) mutates the table and returns nothing
+
+    table = full_sdata["table"]
+    assert table.obsm["spatial"].shape == (table.n_obs, 2)
+    assert "area" in table.obs
+
+    inst = table.obs["instance_id"].to_numpy()
+    finite = np.isfinite(table.obsm["spatial"]).all(axis=1)
+    assert finite[inst != 0].all()  # every non-background label got a centroid
+    assert not finite[inst == 0].any()  # background row stays NaN
+
+    # coordinates must match the element-level Points (global transform is the identity here)
+    pts = get_centroids(full_sdata["labels2d"], coordinate_system="global").compute()
+    _assert_obsm_matches_points(table, pts)
+
+    # area must equal the pixel counts of the corresponding labels
+    ids, counts = np.unique(np.asarray(full_sdata["labels2d"].data), return_counts=True)
+    count_of = dict(zip(ids, counts, strict=True))
+    area = table.obs["area"].to_numpy()
+    for row, label_id in enumerate(inst):
+        if label_id != 0:
+            assert area[row] == count_of[label_id]
+
+
+def test_get_centroids_sdata_persist_rejects_non_identity(full_sdata):
+    # obsm["spatial"] stores intrinsic coordinates only (we don't persist the coordinate system), so
+    # persist_as='adata' refuses a non-identity transform instead of silently storing transformed coords.
+    set_transformation(full_sdata["labels2d"], affine, "aligned")
+    with pytest.raises(ValueError, match="intrinsic coordinates"):
+        get_centroids(full_sdata, "labels2d", coordinate_system="aligned", persist_as="adata")
+
+
+def test_get_centroids_sdata_persist_intrinsic_matches_identity(full_sdata):
+    # coordinate_system=None (intrinsic) equals a coordinate system whose transform is the identity.
+    get_centroids(full_sdata, "labels2d", coordinate_system="global", persist_as="adata")
+    global_spatial = full_sdata["table"].obsm["spatial"].copy()
+    # overwrite=True because the first write already populated these rows.
+    get_centroids(full_sdata, "labels2d", coordinate_system=None, persist_as="adata", overwrite=True)
+    intrinsic_spatial = full_sdata["table"].obsm["spatial"]
+    finite = np.isfinite(global_spatial).all(axis=1)
+    assert np.allclose(global_spatial[finite], intrinsic_spatial[finite])
+
+
+def test_get_centroids_sdata_persist_inplace_false_returns_copy(full_sdata):
+    # inplace=False copies only the target table, writes into the copy, and leaves the sdata untouched.
+    out = get_centroids(full_sdata, "labels2d", return_area=True, persist_as="adata", inplace=False)
+    assert isinstance(out, AnnData)
+    assert out is not full_sdata["table"]
+    assert "spatial" in out.obsm and "area" in out.obs
+    assert "spatial" not in full_sdata["table"].obsm  # original table not modified
+
+
+def test_spatialdata_get_centroids_method(full_sdata):
+    # the method mirrors the module-level function for both persistence modes.
+    pts = full_sdata.get_centroids("labels2d", coordinate_system="global")
+    expected = get_centroids(full_sdata["labels2d"], coordinate_system="global")
+    assert np.allclose(pts.compute().to_numpy(), expected.compute().to_numpy())
+
+    assert full_sdata.get_centroids("labels2d", persist_as="adata") is None
+    assert "spatial" in full_sdata["table"].obsm
+
+
+def test_get_centroids_sdata_persist_instance_key_mismatch_raises(full_sdata):
+    # instance ids stored with a dtype that doesn't match the element's integer labels must fail
+    # loudly instead of silently writing NaN coordinates into obsm["spatial"].
+    table = full_sdata["table"]
+    table.obs["instance_id"] = table.obs["instance_id"].astype(str)
+    with pytest.raises(ValueError, match="No instance id annotating"):
+        get_centroids(full_sdata, "labels2d", persist_as="adata")
+
+
+def test_get_centroids_sdata_persist_refuses_dim_mismatch(full_sdata):
+    # This mismatch arises when one table annotates multiple elements of different dimensionality: e.g. a
+    # 3D element already wrote width-3 coordinates into obsm["spatial"], and we now persist 2D centroids
+    # for `labels2d`. The existing width-3 array (simulated here) must not be silently overwritten, as
+    # that would wipe the other element's coordinates sharing the table.
+    table = full_sdata["table"]
+    table.obsm["spatial"] = np.zeros((table.n_obs, 3))
+    with pytest.raises(ValueError, match="refusing to overwrite"):
+        get_centroids(full_sdata, "labels2d", persist_as="adata")
+
+
+def test_get_centroids_sdata_no_table_raises(full_sdata):
+    # points_0 is not annotated by any table -> the error points the user to persist_as='Points'.
+    with pytest.raises(ValueError, match="persist_as='Points'"):
+        get_centroids(full_sdata, "points_0", persist_as="adata")
+
+
 def test_get_centroids_invalid_element(images):
     # cannot compute centroids for images
-    with pytest.raises(ValueError, match="Expected a `Labels` element. Found an `Image` instead."):
+    with pytest.raises(ValueError, match="Centroids are not supported for elements modeled by Image2DModel"):
         get_centroids(images["image2d"])
 
     # cannot compute centroids for tables
@@ -190,5 +340,40 @@ def test_get_centroids_invalid_element(images):
 
 
 def test_get_centroids_invalid_coordinate_system(points):
-    with pytest.raises(AssertionError, match="No transformation to coordinate system"):
+    with pytest.raises(ValueError, match="No transformation to coordinate system"):
         get_centroids(points["points_0"], coordinate_system="invalid")
+
+
+def test_get_centroids_numerically_identity_transform_is_intrinsic(full_sdata):
+    # a mathematically-identity transform that is not the `Identity` class (e.g. Scale([1, 1])) must be
+    # accepted as intrinsic, both for return_area and for persist_as='adata'.
+    set_transformation(full_sdata["labels2d"], Scale([1.0, 1.0], axes=("x", "y")), "noop")
+    # persist_as='adata' must not raise "cannot apply the non-identity transform"
+    get_centroids(full_sdata, "labels2d", coordinate_system="noop", persist_as="adata")
+    assert "spatial" in full_sdata["table"].obsm
+    # return_area on the Points path must not raise "requires intrinsic coordinates"
+    out = get_centroids(full_sdata["labels2d"], coordinate_system="noop", return_area=True)
+    assert "area" in out.columns
+
+
+def test_get_centroids_points_area_non_identity_reports_points_not_intrinsic(points):
+    # a Points element with a non-identity transform + return_area must report the points-specific error,
+    # not the misleading "requires intrinsic coordinates" message.
+    set_transformation(points["points_0"], affine, "aligned")
+    with pytest.raises(ValueError, match="not supported for points"):
+        get_centroids(points["points_0"], coordinate_system="aligned", return_area=True)
+
+
+def test_get_centroids_sdata_persist_overwrite_guard(full_sdata):
+    # a pre-existing obs['area'] (possibly unrelated) must not be silently clobbered; overwrite=True allows it.
+    table = full_sdata["table"]
+    table.obs["area"] = np.arange(table.n_obs, dtype=float)
+    with pytest.raises(ValueError, match="overwrite=True"):
+        get_centroids(full_sdata, "labels2d", return_area=True, persist_as="adata")
+    # with overwrite=True it proceeds and writes the pixel-count area
+    get_centroids(full_sdata, "labels2d", return_area=True, persist_as="adata", overwrite=True)
+    assert "spatial" in table.obsm
+
+    # obsm["spatial"] is likewise guarded once populated
+    with pytest.raises(ValueError, match="overwrite=True"):
+        get_centroids(full_sdata, "labels2d", persist_as="adata")