Decode x/y coordinates for compatibility with `rioxarray` and `odc-geo`

[peanutfun]

Is your feature request related to a problem? Please describe.

When using cfgrib to create an xarray Dataset from a GRIB file, the lat/lon coordinates for each x/y point of the data are decoded. However, the x/y coordinates are not decoded. These dimensions of the original dataset are not indexed and have no coordinates. As such, this dataset representation is incompatible with other geospatial data packages like rioxarray and odc-geo, which store explicit CRS information and matching spatial coordinates. They thus allow for resampling or remapping to different CRSs and/or grids.

Additionally, it is cumbersome to assign CRS information because the projString variable is not automatically read.

>>> import xarray as xr
>>> FILE = ".../file.grib"  # An EFAS v4 file from https://ewds.climate.copernicus.eu/datasets/efas-forecast
>>> ds = xr.open_dataset(
>>>     FILE, engine="cfgrib", chunks="auto", backend_kwargs={"read_keys": ["projString"]}
>>> )
>>> ds.coords  # x/y have no coordinates!
Coordinates:
  * number      (number) int64 400B 1 2 3 4 5 6 7 8 ... 43 44 45 46 47 48 49 50
  * time        (time) datetime64[ns] 32B 2024-01-01 ... 2024-01-02T12:00:00
  * step        (step) timedelta64[ns] 56B 00:00:00 06:00:00 ... 1 days 12:00:00
  * surface     (surface) float64 8B 0.0
    latitude    (y, x) float64 8MB ...
    longitude   (y, x) float64 8MB ...
    valid_time  (time, step) datetime64[ns] 224B ...

>>> import pyproj
>>> pyproj.CRS.from_user_input(ds["dis06"].attrs["GRIB_projString"])
<Projected CRS: +proj=laea +lon_0=10.000000 +lat_0=52.000000 +a=63 ...>
Name: unknown
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: unknown
- method: Lambert Azimuthal Equal Area
Datum: unknown
- Ellipsoid: unknown
- Prime Meridian: Greenwich

Describe the solution you'd like

When opening a GRIB file, I want its x/y coordinates decoded and I want to retrieve the information on the data CRS.

>>> ds = xr.open_dataset(FILE, engine="cfgrib", chunks="auto")
>>> "x" in ds.coords
True
>>> "y" in ds.coords
True
>>> import rioxarray
>>> ds.rio.crs
CRS.from_wkt('PROJCS[...]')

Describe alternatives you've considered

libgdal-grib is a package that can decode grib files with the rasterio backend of xarray. It is able to decode the x/y coordinates, but not the other dimensions, which are convoluted into the usual band variable of GeoTIFF-like data. However, the x/y coordinates is exactly what I am looking for. Also, libgdal-grib cannot be installed with Python-only package managers like pip.

>>> import rioxarray
>>> ds_rio = xr.open_dataset(FILE, engine="rasterio", chunks="auto")
>>> ds_rio.coords  # ONLY x/y are decoded
Coordinates:
  * band         (band) int64 11kB 1 2 3 4 5 6 ... 1395 1396 1397 1398 1399 1400
  * x            (x) float64 8kB -1.819e+06 -1.814e+06 ... 3.171e+06 3.176e+06
  * y            (y) float64 8kB 2.287e+06 2.282e+06 ... -2.453e+06 -2.458e+06
    spatial_ref  int64 8B ...
>>> ds_rio.rio.crs
CRS.from_wkt('PROJCS["unnamed",GEOGCS["Coordinate System imported from GRIB file",DATUM["unnamed",SPHEROID["GRS80",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_center",52],PARAMETER["longitude_of_center",10],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Metre",1],AXIS["Easting",EAST],AXIS["Northing",NORTH]]')
>>> import odc.geo.xr
>>> ds_rio.odc.geobox

Additional context

No response

Organisation

No response