@@ -4453,57 +4453,72 @@ def _stream_label_centroid_stats(data: Any) -> tuple[ArrayLike, ArrayLike, Array
 
     Streams the raster block by block — one chunk in memory at a time for a dask array, a
     bounded row-block at a time for a numpy array — accumulating per-label ``count`` (= area),
-    ``sum_x`` and ``sum_y``. The reduction is additive, so it is exact across block boundaries
-    and uses only O(n_labels) memory regardless of raster size. Background label 0 is excluded.
+    ``sum_x`` and ``sum_y``. Labels are relabelled to a dense ``0..k-1`` range, so memory is
+    O(number of distinct labels), independent of the raster size *and* of the label-id magnitude
+    (sparse/global ids do not blow it up). The reduction is additive, so it is exact across block
+    boundaries. Background label 0 is excluded.
     """
     n_rows, n_cols = data.shape
-    is_dask = hasattr(data, "chunks") and not isinstance(data, np.ndarray)
-    if is_dask:
-        n_labels = int(data.max().compute()) + 1
+    if hasattr(data, "chunks"):  # dask
         block_slices = slices_from_chunks(data.chunks)
     else:
         data = np.asarray(data)
-        n_labels = int(data.max()) + 1
         # bound the per-block coordinate-weight arrays to ~8M pixels
         step = max(1, min(n_rows, (8 << 20) // max(1, n_cols)))
         block_slices = [(slice(r0, min(r0 + step, n_rows)), slice(0, n_cols)) for r0 in range(0, n_rows, step)]
 
-    count = np.zeros(n_labels, dtype=np.float64)
-    sum_x = np.zeros(n_labels, dtype=np.float64)
-    sum_y = np.zeros(n_labels, dtype=np.float64)
-    for row_sl, col_sl in block_slices:
+    def _load(row_sl: slice, col_sl: slice) -> np.ndarray:
         block = data[row_sl, col_sl]
         block = np.asarray(block.compute() if hasattr(block, "compute") else block)
+        # label rasters are integer-valued even when stored as float; cast so np.unique/searchsorted
+        # stay integer (and the dense bincount indices below are always int).
+        return block.astype(np.int64) if block.dtype.kind == "f" else block
+
+    # Pass 1: the sorted set of present label values -> dense relabelling (keeps memory O(n_labels)).
+    uniq = np.zeros(0, dtype=np.int64)
+    for row_sl, col_sl in block_slices:
+        uniq = np.union1d(uniq, np.unique(_load(row_sl, col_sl)))
+    k = uniq.size
+
+    # Pass 2: additive per-(dense-)label count / sum_x / sum_y.
+    count = np.zeros(k)
+    sum_x = np.zeros(k)
+    sum_y = np.zeros(k)
+    for row_sl, col_sl in block_slices:
+        block = _load(row_sl, col_sl)
         block_rows, block_cols = block.shape
-        flat = block.reshape(-1)
+        idx = np.searchsorted(uniq, block.reshape(-1))  # dense 0..k-1 indices (always int)
         cols = np.tile(np.arange(col_sl.start, col_sl.start + block_cols, dtype=np.float64), block_rows)
         rows = np.repeat(np.arange(row_sl.start, row_sl.start + block_rows, dtype=np.float64), block_cols)
-        count += np.bincount(flat, minlength=n_labels)
-        sum_x += np.bincount(flat, weights=cols, minlength=n_labels)
-        sum_y += np.bincount(flat, weights=rows, minlength=n_labels)
+        count += np.bincount(idx, minlength=k)
+        sum_x += np.bincount(idx, weights=cols, minlength=k)
+        sum_y += np.bincount(idx, weights=rows, minlength=k)
 
-    labels = np.flatnonzero(count)
-    labels = labels[labels != 0]  # drop background and absent labels
-    return labels, sum_x[labels] / count[labels], sum_y[labels] / count[labels], count[labels]
+    keep = uniq != 0  # drop background; every kept label has count >= 1
+    return uniq[keep], sum_x[keep] / count[keep], sum_y[keep] / count[keep], count[keep]
 
 
 def _compute_element_measurements(sdata: SpatialData, element_name: str) -> pd.DataFrame:
     """One row per instance with intrinsic ``["x", "y", "area"]``, indexed by instance id.
 
-    Shapes use shapely's vectorized centroid; circles (``Point`` + ``radius``) have ``area = pi*r**2``
-    (shapely ``.area`` is 0 for them), polygons use the geometric area. 2D labels use the streaming
-    bincount aggregator (``area`` = pixel count) — it holds one chunk plus O(n_labels) accumulators,
-    so it scales to Xenium-size masks where a whole-array ``regionprops`` table would run out of
-    memory. Area meaning differs across element types but is consistent within one element.
+    Shapes use shapely's vectorized centroid; circles (``Point`` geometry + ``radius``) have
+    ``area = pi*r**2`` (shapely ``.area`` is 0 for them), polygons use the geometric area. 2D labels
+    use the streaming bincount aggregator (``area`` = pixel count) — it holds one chunk plus
+    O(n_labels) accumulators, so it scales to Xenium-size masks where a whole-array ``regionprops``
+    table would run out of memory. Area meaning differs across element types but is consistent within
+    one element.
     """
     element = sdata[element_name]
     model = get_model(element)
     if model is ShapesModel:
-        centroids = element.geometry.centroid
-        if "radius" in element.columns:  # circles are Points; shapely .area is 0, so use pi * r**2
+        geometry = element.geometry
+        centroids = geometry.centroid
+        # Dispatch on geometry TYPE, not a column name: circles are Point geometries (shapely .area
+        # is 0 for them) with a radius -> pi*r**2; everything else uses the true geometric area.
+        if (geometry.geom_type == "Point").all():
             area = np.pi * np.asarray(element["radius"], dtype=float) ** 2
         else:
-            area = element.geometry.area.to_numpy()
+            area = geometry.area.to_numpy()
         xy = {"x": centroids.x.to_numpy(), "y": centroids.y.to_numpy(), "area": area}
         return pd.DataFrame(xy, index=element.index)
     if model is Labels2DModel:
@@ -4535,6 +4550,26 @@ def _valid_spatial_obsm(arr: ArrayLike, n_obs: int) -> bool:
     return bool(arr.ndim == 2 and arr.shape == (n_obs, 2))
 
 
+def _obsm_region_finite(table: AnnData, key: str, mask: ArrayLike) -> bool:
+    """Whether ``obsm[key]`` already holds finite coords for every ``mask`` row (already populated)."""
+    if key not in table.obsm:
+        return False
+    arr = np.asarray(table.obsm[key])
+    if not _valid_spatial_obsm(arr, table.n_obs):
+        return False
+    region = arr[mask].astype(float)
+    return bool(region.size and np.isfinite(region).all())
+
+
+def _check_obs_numeric(table: AnnData, key: str) -> None:
+    """Raise if ``obs[key]`` exists but is non-numeric, before any mutation (avoids half-writes)."""
+    if key in table.obs and not is_numeric_dtype(table.obs[key]):
+        raise ValueError(
+            f"Cannot write measurements into obs[{key!r}]: the existing column is "
+            f"{table.obs[key].dtype} (not numeric). Pass a different key or drop the column first."
+        )
+
+
 def _write_obs_region(table: AnnData, mask: ArrayLike, key: str, values: ArrayLike) -> None:
     """Write ``values`` into ``obs[key]`` at ``mask`` rows; other rows stay/are NaN."""
     if key in table.obs:
@@ -4579,14 +4614,46 @@ def _measure_into_table(
 ) -> None:
     """Compute and write the requested measurements for one element into its annotating table.
 
-    Overwrites existing values for the element's rows. A table may annotate several elements, so
-    only the rows belonging to ``element_name`` are touched.
+    Only the rows belonging to ``element_name`` are touched (a table may annotate several elements).
+    Centroids already present for those rows (e.g. reader-provided ``obsm["spatial"]``) are not
+    overwritten; ``area``/``diameter`` overwrite our own columns. All targets are validated before
+    the first write, so a bad column never leaves the table half-written.
     """
     table = sdata.tables[table_name]
     instance_key, mask = _region_mask_and_keys(table, element_name)
     if not mask.any():
         raise ValueError(f"Table {table_name!r} does not annotate element {element_name!r} (no matching rows).")
-    meas = _compute_element_measurements(sdata, element_name).reindex(table.obs[instance_key].to_numpy()[mask])
+
+    # #1: never clobber centroids already populated for this element's rows (reader/prior-call coords).
+    if centroids and _obsm_region_finite(table, obsm_key, mask):
+        warnings.warn(
+            f"obsm[{obsm_key!r}] is already populated for element {element_name!r}; not overwriting "
+            f"its centroids. Remove it or use a different `obsm_key` to recompute.",
+            UserWarning,
+            stacklevel=3,
+        )
+        centroids = False
+    if not (centroids or area or diameter):
+        return
+
+    keys = table.obs[instance_key].to_numpy()[mask]
+    meas = _compute_element_measurements(sdata, element_name).reindex(keys)
+    # #2: instance ids annotated in the table but absent from the element reindex to NaN -> warn.
+    missing = int(meas[["x", "y"]].isna().any(axis=1).sum())
+    if missing:
+        warnings.warn(
+            f"{missing}/{len(keys)} instances annotated for {element_name!r} have no match in the "
+            f"element (instance-id dtype mismatch, e.g. str vs int?); writing NaN for them.",
+            UserWarning,
+            stacklevel=3,
+        )
+
+    # #4: validate every obs target up front so an existing non-numeric column raises before any write.
+    if area:
+        _check_obs_numeric(table, area_key)
+    if diameter:
+        _check_obs_numeric(table, diameter_key)
+
     if centroids:
         _write_obsm_region(table, mask, obsm_key, meas[["x", "y"]].to_numpy())
     if area:
@@ -4645,16 +4712,18 @@ def measure_obs(
     element and writes them, squidpy-style, into the annotating :class:`~anndata.AnnData` table:
     centroids go to ``table.obsm[obsm_key]`` (an ``(n_obs, 2)`` array, the canonical
     ``obsm["spatial"]``), area and diameter to ``table.obs``. Values are stored in the element's
-    *intrinsic* coordinates/units (coordinate-system independent). Labels area is the pixel count;
-    shapes area is ``geometry.area`` (``pi*r**2`` for circles); equivalent diameter is
-    ``2 * sqrt(area / pi)``. Persisting them once lets later renders (and downstream tools such as
-    squidpy) reuse them instead of recomputing.
-
-    Existing values for the measured rows are overwritten; pass ``centroids=False`` to keep a
-    pre-existing ``obsm[obsm_key]`` (e.g. from a reader). Per-cell measurements need a table to
-    write into, so an element without an annotating table cannot be measured (this raises). The
-    label path never densifies the raster — it streams it block by block, scaling to Xenium-size
-    masks.
+    *intrinsic* pixel coordinates/units (which align directly with the element's own raster/geometry).
+    Labels area is the pixel count; shapes area is ``geometry.area`` (``pi*r**2`` for circles);
+    equivalent diameter is ``2 * sqrt(area / pi)``. Persisting them once lets later renders (and
+    downstream tools such as squidpy) reuse them instead of recomputing.
+
+    Centroids already present for an element's rows (e.g. a reader-provided ``obsm[obsm_key]``) are
+    **not** overwritten — a warning is emitted and that element's centroid write is skipped (remove
+    the key or use a different ``obsm_key`` to recompute); ``area``/``diameter`` overwrite our own
+    columns. Instances annotated in the table but absent from the element are written as NaN with a
+    warning. Per-cell measurements need a table to write into, so an element without an annotating
+    table cannot be measured (this raises). The label path never densifies the raster — it streams
+    it block by block with memory O(n_labels), scaling to Xenium-size masks.
 
     Parameters
     ----------
 
@@ -9,7 +9,7 @@
 from anndata import AnnData
 from shapely.geometry import Point
 from spatialdata import SpatialData, get_centroids
-from spatialdata.models import PointsModel, ShapesModel, TableModel
+from spatialdata.models import Labels2DModel, PointsModel, ShapesModel, TableModel
 
 import spatialdata_plot
 from spatialdata_plot.pl import measure_obs
@@ -511,6 +511,19 @@ def _add_shapes_table(sdata: SpatialData, element: str = "blobs_polygons", name:
     return sdata
 
 
+def _labels_sdata(arr: np.ndarray, name: str = "lab", table: str = "t") -> SpatialData:
+    """Build a SpatialData with a single 2D-labels element annotated by a table."""
+    ids = np.unique(arr)
+    ids = ids[ids != 0].astype(int)
+    adata = AnnData(np.zeros((len(ids), 1), dtype=np.float32))
+    adata.obs["instance_id"] = ids
+    adata.obs["region"] = name
+    return SpatialData(
+        labels={name: Labels2DModel.parse(arr, dims=("y", "x"))},
+        tables={table: TableModel.parse(adata, region_key="region", instance_key="instance_id", region=name)},
+    )
+
+
 class TestMeasureObs:
     """`measure_obs` materializes centroids/area/equivalent diameter into the annotating table."""
 
@@ -572,12 +585,15 @@ def test_inplace_false_leaves_original_untouched(self, sdata_blobs: SpatialData)
         assert "spatial" in out["table"].obsm
         assert "spatial" not in sdata_blobs["table"].obsm  # original not mutated
 
-    def test_recompute_overwrites(self, sdata_blobs: SpatialData) -> None:
-        measure_obs(sdata_blobs, "blobs_labels")
+    def test_existing_centroids_not_clobbered(self, sdata_blobs: SpatialData) -> None:
+        # #1: a populated obsm["spatial"] (reader- or prior-call-provided) is not overwritten.
         table = sdata_blobs["table"]
-        table.obsm["spatial"][0] = [999.0, 999.0]  # corrupt one row
-        measure_obs(sdata_blobs, "blobs_labels")  # second call overwrites with the real centroid
-        assert tuple(table.obsm["spatial"][0]) != (999.0, 999.0)
+        sentinel = np.arange(table.n_obs * 2, dtype=float).reshape(table.n_obs, 2)
+        table.obsm["spatial"] = sentinel.copy()
+        with pytest.warns(UserWarning, match="already populated"):
+            measure_obs(sdata_blobs, "blobs_labels")
+        np.testing.assert_array_equal(table.obsm["spatial"], sentinel)  # centroids untouched
+        assert "area" in table.obs  # area/diameter still written
 
     def test_centroids_false_keeps_existing_obsm(self, sdata_blobs: SpatialData) -> None:
         table = sdata_blobs["table"]
@@ -593,6 +609,51 @@ def test_incompatible_obsm_shape_raises(self, sdata_blobs: SpatialData) -> None:
         with pytest.raises(ValueError, match="Refusing to overwrite"):
             measure_obs(sdata_blobs, "blobs_labels")
 
+    def test_unmatched_instance_ids_warn_and_write_nan(self, sdata_blobs: SpatialData) -> None:
+        # #2: instance ids that don't match the element (e.g. str vs int) -> warn + NaN, not silent.
+        table = sdata_blobs["table"]
+        table.obs["instance_id"] = table.obs["instance_id"].astype(str)
+        with pytest.warns(UserWarning, match="no match"):
+            measure_obs(sdata_blobs, "blobs_labels")
+        assert np.isnan(table.obsm["spatial"]).all()
+
+    def test_float_dtype_labels_supported(self, sdata_blobs: SpatialData) -> None:
+        # #3: a float-typed (but integer-valued) labels raster must not crash np.bincount.
+        arr = np.asarray(sdata_blobs["blobs_labels"].data).astype(np.float32)
+        sd = _labels_sdata(arr)
+        measure_obs(sd, "lab", table_name="t")
+        assert np.isfinite(sd["t"].obsm["spatial"]).all()
+
+    def test_existing_nonnumeric_column_raises_before_any_write(self, sdata_blobs: SpatialData) -> None:
+        # #4: a non-numeric collision raises BEFORE obsm is mutated (no half-written table).
+        table = sdata_blobs["table"]
+        table.obs["area"] = pd.Categorical(["lo"] * table.n_obs)
+        with pytest.raises(ValueError, match="not numeric"):
+            measure_obs(sdata_blobs, "blobs_labels")
+        assert "spatial" not in table.obsm  # atomic: nothing written
+
+    def test_sparse_high_label_ids(self, sdata_blobs: SpatialData) -> None:
+        # #5: sparse/high label ids (max id >> n_labels) are measured correctly (dense relabelling).
+        arr = np.asarray(sdata_blobs["blobs_labels"].data)
+        hi = (arr.astype(np.int64) * 1000)  # ids become 1000, 2000, ... ; max id is huge, few labels
+        measure_obs(sd_hi := _labels_sdata(hi), "lab", table_name="t")
+        measure_obs(sd_lo := _labels_sdata(arr.astype(np.int64)), "lab", table_name="t")
+        # relabelling values does not move pixels -> identical centroid set
+        np.testing.assert_allclose(
+            np.sort(sd_hi["t"].obsm["spatial"], axis=0), np.sort(sd_lo["t"].obsm["spatial"], axis=0)
+        )
+
+    def test_polygon_with_radius_column_uses_geometric_area(self, sdata_blobs: SpatialData) -> None:
+        # #7: dispatch on geometry type, not the "radius" column name -> polygons use geometry.area.
+        gdf = sdata_blobs["blobs_polygons"].copy()
+        gdf["radius"] = 5.0  # incidental column; must NOT trigger the circle (pi*r**2) branch
+        sdata_blobs["pr"] = ShapesModel.parse(gdf)
+        _add_shapes_table(sdata_blobs, "pr", name="pt")
+        measure_obs(sdata_blobs, "pr", table_name="pt")
+        order = sdata_blobs["pt"].obs["instance_id"].to_numpy()
+        expected = gdf.geometry.area.to_numpy()[[list(gdf.index).index(i) for i in order]]
+        np.testing.assert_allclose(sdata_blobs["pt"].obs["area"].to_numpy(), expected, rtol=1e-9)
+
     def test_flags_select_outputs(self, sdata_blobs: SpatialData) -> None:
         measure_obs(sdata_blobs, "blobs_labels", area=False, diameter=False)
         table = sdata_blobs["table"]