DataArray Inheritance from Parent Dataset for Zarr Conventions

[emmanuelmathot]

Problem

Currently, when a Dataset has Zarr conventions declared at the dataset level, individual DataArrays within that Dataset cannot inherit the conventions. This leads to inconsistent behavior:

# Dataset has conventions at group level
ds = xr.Dataset({"var1": (["y", "x"], data)})
ds.attrs = {
    "zarr_conventions": [{"name": "proj:", ...}],
    "proj:code": "EPSG:4326"
}

print(ds.rio.crs)        # ✅ Works: reads from dataset attrs 
print(ds["var1"].rio.crs) # ❌ Returns None: DataArray doesn't inherit

According to Zarr conventions, DataArrays should inherit group-level metadata when they don't have their own conventions declared.

Root Cause

• xarray doesn't maintain parent-child relationships between Dataset and DataArray
• In Implement Zarr spatial and proj conventions support #883 zarr.read_crs(obj) only examines obj.attrs, not potential parent Dataset attributes
• XRasterBase.crs property has no context about inheritance hierarchy

Proposed Solutions

Option 1: Explicit Parent Registration

Add explicit inheritance support to XRasterBase:

class XRasterBase:
    def set_parent_dataset(self, parent_dataset: xarray.Dataset) -> 'XRasterBase':
        """Enable inheritance from parent Dataset"""
        self._parent_dataset = parent_dataset
        return self
        
    def _read_with_inheritance(self, read_func):
        """Try object-level first, then parent inheritance"""
        result = read_func(self._obj)
        if (result is None and 
            isinstance(self._obj, xarray.DataArray) and 
            self._parent_dataset is not None):
            result = read_func(self._parent_dataset)
        return result

# Usage:
da = ds["var1"]
da.rio.set_parent_dataset(ds)  # Enable inheritance
print(da.rio.crs)  # Now inherits EPSG:4326

Pros: Explicit, backward compatible, extensible to all metadata types
Cons: Requires manual setup

Option 2: Enhanced Dataset Accessor

Automatically set up inheritance when accessing DataArrays from Dataset:

class RasterDataset(XRasterBase):
    def __getitem__(self, key):
        """Auto-enable inheritance for DataArray access"""
        result = self._obj[key]
        if isinstance(result, xarray.DataArray):
            result.rio.set_parent_dataset(self._obj)
        return result
        
    def get_array_with_inheritance(self, var_name: str) -> 'RasterArray':
        """Convenience method for inheritance-enabled DataArray access"""
        return self._obj[var_name].rio.set_parent_dataset(self._obj)

Pros: More automatic, better UX
Cons: Requires careful integration with xarray accessors

Implementation Considerations

1. Inheritance Priority: Array-level attributes take precedence over dataset-level
2. Convention Consistency: Same inheritance behavior for both CF and Zarr conventions
3. Metadata Scope: Extend to CRS, transform, and spatial dimensions
4. Performance: Minimal overhead for non-inheritance cases

Breaking Changes

None - inheritance would be opt-in and backward compatible.

Alternative Approaches Considered

• Stack inspection for parent detection
• Global context managers

[snowman2]

I don't think there is a good solution that will be stable. The properties on accessors are very lossy. Once you do an xarray operation, the connection is lost.

It is unfortunate that zarr chose to use attrs. Those are very lossy and not reliable in the xarray world. The CF convention uses a coordinate as it is shared between Dataset and DataArray and is not as easily lost.

I suggest using a coordinate to store the Zarr information to enable persisting the data between Dataset and DataArray.

See:

• Store band specific attributes in coordinates #81
• ENH: Add ability to load bands as data_vars #296
• Copied DataArray loses its set CRS #532
• CRS set on Dataset does not persist for DataArrays obtained by selecting data variable  #356
• Expose sample on the rioxarray accessor #800

[emmanuelmathot]

You're right that attributes are lossy in xarray, it's a well-known pain point. The zarr geospatial conventions have been months (years including geozarr) in the making with extensive community discussions. Switching to coordinates would break that consensus and interoperability.

Given that context, Could we move forward with Option 1? It would:

• Respect the zarr convention design
• Work within xarray's current limitations
• Be backward compatible and opt-in

We can implement it as a practical solution while continuing to advocate for better attribute support in xarray itself.

[snowman2]

Option 1 is worse than attrs. See #532, #356, #252. The set_ methods are highly discouraged because they are so lossy.

[emmanuelmathot]

I appreciate your clarification regarding the inheritance issue in xarray. It seems that, despite our hopes, it may not be applicable in rioxarray.

[benbovy]

Great to see the ongoing work on Zarr/CF multi-convention support! I'm very supportive of the Zarr geospatial conventions.

Just wanted to chime in briefly.

Switching to coordinates would break that consensus and interoperability.

