Retain time coordinate when converting GRIB to xarray

[JamesVarndell]

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

When converting from a GRIBReader object to an xarray object with to_xarray, squeeze() is called on the resulting xarray Dataset. If the variable in the Dataset is a function of a dimension of length 1, then this dimension is removed as both a dimension and a coordinate, leaving the variable as a function of only the remaining dimensions. This is removing useful metadata - for example, if you have a single timestep of a gridded field, this is removing the time coordinate completely, removing useful metadata.

To see an example of this in action, try the following:

import earthkit.data

data = earthkit.data.from_source(
    'cds',
    'reanalysis-era5-single-levels',
    {
        'product_type': 'reanalysis',
        'variable': '2m_temperature',
        'year': '2024',
        'month': '08',
        'day': '09',
        'time': '18:00',
    },
)

ds = data.to_xarray()

The resulting xarray Dataset has no time coordinate:

<xarray.Dataset> Size: 8MB
Dimensions:    (latitude: 721, longitude: 1440)
Coordinates:
  * latitude   (latitude) float64 6kB 90.0 89.75 89.5 ... -89.5 -89.75 -90.0
  * longitude  (longitude) float64 12kB 0.0 0.25 0.5 0.75 ... 359.2 359.5 359.8
Data variables:
    2t         (latitude, longitude) float64 8MB ...

Describe the solution you'd like

When calling to_xarray() on a GRIBReader, if there is a dimension of length 1, it should be dropped as a dimension but not as a coordinate. For example, using the same example as above, the coordinates and dimensions of the resulting xarray Dataset should look something like this:

<xarray.Dataset> Size: 4MB
Dimensions:     (latitude: 721, longitude: 1440)
Coordinates:
    valid_time  datetime64[ns] 8B 2024-08-09T18:00:00
  * latitude    (latitude) float64 6kB 90.0 89.75 89.5 ... -89.5 -89.75 -90.0
  * longitude   (longitude) float64 12kB 0.0 0.25 0.5 0.75 ... 359.2 359.5 359.8
Data variables:
    2t         (latitude, longitude) float32 4MB dask.array<chunksize=(721, 1440), meta=np.ndarray>

Notice that the valid_time coordinate remains even though it was dropped as a dimension.

Describe alternatives you've considered

No response

Additional context

No response

Organisation

ECMWF