Message arrays using non-supported projection are returned as 1-d array

[corentincarton]

Hi,

I'm working on a python tool that deals with grib and netcdf inputs. Those inputs are defined on a regular grid that lives in the Lambert Azimuthal Equal Area projection. We want to use xarray with both the netcdf and cfgrib engines. We are well aware that the LAEA projection is not supported by eccodes (and therefore cfgrib), but we can reconstruct the coordinates by hand when we import our grib files in our scripts. I just have an issue regarding the dimension of the field that cfgrib returns when the projection is not supported. At the moment, xarray with cfgrib will return the following:

ECCODES ERROR   :  unable to get radius as double (Key/value not found)
ECCODES ERROR   :  Lambert Azimuthal Equal Area only supported for spherical earth.
ECCODES ERROR   :  unable to create iterator
ecCodes provides no latitudes/longitudes for gridType='lambert_azimuthal_equal_area'
>>> print(ds)
<xarray.Dataset>
Dimensions:     (step: 17, values: 4180000)
Coordinates:
    time        datetime64[ns] ...
  * step        (step) timedelta64[ns] 08:00:00 09:00:00 ... 1 days 00:00:00
    surface     int64 ...
    valid_time  (step) datetime64[ns] ...
Dimensions without coordinates: values
Data variables:
    tp          (step, values) float32 ...
Attributes:
    GRIB_edition:            2
    GRIB_centre:             ecmf
    GRIB_centreDescription:  European Centre for Medium-Range Weather Forecasts
    GRIB_subCentre:          0
    Conventions:             CF-1.7
    institution:             European Centre for Medium-Range Weather Forecasts
    history:                 2020-01-27T15:11:36 GRIB to CDM+CF via cfgrib-0....

We can see here that cfgrib returns an xarray dataset with dimension: (step: 17, values: 4180000) But, as we deal here with a regular grid, I would rather expect something like: (step: 17, y: 2200 x: 2000) With the current approach, the user has to figure out what is the leading dimension of the coordinates. Also, it is not straightforward to resize an array in xarray. This means it requires an extra processing step to redefine the dimensions, which can break the dask parallelisation (for instance, calling ds['tp'].value.resize(2200, 1900) will trigger the dask.compute() function).

Would it be possible to have cfgrib returning the correct dimensions, even if we can't reconstruct the coordinates values with eccodes? For instance, using the grid type gridDefinitionTemplateNumber together with the keys numberOfPointsAlongXAxis and numberOfPointsAlongYAxis.

I'm using the following version: eccodes 2.14.0 cfgrib 0.9.7.1

Best, Corentin

---

Computer Scientist Development Section ECMWF

[alexamici]

@corentincarton this is an interesting enhancement and I agree that it looks feasible. Can you share the URL to a GRIB file with these characteristics?

BTW @shahramn do you know if updating to a newer version of ecCodes solves the problem?

[corentincarton]

@alexamici here is a file you can use for your tests: https://drive.google.com/open?id=1Ud-j-HLJH_GdEQ2efGjo_ITnBtPL3Ufi

This projection will not be supported by eccodes any time soon.

[matteodefelice]

Any news on this? (Sorry for the cross-posting, see https://github.com/ecmwf/cfgrib/issues/28 )

[alexamici]

@corentincarton and @matteodefelice I've added lambert_azimuthal_equal_area to the semi-supported gridTypes. No lat lon coordintaes, but at least you get the 2-D array.

Will get into the upcoming 0.9.8.3 release

[shahramn]

The Geoiterator (Point Iterator) for Lambert azimuthal equal area projection on an ellipsoid has been implemented in ecCodes and will be released in the next version (v2.24.0, Winter 2021)