Utilities for writing/updating SpatialData Zarr stores

We can implement methods like create_zarr and create_zarr_group until they are available in Rarr
Huber-group-EMBL/Rarr#18

@HelenaLC

@HelenaLC

Is checking for dask index necessary here ? is it failing for any of the Spatialdata.data examples without ?

https://github.com/HelenaLC/SpatialData/blob/1c143790a1adcdf8b9013009df5fdadb844aacc1/R/PointFrame.R#L108-L123

Did also a quick check on indexing, the regular indexing seem also be doing stuff lazily.

> library(SpatialData.data)
> x <- SpatialData.data::MouseBrainMERFISH()
> object.size(x@points$single_molecule[1:100000,]@data)
488 bytes
> object.size(x@points$single_molecule@data[1:100000,])
488 bytes

@kollo97

@kollo97 considerably improved the writing support. PR incoming ...

I'm not sure as I didn't really understand what this is/if it's needed. At some point I was even hiding it from the user as I find it confusing. Shall we do that? Ie, not protect it from dropping and also hide it from the user (if it exists)?

ok, the question was old and __null_dask_index__ is nowhere to be found anymore, so everything works as intended no problem ?

Well, it is to be found when you read a point element. Just not sure if maybe we should "hide" it from the user, because it's not helping anyone...

Blobs dataset:

> point(sd)
class: PointFrame
count: 200           
data(4): genes instance_id __null_dask_index__ geometry

Also related to versioning of sd zarr stores. We should not simply call versions of sd v2 and v3 since SpatialData versions are defined by combinations of element versions, see here:
https://github.com/scverse/spatialdata/blob/5ead6cc3d25226b46561baa64d02b92ddeadbc10/src/spatialdata/_io/format.py#L309-L346

ContainerFormatValidElements = {
    SpatialDataContainerFormatV01().__str__(): [
        RasterFormatV01().__str__(),
        RasterFormatV02().__str__(),
        ShapesFormatV01().__str__(),
        ShapesFormatV02().__str__(),
        PointsFormatV01().__str__(),
        TablesFormatV01().__str__(),
    ],
    SpatialDataContainerFormatV02().__str__(): [
        RasterFormatV03().__str__(),
        ShapesFormatV03().__str__(),
        PointsFormatV02().__str__(),
        TablesFormatV02().__str__(),
    ],
}

Here, the version argument should be either 0.1 or 0.2 which create v2 and v3 zarr stores, respectively.

>>> import spatialdata as sd
>>> from spatialdata._io.format import (
>>>   SpatialDataContainerFormatV01,
>>>   SpatialDataContainerFormatV02
>>> )
>>> sdata = sd.datasets.blobs()
>>> sdata.write("temp.zarr", sdata_formats = SpatialDataContainerFormatV01())
(dummy_sd_env) LANDTHALER57:~ amanuky$ tail temp.zarr/labels/blobs_labels/.zattrs 
      ]
    }
  ],
  "image-label": {
    "version": "0.4"
  },
  "spatialdata_attrs": {
    "version": "0.2"
  }
}
(dummy_sd_env) LANDTHALER57:~ amanuky$ tail temp.zarr/points/blobs_points/.zattrs 
        ]
      }
    }
  ],
  "spatialdata_attrs": {
    "instance_key": "instance_id",
    "feature_key": "genes",
    "version": "0.1"
  }
}

By golly. Que confusion...

@HelenaLC

Hi @HelenaLC, @Artur-man,

Following on from #115, we'd like to help get table writing into spatialdataR, and we have a proposal for how.

Summary

The problem. spicyR now reads SpatialData stores a block of rows at a time, so we needed to write stores of hundreds of millions of cells from R. Two things stop that today:

spatialdataR can't yet write tables in Zarr v3. The writeTable() in Write Support #247 called anndataR::write_zarr() and refused v3.
anndataR (1.3.2) does write Zarr v3, but it puts every array in a single, uncompressed chunk. That ~~fails on large tables~~ needs about six times the table's memory to write (correction), and stores written that way can't be read in blocks.

The fix is chunked writing: cut each array along the observations and write it one chunk at a time. We've built and tested this in spicyR (details below). Writing 100 million cells takes about 2 minutes, with memory near the size of the input.

Where it could live. Two options:

