issue_comments

I've proposed adding CRSIndex to rioxarray to experiment with propagating CRS info in existing workflows: https://github.com/corteva/rioxarray/issues/588

could be moved into some PandasMetaIndex helper class that would encapsulate one or more PandasIndex

Agreed.

reuse the CRS-related logic with other kinds of index structures (like kd-trees). I've been thinking a bit about the general issue of flexible geospatial xarray indexes but I'm not sure yet how best it could be solved.

Yeah I don't think the PandasMetaIndex is a good pathway for CRSIndex, since we'd want other underlying tree-like structures.

be supported with https://github.com/pydata/xarray/pull/6800

Hadn't seen that. That would be great!

the variables argument of Index.create_variables should be optional or not.

I actually didn't understand what I was supposed to do with variables :) . It would be nice to document in PandasIndex

There's some discussion in https://github.com/pydata/xarray/issues/4366 about adding a new .drop_indexes() (or drop_xindexes()?) method.

My takeaway was that on the API side, we should prioritize adding set_xindex and drop_indexes. These are critical for experiments.

Re: reduction to scalar, I agree it seems tricky.

That's great @dcherian! Some comments (notably regarding your notes in your linked notebook):

A lot of boilerplate code in your CRSIndex class (and in RasterIndex) could be moved into some PandasMetaIndex helper class that would encapsulate one or more PandasIndex. Then CRSIndex could subclass it and just implement the CRS-related logic. This would be useful for many other indexes I think.

But as you suggest it, it would be nice if we could also reuse the CRS-related logic with other kinds of index structures (like kd-trees). I've been thinking a bit about the general issue of flexible geospatial xarray indexes but I'm not sure yet how best it could be solved.

Can we support passing a CRS object directly using **kwargs? Not at the moment. This is the data model: Indexes are constructed from a subset of Coordinate variables. all necessary information should be in a coordinate variable

This should be supported with #6800 (I need to re-submit a PR targeting main, see #6849). For the case of CRSIndex, I think it might be useful to support both ways (i.e., either set the CRS via **kwargs or get it from a coordinate variable).

Index.create_variables. set_xindex tries to pass a variables kwarg but other methods don't. Seems like a bug.

Yes this should be clarified, i.e., whether the variables argument of Index.create_variables should be optional or not. I think we can make it required but I need to check.

How do we delete an existing index?

There's some discussion in #4366 about adding a new .drop_indexes() (or drop_xindexes()?) method.

error message with join='exact' is very generic.

That's a tricky one to improve.

GroupBy doesn't propagate index.

Maybe related to #6836 ?

spatial_ref is confusing. A scalar index?

I agree. The index should probably be dropped in that case (i.e., reduction of both the x and y dimensions), leaving the spatial_ref coordinate as a non-indexed coordinate. But this may be tricky to handle in general (i.e., for any index). Ideally this should also be consistent with the behavior of the index for other operations like scalar selection.

Here is an alpha version of a CRSIndex heavily drawing on @benbovy's RasterIndex https://github.com/dcherian/crsindex/blob/main/crsindex.ipynb

As a non-expert, I very arbitrarily chose to propagate CRS info using a spatial_ref scalar variable. When assigned to a dataset we get

so spatial_ref is bolded but is not a dimension, and is associated with an index.

Is there a way in xarray to associate these three arrays in a DataArray so that slicing is handled automatically but also not put the arrays in the coordinates?

Not yet, unfortunately, but this is what https://github.com/pydata/xarray/pull/2302 is trying to solve.

I could always add this logic myself to geoxarray's version of to_netcdf.

I think this would be the preferred approach.

But I guess that is intended behavior and if the crs is a coordinate then joining things from different projections would not be allowed and raise an exception. However that is exactly what satpy wants/needs to handle in some cases (satellite datasets at different resolutions, multiple 'regions' of from the same overall instrument, two channels from the same instrument with slightly shifted geolocation, etc).

I think it would make more sense to think about using multiple xarray.Dataset objects for these use cases, possibly in some sort of hierarchical collection. The notion of xarray.Dataset is pretty closely tied to a single grid. The discussion in https://github.com/pydata/xarray/issues/1092 is definitely worth reading.

Correct me if I'm wrong, but from the xarray side it would already be enough if there is a way in xarray to have a special crs attribute which is preserved by all operations. What the crs really is shouldn't bother xarray, and it can be specific to the downstream library [1]. With a crs and xarray's coordinates the geoloc can be sorted out in almost all cases, right?

The solution currently is to store crs as a scalar coordinate as mentioned above.

What could also be a possibility (without knowing how the internals would look like) is to have a registry of "special attributes" names which would always be preserved by xarray's operations. This registry would live in xarray and can be updated by downstream libraries and/or accessors. (xref: https://github.com/pydata/xarray/issues/1614 ).

[1] : my preference goes for a simple PROJ4 string

This has been discussed from time to time, although not in such an extensive way: thanks!

For me, the long list of discussion points that you mention is a clear argument for a dedicated package which will have to solve all these issues. A good example of difficulty to tackle is the fact that geo libraries often do not use the same standards (see the incompatibility between pyproj and cartopy.crs objects as a striking example).

Regarding the xarray side: I also brought this up for my own geo-lib a while ago, and it seems that the cleanest solution to carry a crs info is to store it as (scalar) coordinate instead of an attribute, which are lost in arithmetic computations (see @shoyer 's comment).

geopandas would be a great template for a geoxarray package ;-)