Unexpected Coordinate Values with Image2DModel

Help spatialdata
1

Hi, I’ve been learning spatialdata for a couple of days now and I’m trying to apply it to some MIBI data represented as a collection of n single channel tiffs.

Currently for each FOV we store it’s channels as individual TIFF files (512 \times 512), for example:

mibi_cohort/
└── fov0/
    ├── chan_0.tiff
    ├── chan_1.tiff
    ├── chan_2.tiff
    ├── ...
    └── chan_n.tiff

When loading the data, I’m loading all the images together with tifffile and reading / opening them as a Zarr Store / Array. Then I convert it to an Xarray.DataArray, and then I convert it to a SpatialData object with Image2DModel.

This all works, however the x and y coordinates are now floating points starting at 0.5 to 511.5 incrementing by 0.5.

Here is a minimum working example

# Imports

from pathlib import Path

import spatialdata as sd
from spatialdata.models import Image2DModel, C, Y, X

import tifffile
import numpy as np
import zarr

import xarray as xr

# Creating some example single channel TIFFS
ex_fov_path = Path("./data/fov_ex")
ex_fov_path.mkdir(exist_ok=True, parents=True)

rng = np.random.default_rng()

rng_data = rng.random(size=(3, 512, 512))

channels = ["chan_0", "chan_1", "chan_2"]

# Populate an array (Channels, Y, X) with random data and save as TIFFs

for data, channel in zip(rng_data, channels):
    tifffile.imwrite(ex_fov_path / f"{channel}.tiff", data=data)


# 1. Load the tiffs
# 2. Convert to `Zarr` Store and the open `Zarr` Store, 
# 3. Convert to `xr.DataArray`
# 4. Convert to `spatialdata.models.Image2DModel`
# 5. Convert to SpatialData Object

def load_folder_of_tiffs(p: Path):
    tiff_paths = list(p.glob("*.tiff"))
    channel_names: list[str] = [tp.stem for tp in tiff_paths]
    zarr_store: tifffile.ZarrFileSequenceStore = tifffile.imread(
        tiff_paths, aszarr=True
    )
    zarr_array = zarr.open(zarr_store)

    fov = xr.DataArray(
        data=zarr_array,
        dims=(C, Y, X),
        coords={
            C: channel_names,
            Y: np.arange(512, dtype=np.int16),
            X: np.arange(512, dtype=np.int16),
        },
    )

    fov_sd = sd.SpatialData(
        images={"fov": Image2DModel.parse(data=fov, scale_factors=None)}
    )

    return fov_sd

Run the function with the example directory

fov = load_folder_of_tiffs(ex_fov_path)
fov.images["fov"]

Here is the html rendered output of the DataArray.

Example Notebook Output
Example Notebook Output724×408 16.6 KB

I expected the x and y coordinates to be integers covering [0,511].

Thank you!

2

Thank you @srivarra for reaching out and for the reproducible example. The behavior you encountered is expected and I will elaborate on the reasons, how to circumvent it and on the limitations.

TLDR; if you need positional indexing you can use fov.isel(y=1) or fov.isel(y=slice(1, 2)).

Saving and reading

In SpatialData we use the NGFF specification to save and load raster tensors to Zarr. The in-memory representation of this Zarr storage is similar to what xarray.DataArray offers, but with some differences. One of these differences is that NGFF doesn’t save the z, y, x coordinates explicitly (it saves only c) and it has a more laborious (but more powerful) way to deal with non-regular grids, which includes also more complex warps (see specification here). One note: the specification that I linked is still in draft mode and we currently don’t support non-linear transformations, but we have one master student, Tobias, that started working on this.

Considering the above, when we save an xarray.DataArray object to disk we discard the z, y, x coordinates, and we reconstruct them upon loading the object. By design in doing this we reconstruct uniformly spaced coordinates.

Centered vs non-centered coordinates

Before even explaining why we use centered coordinates as you saw, here are two limitations of this approach:

  1. if you have coordinates that are non-linearly spaced, like 0, 1, 10, these will be lost and there is no way to save and load them for the moment.
  2. if you have two tensors with shape 1x3x3 (1 channel image) with uniformly spaced coordinates, one with 0-based coordinates 0, 1, 2, and one with centered coordinates 0.5, 1.5, 2.5, when you save them and load them they will be reconstructed in only one way. As I said currently we use centered coordinates. So in both cases you will get 0.5., 1.5, 2.5.

The reason why we opted for centered coordinates is for when we deal with multiscale images, which are objects that are represented with xarray.DataTree and store multiple xarray.DataArray images with decreasing resolutions and aligned coordinates.