A. Upstream first (our preference). Chunked writing goes into anndataR, where it helps every R package that writes AnnData. spatialdataR's writeTable() then stays a thin wrapper: anndataR::write_zarr() plus the SpatialData table attributes. You don't take on any AnnData encoding code. The cost is waiting for anndataR to accept the change.
B. spatialdataR first. We contribute a self-contained chunked table writer to spatialdataR now (about 300 lines, written to match the writeTable(x, name, path, ...) interface from Write Support #247), and replace it with a call to anndataR once anndataR chunks its writes. You get Zarr v3 tables straight away, but spatialdataR carries its own copy of the AnnData encoding in the meantime.

Either way, we'll do the work: the anndataR PR for A, or the spatialdataR PR for B. We'll also coordinate with whoever is splitting up #247.

Two smaller things, under either option:

Rarr. Writing chunk by chunk through Rarr is slower than it needs to be, and Rarr can't write the string type anndata uses. We'll propose both changes to Rarr (details below).
Later, SpatialExperiment to SpatialData. A conversion that builds the table, the cells and the images from a SpatialExperiment and writes them with your writers. It would complement Converting images, labels and shapes to SpatialExperiment / SpatialFeatureExperiment #251, which covers the opposite direction. We'd only start on this once the table writer exists, and we'd like your view on the layout first (questions at the end).

What we'd like from you:

A or B?
Is anyone already working on Zarr v3 tables in the PRs replacing Write Support #247?
For the SpatialExperiment conversion: the layout questions at the end.

Evidence

Single-chunk writing doesn't scale. anndataR's zarr_write_compressed() calls Rarr::write_zarr_array() with chunk_dim = dims, so each array is a single chunk, and compression defaults to "none". ~~When spicyR wrote its tables this way, it failed at 100 million cells while writing the obs index.~~ (That was with older anndataR and Rarr; current anndataR writes 100 million cells, with a 51 GB peak: see the correction.) Reading any slice of a column means decompressing the whole column. (The writeImage() in #247 also writes each scale as one chunk, chunk_dim = dim(arr); the same reasoning applies to large images.)

Chunked writing does. We cut every array along its first axis: 2^20 rows for a column, and whole rows for a matrix (about two million values per chunk). Each chunk is compressed with zstd and written before the next is made. On 100 million cells (two categorical columns and two coordinates; 2.4 GB as a data.frame):

	time	peak memory (whole process)	store size
chunks encoded by our own code	85 s	3.4 GB	1.6 GB
the same chunks written through Rarr	121 s	3.4 GB	1.6 GB
through Rarr, plus a points element per image	197 s	3.4 GB	3.4 GB

Memory is the input plus about a gigabyte, whatever the number of cells. With stores chunked this way, spicyR analyses one billion cells in 500 images on one core in under 2 GB.

The cost of going through Rarr. create_empty_zarr_array() followed by update_zarr_array() per chunk produces exactly the same compressed chunks as our own encoder. It's 43% slower, though: 36 extra seconds at 100 million cells, or roughly 6 minutes at a billion, for a write done once per dataset. The reason is that update_zarr_array() works out which chunk each element of index falls in (split() and as.factor() in R), even when the index covers whole chunks. Profiling 16 chunks of 2^20 values:

For doubles, that step is about 30% of the time (2.7 s against 2.2 s for our encoder).
For small integer columns such as category codes, where compression is cheap, it's 1.7 s against 0.4 s.

We kept Rarr anyway, so that the code stays close to yours. A fast path in update_zarr_array() for indices that cover whole chunks would remove the overhead, and we'll propose it to Rarr.

Strings. Rarr writes character arrays as fixed-length null_terminated_bytes, while anndata and Python's zarr use vlen-utf8. So we encode string arrays (indices, string columns, categories) ourselves. This is the one part not written through Rarr, and we'll also propose vlen-utf8 to Rarr.

Our working version

R/writeSpatialData.R on spicyR devel, called by spicyR::writeCellStore(). It's structured like #247 (.writeTable(), .writeImage(), .writePoints()), so the table part could move with little change. From a SingleCellExperiment or SpatialExperiment it writes the following.

The table (anndata encoding 0.1.0):

X: the first assay. A dgCMatrix is written as csr_matrix straight from its slots; any other matrix, including a DelayedMatrix, is written dense, a block of cells at a time.
layers: the other assays.
obs: the colData, with categoricals and nullable columns. The image column is the region key, and instance_id numbers the cells.
var: the rowData.
obsm: spatialCoords() as "spatial", and the reducedDims.

The cells as points. One element per image, named after the image, with x, y and an int64 instance_id that links each point to its row of the table; the table annotates these elements. Each chunk of rows is written to parquet with duckdb, and one duckdb pass then splits them into a single file per image, so memory stays at one chunk however the cells are ordered.

Images. Each image in imgData() becomes an image element named <sample_id>_<image_id>, at scale 1 / scaleFactor to its sample's coordinate system. The points of a sample's cells sit in the same coordinate system.

Versions. Zarr v3; container 0.2; points 0.2; images OME 0.5-dev-spatialdata with spatialdata_attrs 0.3, as spatialdata 0.7.3 writes them.

Limitations

The image writer handles RGB images only (it works from imgRaster()), at one scale, and reads each image whole. Only the cell data is chunked from end to end.
The first assay comes back as X, so its name is lost.
metadata(), altExps() and the column and row pairings aren't written.

Tests (test-writeSpatialData.R) write a SpatialExperiment with:

sparse and dense assays, rowData and a PCA;
cell names and missing values;
six images across two samples, each sample with an image.

They use small chunks so every array has several, and check that readSpatialData() and Python's spatialdata 0.7.3 return every value, point and pixel exactly. Python runs with UserWarning and anndata's OldFormatWarning turned into errors. Requirements: Rarr ≥ 2.1.46, plus duckdb and DBI for the points.

Questions for the SpatialExperiment conversion

These are the choices we made for spicyR. We'd follow your conventions instead:

The cells: points per image (what we do), shapes such as circles, or no element?
Coordinate systems: one per sample, or one per image?
Element names: <image> for points and <sample_id>_<image_id> for images?
Images with more than three channels (IMC, CODEX): how would you like those written?

A separate bug

readPoint() reads only points.parquet/part.0.parquet (R/read.R:80). Python's spatialdata writes one part per dask partition, so points built from a partitioned dask frame lose rows on reading. A 5-million-point element written from a dask frame with 4 partitions reads back as 1,250,000 points, with no warning. We'll open a separate issue with that example.

Follow-ups: #287 (readPoint() reads only the first parquet part), Huber-group-EMBL/Rarr#260 (fast path in update_zarr_array(), with a tested patch) and Huber-group-EMBL/Rarr#261 (data_type = "|O" writes unreadable string arrays, with a tested fix).

@ellispatrick

Dear @ellispatrick,

Thanks so much for giving a heads up. We are normally working towards writeTable and all other write utilities that you mentioned but I was not able to get to it yet due to other projects.

In fact, if you can send the PRs you mentioned, and we can finally start. Thus, I will review those PRs that would be super cool.

One thing to be careful of is the management of SpatialData versioning since a Zarr v3 table should not be written to a SpatialData store of version 0.1 (associated with Zarr v2). I was getting prepared to make an update to manage versions in spatialdataR so that the writing happens only for the SD version 0.2 (Zarr v3). I will send a PR for that at some point and merge will yours later to solidify the version controlling while writing the table.

Looking forward to the PR, and thanks for your contributions in advance!

Hmm, looking at this now, some of these might already being handled by anndataR::write_zarr in BioC devel (Zarr v3 writing may not be in Bioc release now, I am not quite sure, eventhough I wrote the PR hehe)
https://github.com/SydneyBioX/spicyR/blob/fc1719c1c463f96090e82a905a67dfb5f61fa7c3/R/writeSpatialData.R#L251-L266

You guys wanna have a chat ?

@HelenaLC

Answers to the questions:

I think we should go with option A first, I can also check once the anndataR PR is set. Do you guys have a working version of your writeTable somewhere, we can perhaps start B in the same time if I can take a look.
It would be me but would love to collaborate on that if you guys are up for it :) as I mentioned in the previous comment, we need a bit auxiliary tools for version control. We can do it while working on the PRs.
SpatialExperiment conversion I do not know @HelenaLC ? You have more opinions on this.

@Artur-man

Thanks @Artur-man, that's great. Our working version is .writeTable() in R/writeSpatialData.R (spicyR devel); writeCellStore() calls it. It only ever creates new stores, always SpatialData 0.2 with Zarr v3, so it never writes into an existing 0.1 store. The version check you describe is what's needed once spatialdataR writes tables into existing stores, and we'll build on your versioning PR.

We'll start with the anndataR side (A) and, in parallel, open a draft PR here for writeTable() (B), in the shape of #247 and limited to Zarr v3 / SpatialData 0.2. We'll keep it as a draft until your versioning PR is in, and rebase onto it.

A correction to the evidence in our comment above. We re-ran anndataR's write_zarr() from current main (anndataR 1.3.3, Rarr 2.1.46), and it does write 100 million cells. The failure we described came from the older anndataR and Rarr that spicyR used when we first wrote this code, so that sentence is wrong for the current versions; apologies.

What still holds is the cost. Writing cells with their names, a factor column and a numeric column:

cells	time	peak memory	chunks per array	store (uncompressed, the default)
1 million	7 s	0.9 GB	1	23 MB
10 million	73 s	5.9 GB	1	238 MB
50 million	377 s	28 GB	1	1.3 GB
100 million	730 s	51 GB	1	2.5 GB

The table itself takes about 8.4 GB in R at 100 million cells, so writing peaks at about six times its size. Every array is one chunk, so the store can't be read a block of rows at a time. That's the case we'll make in the anndataR issue, with measurements of chunked writing alongside.

If the default compression in writing anndataR tables is also "none", thats a problem. We should default it to zstd like Rarr does:
https://huber-group-embl.github.io/Rarr/reference/write_zarr_array.html

@Artur-man

Thanks @Artur-man. Agreed on zstd: with chunking it costs almost nothing and makes the store about 4.7 times smaller (2.5 GB → 532 MB at 100 million cells). We've raised it in scverse/anndataR#541, together with chunking, and opened scverse/anndataR#543 for the chunking itself. You're right that anndataR::write_zarr() on devel already writes Zarr v3 tables, which is why we'd rather add the chunking there (option A) than keep a writer of our own. Most of the remaining write time turned out to be Rarr's string encoder, so we've proposed a C encoder in Huber-group-EMBL/Rarr#262.

A chat would be great; it would be good to know what your plans are.

A correction to the memory figures in my comment above. The peak memory there was for the whole R process, which includes building the test table, and creating 100 million cell names itself takes a lot of memory. Measuring only what writing adds (R's peak memory during write_zarr(), above what the table already uses; anndataR 1.3.3, Rarr 2.1.46), on the same table of names, a factor column and a numeric column:

cells	table in R	added by writing
1 million	0.1 GB	0.6 GB
10 million	1.0 GB	4.7 GB
50 million	4.9 GB	23.8 GB
100 million	9.8 GB	41.4 GB

So writing needs about four times the table's size on top of it, not six. With chunking (scverse/anndataR#543), writing 100 million cells adds 16.1 GB, and with the C string encoder (Huber-group-EMBL/Rarr#265) as well, 4.4 GB.

A chat would be great; it would be good to know what your plans are.

I will DM you for this on Zulip.

@Artur-man

Thanks @Artur-man, that's great. Our working version is .writeTable() in R/writeSpatialData.R (spicyR devel); writeCellStore() calls it. It only ever creates new stores, always SpatialData 0.2 with Zarr v3, so it never writes into an existing 0.1 store.

#288 is now merged, which can be used to check the format version of the store before writing tables and other. Ultimately (and as you also mentioned) we should not write an SCE as table if the container/store is not version 0.2.

The important point here is that, there is no version SpatialData version 0.2 written anywhere in the store (metadata). We must deduce this by checking all other elements. The current version schema for sd version 0.2 includes version 0.3 of shapes, images, labels and 0.2 of points. We will simply ignore all other combos of versions (hoping that versioning will be simplified in the future).

If it is going to make everybody's job easier, we can first limit writeTable to updating existing tables, thus would require only checking the table version, which is supposed to be 0.2.

Navigation Menu

Repository navigation

Utilities for writing/updating SpatialData Zarr stores #115

1 remaining item

Summary

Evidence

Our working version

Questions for the SpatialExperiment conversion

A separate bug

Metadata

Assignees

Labels

Projects

Milestone

Relationships

Development

Utilities for writing/updating SpatialData Zarr stores #115

Description

Activity

Artur-man commented on May 30, 2025

Artur-man commented on Apr 24, 2026

HelenaLC commented on May 3, 2026

Artur-man commented on May 3, 2026

HelenaLC commented on May 3, 2026

1 remaining item

Artur-man commented on May 6, 2026

HelenaLC commented on May 6, 2026

ellispatrick commented on Oct 22, 2026

Summary

Evidence

Our working version

Questions for the SpatialExperiment conversion

A separate bug

Artur-man commented on Oct 22, 2026

Artur-man commented on Oct 22, 2026

ellispatrick commented on Oct 22, 2026

ellispatrick commented on Oct 22, 2026

Artur-man commented on Oct 22, 2026

ellispatrick commented on Oct 22, 2026

ellispatrick commented on Oct 22, 2026

Artur-man commented on Oct 23, 2026

Artur-man commented on Oct 29, 2026

Metadata

Metadata

Assignees

Labels

Projects

Milestone

Relationships

Development

Issue actions