Create auxiliary (index) coordinates with a different CRS

[benbovy]

It would be nice to have a simple and common way of assigning auxiliary coordinates, sharing the same dimension(s) but with an alternative CRS, to a dataset / dataarray that already has spatial coordinates with a CRS-aware index. It would be even nicer if the new auxiliary coordinate(s) have an index too. Such feature would fit well in XProj IMO.

---

Example with a raster data cube (rasterix):

>>> ds_raster

<xarray.Dataset>; Size: 120B
Dimensions:      (y: 4, x: 2)
Coordinates:
    spatial_ref  int64 8B 0
  * x            (x) float64 16B ...
  * y            (y) float64 32B ...
Data variables:
    foo          (y, x) float64 64B 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
Indexes:
  ┌ x        RasterIndex
  └ y


>>> ds_raster.proj.create_aux_coords(("lon", "lat"), from=("x", "y"), crs=4326)

<xarray.Dataset>; Size: 168B
Dimensions:      (y: 4, x: 2)
Coordinates:
    spatial_ref  int64 8B 0
  * x            (x) float64 16B ...
  * y            (y) float64 32B ...
  * lon          (y, x) float64 16B ...
  * lat          (y, x) float64 32B ...
Data variables:
    foo          (y, x) float64 64B 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
Indexes:
  ┌ x        RasterIndex
  └ y
  ┌ lon      RasterIndex
  └ lat

Example with a vector data cube (xvec):

>>> ds_vec

<xarray.Dataset (county: 3085)> Size: 35kB
Coordinates:
  * county   (county) object 25kB POLYGON ((...
Data variables:
    foo      (county) float64 10kB ...
Indexes:
    county   GeometryIndex (crs=EPSG:2263)


>>> ds_vec.proj.create_aux_coords("county2", from="county", crs=4326)

<xarray.Dataset (county: 3085)> Size: 60kB
Coordinates:
  * county   (county) object 25kB POLYGON ((...
  * county2  (county) object 25kB POLYGON ((...
Data variables:
    foo      (county) float64 10kB ...
Indexes:
    county   GeometryIndex (crs=EPSG:2263)
    county2  GeometryIndex (crs=EPSG:4326)

---

For the implementation we could rely on xproj.ProjIndexMixin._proj_to_crs(), which returns a new Index instance from an existing CRS-aware index.

The API in the examples above is just a suggestion. It would be nice to also create a new spatial reference, scalar coordinate (with a CRSIndex) alongside the new auxiliary coordinates, perhaps with something like this?

ds.proj.create_aux_coords(..., spatial_ref="spatial_ref2")

@dcherian @scottyhq @martinfleis - Thoughts?

[dcherian]

IIUC The broader use case here is to allow indexing with coordinates from a different CRS.

Perhaps a more general solution is to have xproj provide a wrapper around pyproj transformers, and any CRS-aware index can call that to transform the inputs as needed.

[benbovy]

Indexing with coordinates from a different CRS is indeed one use case, which could also be supported by alternative approaches like the wrapper you suggest or via a custom option passed to .sel() (pydata/xarray#7099).

Besides that, I think that other potential use cases would benefit from the feature proposed here, e.g., align objects having different CRS but also keep their original CRS, visualization, save a dataset with various auxiliary coords / CRSs...