Say that you have an image with resolution 999x999 spanning between the coordinates 0 and 1000. With centered coordinates the xarray coordinates of the image are 0.5, 1.5, … 999.5. If you parse the image with ImageModel.parse(image, scale_factors=[2, 2, 2, 2]) (which leads to 5 total scales, the smallest having a 16x downscale), the smallest scale will have this coordinates

└── DataTree('scale4')
        Dimensions:  (c: 3, y: 62, x: 62)
        Coordinates:
          * c        (c) int64 0 1 2
          * y        (y) float64 8.056 24.17 40.28 56.4 ... 942.6 958.7 974.8 990.9
          * x        (x) float64 8.056 24.17 40.28 56.4 ... 942.6 958.7 974.8 990.9
        Data variables:
            image    (c, y, x) float64 dask.array<chunksize=(3, 62, 62), meta=np.ndarray>

Here, if you do some slicing operations based on physical coordinates, having non-centered coordinates would lead to some asymmetry in what is selected and what not.

Also if you suppose to do operations like queries (subset) or aggregations (e.g. summing pixels by raster masks), I think it’s natural to interpret the weigh of the pixel as in the center of the pixel. This comes automatic if we use centered coordinates.

Positional slicing

Please notice that the above implementation doesn’t impede the users from doing positional based indexing slicing operations for instance fov.isel(y=1), which would be equivalent to
fov.sel(y=1.5).

Similarly one can use fov.sel(y=slice(2, 3)), which is equivalent to the more natural fov.sel(y=(1, 2)).

Anyway when dealing with physical coordinates, sel (non-positional) indexing is probably more common.

Asking for feedback

To conclude, I want to anticipate that soon we will refactor the way coordinate transformations interact with xarray coordinates, improving the ergonomics of the library. For instance if one applies a scale operation to the coordinates (i.e. only the coordinates change but the shape of the tensor doesn’t change), we will have a function that saves this as a tensor with a NGFF Scale transformation, and then upon loading the scale transformation will be converted to xarray coordinates.

In doing so it will be crucial to have uniformly spaced coordinates, but it doesn’t change much if we use centered vs non-centered coordinates. So we kindly ask you for feedback on this. In the end the argumentations above were more like suggestions for us to use centered coordinates, but we could also go for non-centered ones if many users find this more natural.

The feedbacks we got so far is that centered coordinates are fine to use, but that we need to improve the documentation making it clear that we use them, and showing examples (isel, sel) on how to deal with those.

3

Unrelated comments.

  1. I think that zarr_store: tifffile.ZarrFileSequenceStore = is not needed as the return of tifffile.imread() should already be a Zarr store.

  2. I have never used the approach above: when using small images, we just create the xarray object from the in-memory object, and for large images we use imread() from dask_image to construct the xarray object and then we save this to Zarr. I think that in the case of large images your method would not improve performance because the chunking would still be channel-wise and not across the spatial axes. So I think that one would still have to save to Zarr and reload the new chunked object. But please correct me if I am wrong.

4

Hi @LucaMarconato thank you for help!

  1. think that zarr_store: tifffile.ZarrFileSequenceStore = is not needed as the return of tifffile.imread() should already be a Zarr store.

Often (unfortunately) in VS Code the Python LSP Pylance can get confused with functions that can return multiple types, so in order for it to be useful I just add the type anyway. It’s basically just decor.

  1. I have never used the approach above: when using small images, we just create the xarray object from the in-memory object, and for large images we use imread() from dask_image to construct the xarray object and then we save this to Zarr. I think that in the case of large images your method would not improve performance because the chunking would still be channel-wise and not across the spatial axes. So I think that one would still have to save to Zarr and reload the new chunked object. But please correct me if I am wrong.

I agree that just loading in small images in memory would be easier. The chunking across the spatial axis sounds great with Zarr, and I didn’t even know that there was a dask_image package. I’ll give that a try with some larger data.
I believe you are correct with how the current method would chunk the data from looking at these lines https://github.com/cgohlke/tifffile/blob/e680307a1cbafd77151b37c30c5f4c68c24328f4/tifffile/tifffile.py#L13146-L13168

Thank you!

5

Thank you, this is a very informative overview on how spatialdata saves data!

Centered coordinates for multiscale images make sense, and sounds super interesting to try. I think I’ll try to use this for a larger data test!

In terms of feedback w.r.t. coordinates

  • I agree that using centered coordinates is best
  • Using more of the NGFF builtins like Scale is a good idea

I agree when viewing a pixel as an object with a “mass” (i.e. brightness), utilizing it’s center as the basis of a coordinate system makes complete sense. And a section in the docs explaining this choice would be great!

6

Thanks a lot for your feedback and the comment regarding ZarrFileSequenceStore!

For checking out the multiscale images you can look at the code in the spatialdata-io repository, where we mostly use that (see for instance xenium.py).