Adding a new array to a dataset

[benjeffery]

It's a common use case to want to add an array to a dataset - in my case adding some variant masks based on duplicated sites and quality metrics. However, I don't want to have to rewrite the whole dataset every time I add an array.
pydata/xarray#6700 suggests this is not possible, and https://stackoverflow.com/questions/58042559/adding-new-xarray-dataarray-to-an-existing-zarr-store-without-re-writing-the-who suggests some clever use of the 'a' mode in to_zarr with only the new arrays in a dataset, and specified regions:
ds2.to_zarr("test.zarr/", mode="a", region={"x": slice(0, ds2.x.size), "z": slice(0, ds2.z.size)})

The other option is to write to the zarr files directly, but that seems a bit hacky.

[tomwhite]

You should be able to do it with the 'a' mode (no need to specify regions): https://github.com/pystatgen/sgkit/issues/347#issuecomment-730400413

[jeromekelleher]

Issue 1000!

[benjeffery]

Dammit, missed the chance for a party thread. tskit-dev/tskit#1000

[jeromekelleher]

2**10 is still available! 1000 is just a reminder that the big one is coming up.

[benjeffery]

I ended up with this - which seems to work but was non-obvious initially:

    ds.update({"variant_mask": xr.DataArray(data=selected_variants,dims=["variants"],name="variant_mask")})
    sgkit.save_dataset(ds.drop_vars(set(ds.data_vars) - {"variant_mask"}), zarr_path, mode="a")

This is a very common pattern - I wonder if we should have a function such as:
sgkit.update_dataset({"name":array}) which overwrites by default?

[tomwhite]

+1 to update_dataset. Worth thinking about the API a bit - it would be nice to be able to give it a dataset that has been returned from an sgkit method with merge=False.

[benjeffery]

Yes, was thinking the same thing - the usual workflow in sgkit is to add some arrays using sgkit methods, so merge=False could be the default, however I think that is surprising and I would expect merge=True to be the default for an updating method.
So my workflow would be:

ds = ds.update({....})
sg.update_dataset(ds, path, merge=False)

Which seems more obvious than what I have now.

[tomwhite]

I meant merge=False on the sgkit method, not on update_dataset. E.g.

ds = sg.divergence(ds, merge=False) # ds now only has stat_divergence variable
sg.update_dataset(ds, path)