Merging CRSs

[benbovy]

A useful feature (regardless of the single vs. multi CRS model discussed in #2) would be to merge all the CRSs found in a dataset or dataarray.

Let's take the example below where two datasets differ only by the name of the CRS-indexed coordinate:

>>> ds = xr.Dataset(coords={"lat": [1, 2, 3], "lon": [1, 2, 3]})
>>> ds_wgs84 = ds.proj.assign_crs(crs_wgs84="epsg:4326")
>>> ds_spatial_ref = ds.proj.assign_crs(spatial_ref="epsg:4326")

>>> ds_wgs84
<xarray.Dataset> Size: 56B
Dimensions:    (lat: 3, lon: 3)
Coordinates:
  * lat        (lat) int64 24B 1 2 3
  * lon        (lon) int64 24B 1 2 3
  * crs_wgs84  int64 8B 0
Data variables:
    *empty*
Indexes:
    crs_wgs84  CRSIndex (crs=EPSG:4326)

>>> ds_spatial_ref
<xarray.Dataset> Size: 56B
Dimensions:      (lat: 3, lon: 3)
Coordinates:
  * lat          (lat) int64 24B 1 2 3
  * lon          (lon) int64 24B 1 2 3
  * spatial_ref  int64 8B 0
Data variables:
    *empty*
Indexes:
    spatial_ref  CRSIndex (crs=EPSG:4326)

Merging these two datasets with xarray.merge() will keep both CRS coordinates:

>>> merged = xr.merge([ds_wgs84, ds_spatial_ref])
>>> merged
<xarray.Dataset> Size: 64B
Dimensions:      (lat: 3, lon: 3)
Coordinates:
  * lat          (lat) int64 24B 1 2 3
  * lon          (lon) int64 24B 1 2 3
  * crs_wgs84    int64 8B 0
  * spatial_ref  int64 8B 0
Data variables:
    *empty*
Indexes:
    crs_wgs84    CRSIndex (crs=EPSG:4326)
    spatial_ref  CRSIndex (crs=EPSG:4326)

Getting the unique CRS from the merged dataset won't work despite all CRS are the same:

>>> merged.proj.crs
ValueError: found more than one coordinate with a CRSIndex in Dataset or DataArray

For convenience we could provide a proj.merge_crs() method such that:

>>> merged = merged.proj.merge_crs("spatial_ref")
>>> merged
<xarray.Dataset> Size: 56B
Dimensions:      (lat: 3, lon: 3)
Coordinates:
  * lat          (lat) int64 24B 1 2 3
  * lon          (lon) int64 24B 1 2 3
  * spatial_ref  int64 8B 0
Data variables:
    *empty*
Indexes:
    spatial_ref  CRSIndex (crs=EPSG:4326)
>>> merged.proj.crs
<Geographic 2D CRS: EPSG:4326>
...

We could first restrict .proj.merge_crs("coord_name") to checking that all CRS indexes found in the dataset are equal and if yes merge all CRS coordinates into one (with the name given as argument). Later on we could expand the functionality with optional arguments (e.g., trigger re-projection).