This looks like a nice use case for the forthcoming Xarray's custom index feature.

How I see CRS/WCS-aware Xarray datasets with custom indexes:

• A set of coordinate(s) and their attributes hold data or metadata relevant for public use and that could be easily (de)serialized

• A custom index (CRSIndex or WCSIndex) provides CRS/WCS-aware implementations of common Xarray operations such as alignment (merge/concat) and data selection (sel), via Xarray.Index's equals, union, intersection and query methods added in #5102 and #5322 (not yet ready for use outside of Xarray). Such custom index may also be used to hold some data that is tricky to propagate by other means, e.g., some internal information like "functional" coordinate parameters or a crs object. Xarray indexes should definitely provide more flexibility than coordinate data or attributes or accessor attributes for propagating this kind of information.

• Xarray accessors may be used to extend Dataset/DataArray public API. They could use the information stored in the CRSIndex/WCSIndex, e.g., add a crs read-only property that returns the crs object stored in CRSIndex, or add some some extract_crs_parameters method to extract the parameters and store them in Dataset/DataArray attributes similarly to what @djhoese suggests in his comment above.

For this use case a possible workflow would then be something like this:

```python

create or open an Xarray dataset with x, y, z "pixel" (possibly lazy) coordinates

and set a WCS index

dataset = ( xr.Dataset(...) .set_index(['x', 'y', 'z'], WCSIndex, wcs_params={...}) )

select data using pixel coordinates

dataset.sel(x=..., y=..., z=...)

select data using world coordinates (via the "astro" accessor,

which may access methods/attributes of the WCS index)

dataset.astro.sel_world(x=..., y=..., z=...)

return a new dataset where the x,y,z "pixel" coordnates are replaced by the "world" coordinates

(again using the WCS index, and propagating it to the returned dataset)

world_dataset = dataset.astro.pixel_to_world(['x', 'y', 'z'])

select data using world coordinates

world_dataset.sel(x=..., y=..., z=...)

select data using pixel coordinates (via the "astro" accessor)

world_dataset.astro.sel_pixel(x=..., y=..., z=...)

this could be reverted

pixel_dataset = world_dataset.astro.world_to_pixel(['x', 'y', 'z']) assert pixel_dataset.identical(dataset)

depending on the implementation in WCSIndex, would either raise an error

or implicitly convert to either pixel or world coordinates

xr.merge([world_dataset, another_pixel_dataset]) ```