Would it be acceptable to have a Dataset or DataArray that has both Zarr conventions declared and some extra (index) coordinates? (e.g., #855 that I’m guilty of having left it stale).

It is not a perfect solution since there would be duplicated information and perhaps a bit more to check/maintain in the convention classes. But I don’t think it would break consensus and interoperability: couldn’t one “just” ignore or delete those extra coordinates and continue working with an object that has valid Zarr conventions? Maybe deleting those coordinates could be an option added to the rioxarray methods for writing the CRS or transform?

Xarray coordinates do not only allow better persistence of (meta)data across Xarray operations, they also play a key role in the Xarray user experience since lots of Xarray API & features rely on it (e.g., label indexing, alignment, plotting, etc.). For example, it is very handy to have x/y spatial coordinates attached to a Dataset or DataArray, e.g., for using like ds.sel(x=…, y=…), even if those coordinates are not explicitly part of the Zarr spatial conventions. Scalar coordinates could also be useful when they are baked by a custom index like xproj.CRSIndex for CRS-aware auto-alignment of Xarray objects, even if the CRS is also encoded in the Dataset/DataArray as attributes following the Zarr proj convention.

[RichardScottOZ]

Where would that get looked at, in practice...to and from zarr?

[benbovy]

Where would that get looked at, in practice...to and from zarr?

The ZarrConvention class would still be used to read/write georeferencing information to/from zarr.

Actually, thinking more about it: in addition to the ZarrConvention and CFConvention classes there could be another XarrayIndexesConvention class, which instead of relying on attributes would rely on Xarray coordinates baked by the following indexes:

• rasterix.RasterIndex: wraps coordinate (affine) transformation
• xproj.CRSIndex: wraps CRS

(BTW note the perfect alignment in scope of those indexes with the Zarr spatial & proj conventions)

Such index-based convention has the great advantage of being optimized for Xarray UX:

• It provides a robust & flexible way to propagate georeferencing info across Xarray operations, as opposed to lossy attributes and accessor properties
• It allows operating, e.g., on the affine transform directly via Xarray API and x/y coordinates or dimensions (see for example https://xarray-indexes.readthedocs.io/earth/raster.html#indexing)

After a quick inspection, it looks like the convention protocol added in Rioxarray would already work well for implementing such XarrayIndexesConvention class. Switching from one convention to another seem also rather seamless? This is great if it enables the workflow below:

ds_raster_zarr = xr.open_dataset("store.zarr")

# auto-detect Zarr conventions and convert to Xarray indexes
ds_raster_idx = (
    ds_raster_zarr
    .rio.write_transform(convention=XarrayIndexesConvention)
    .rio.write_crs(convention=XarrayIndexesConvention)
)

# do fancy stuff with Xarray geospatial indexes
ds_result = (
    ds_raster_idx
    .isel(x=..., y=...)
    ...
)

# convert the result to zarr conventions
ds_result_zarr = (
    ds_result
    .rio.write_transform(convention=ZarrConvention)
    .rio.write_crs(convention=ZarrConvention)
)

[emmanuelmathot]

Thank you @benbovy for the input. Would XarrayIndexesConvention inherit from ZarrConventions and extend it? So far, the conventions mechanism allows for 1 convention at a time.

Second, I wonder about the performance at opening of the Xarray indices. One of the main driver to have an alternative to CF convention was to avoid opening other arrays to just get attributes from a 0D array. Even if lossy, attributes have the advantage to be available in the consolidated metadata.

[benbovy]

No XarrayIndexesConvention would be independent, and I think it would work well with the current mechanism of one convention at a time.

The basic idea of XarrayIndexesConvention would be to either look or create RasterIndex or CRSIndex instances like this:

# this is a bit simplified, although quite close to the actual implementation!

class XarrayIndexesConvention:

    @classmethod
    def read_crs(cls, obj):
        for index in obj.xindexes.values():
            if isinstance(index, xproj.CRSIndex):
                return index.crs
        return None

    @classmethod
    def read_transform(cls, obj):
        for index in obj.xindexes.values():
            if isinstance(index, rasterix.RasterIndex):
                return index.transform
        return None

    @classmethod
    def read_spatial_dimensions(cls, obj):
        for index in obj.xindexes.values():
            if isinstance(index, rasterix.RasterIndex):
                return index.xy_dims
        return None

    @classmethod
    def write_crs(cls, obj, crs, **kwargs):
        return obj.set_xindex("crs", xproj.CRSIndex, crs=crs)

    @classmethod
    def write_transform(cls, obj, *, transform, **kwargs):
        index = rasterix.RasterIndex.from_transform(transform, **kwargs)
        # Xarray lazy coordinate variables wrapping the affine transform
        spatial_coords = xr.Coordinates.from_xindex(index)
        return obj.assign_coords(spatial_coords)

XarrayIndexConvention is thus neither opening any array nor looking at any attribute. For Rioxarray to detect the georeferencing info using this index-based convention, it needs either that:

• RasterIndex and CRSIndex instances have been created / assigned either manually (e.g., creating a Dataset or DataArray from scratch) or using a 3rd-party Xarray extension
• the georeferencing info has been detected using ZarrConvention or CFConvention (e.g., on a newly opened Dataset) and then explicitly written using XarrayIndexConvention (code example in DataArray Inheritance from Parent Dataset for Zarr Conventions #888 (comment))

Regarding performance, create/assign CRSIndex and RasterIndex objects basically consists in wrapping CRS and Affine Python objects then create lazy or scalar Xarray variables, so it is almost free I would say.

---

Sorry, this deviates from the original topic of this issue. The index-based convention works well with DataArray-Dataset inheritance, though, so assuming we can - and we want to - easily switch between conventions (one at a time) in a processing pipeline it could provide a solution to this issue.

[snowman2]

xarray indexes aren't a convention for storing metadata on a xarray object that can be persisted on disk. They indexes are a way of representing convention metadata on the xarray object.

Instead, to_netcdf and to_zarr should be able to interact with the indexes to convert them to the correct on disk representation.

[benbovy]

xarray indexes aren't a convention for storing metadata on a xarray object that can be persisted on disk.

I still see the model with RasterIndex + CRSIndex (and their associated coordinates) as a convention.

Indeed it is not a convention for writing metadata on "disk" like CF and Zarr conventions. But it is still a conventional representation of Xarray objects that contain georeferencing information, which can be used across any Xarray library that supports it (e.g., gdx, xpublish-tiles... maybe other libraries in the future like odc-geo?), for better interoperability within the Xarray geospatial ecosystem.

Is the Rioxarray convention mechanism strictly limited to reading/writing metadata that is meant to persist on disk? If yes, it would be good to clarify that as I didn't see there were such restriction.

What would be the reasons to prevent generalizing Rioxarray's convention concept to simply "where to find the georeferencing info" and "where to put the georeferencing info" in an Xarray object?

Rioxarray's CFConvention and ZarrConvention are parsing (creating) Xarray attributes and/or coordinate variables... Is there any reason why another convention wouldn't be allowed to parse (create) Xarray indexed coordinates?

If an XarrayIndexesConvention is out of scope or doesn't make sense, how could (should) we tell Rioxarray to get the CRS and/or transform from an Xarray object that we created from scratch with a CRSIndex and/or RasterIndex?

Instead, to_netcdf and to_zarr should be able to interact with the indexes to convert them to the correct on disk representation.

Those functions are in general agnostic of the Xarray indexes. Instead, writing the georeferencing info to variables and attributes using either the CF or Zarr conventions is a prerequisite.

[snowman2]

The CF & Zarr conventions are able to be understood by software outside of xarray. Being able to persist on disk enables that.

The xarray indexes are an abstraction layer constructed using information from either the CF or Zarr conventions.

[snowman2]

Here is an example of what it could look like where you read the convention metadata from the xarray object to build the index:

ds_raster_zarr = xr.open_dataset("store.zarr")
# build xarray indexes from Zarr/CF conventions
ds_raster_idx = (
    ds_raster_zarr
    .set_xindex("crs", xproj.CRSIndex, crs=ds_raster_zarr.rio.crs)
    .assign_coords(xarray.Coordinates.from_xindex(rasterix.RasterIndex.from_transform(ds_raster_zarr.rio.transform()))
)

# do fancy stuff with Xarray geospatial indexes
ds_result = (
    ds_raster_idx
    .isel(x=..., y=...)
    ...
)

[benbovy]

Yes sure one could build the indexes manually.

Then we have this kind of issue:

# returns the wrong transform read from Zarr attr (not updated after slicing)
ds_result.rio.transform 

#  returns the right transform (updated after slicing)
ds_result.xindexes["x"].transform

Also by adding indexes we add their corresponding (spatial + "grid mapping") coordinates, and we're back to

Switching to coordinates would break that consensus and interoperability.

Of course we could let users manually deal with all those inconsistencies... But I thought it would be nice for Rioxarray to provide something to prevent this kind of situation...

It turns out that the convention protocol newly added in rioxarray already provides a nice framework for opt-in integration of rasterix and xproj with Rioxarray, as shown in #888 (comment).

I honestly don't fully get your point in #888 (comment), but I'd understand that you want to restrict Rioxarray conventions to metadata that can be persisted on disk: it is easy for anyone to relate the ZarrConvention and CFConvention classes here to their "official" declaration, which isn't the case for the index-based rioxarray convention proposal.

(sorry again for hijacking this issue, we can move this discussion into a new issue, if anyone would like to continue it).