saving and reading vector data cubes

Start by reproducing the reported write_stars failure with the example vector data cube, then read the CF geometry conventions and compare the shapely_to_cf and cfgeom approaches mentioned. The issue does not define a chosen format or implementation entry point; done would require an agreed interoperable storage approach and working round-trip wrappers for xvec-backed xarray objects.

At the moment, we have a decent set of functionality to work with in-memory vector data cubes. What we do not have is a way of saving them to a disk and reading back again, apart from using pickle which is not very interoperable.

We should find a way of saving to a file format that allows us interoperability with @edzer's stars in RSpatial. @edzer, you were mentioning that NetCDF is able to handle VDC. Does stars implements that? Can we create a data cube indexed by geometry and save it via GDAL to NetCDF? I tried only very briefly to play with the write_stars function to no avail[^1].

The spec for NetCDF way of storing geometries is here but as someone who never worked with the file format, it is not super clear to me how it all works together.

Over in Python, there is some work @dcherian posted on conversion of shapely geometries to CF geometries and back (https://cf-xarray.readthedocs.io/en/latest/generated/cf_xarray.shapely_to_cf.html) [^2]. There is also https://github.com/twhiteaker/cfgeom which may potentially be used (although it is only a reference implementation that is not maintained).

So it seems that NetCDF should be able to store vector data cubes and we may just need to figure out some wrappers for a convenient IO from xvec-backed xarray objects. But I'd welcome some help or at least a guidance from someone more familiar with CF conventions.

The other option seems to be Zarr as discussed here but that is at this stage only an idea and (Geo)Zarr spec is not ready as of now.

The last option is to convert the cube to a long-form dataframe and save it as a GeoParquet but that kind of breaks the point of having a cube in the first place.

[^1]: Getting Error in !all.equal(match(xydims, names(d)), 1:2) : invalid argument type when trying to use write_stars(st, 'test.nc') on the example cube from https://r-spatial.github.io/stars/articles/stars1.html#vector-data-cube-example.
[^2]: only points are implemented so far but that may change