Unable to read additional grib2 attributes

Chrismarsh commented 3 years ago

Broadly related to #89.

I'm noticing that various attributes are missing when I read in an Environment Canada grib file. The backend_kwargs arg approach is non optimal as it requires knowing the attributes before reading the file, making data exploration difficult.

However, my specific problem is that setting backend_kwargs with the known attribute does not seem to work as expected. In the following example I expect forecast_time_units to be present.

Datafile: CMC_hrdps_west_ALBDO_SFC_0_ps2.5km_2021021600_P000-00.grib2 (I am linking to the folder as the daily forecasts are deleted daily).

ds = xr.open_dataset('download_grib/CMC_hrdps_west_ALBDO_SFC_0_ps2.5km_2021021306_P001-00.grib2',
engine='cfgrib',  backend_kwargs={'read_keys': ['forecast_time_units']})

ds.al
Out[25]:
<xarray.DataArray 'al' (y: 485, x: 685)>
[332225 values with dtype=float32]
Coordinates:
    time        datetime64[ns] ...
    step        timedelta64[ns] ...
    surface     int64 ...
    latitude    (y, x) float64 ...
    longitude   (y, x) float64 ...
    valid_time  datetime64[ns] ...
Dimensions without coordinates: y, x
Attributes:
    GRIB_paramId:                    260509
    GRIB_shortName:                  al
    GRIB_units:                      %
    GRIB_name:                       Albedo
    GRIB_cfVarName:                  al
    GRIB_dataType:                   af
    GRIB_missingValue:               9999
    GRIB_numberOfPoints:             332225
    GRIB_typeOfLevel:                surface
    GRIB_NV:                         0
    GRIB_stepUnits:                  1
    GRIB_stepType:                   instant
    GRIB_gridType:                   polar_stereographic
    GRIB_gridDefinitionDescription:  Polar stereographic
    long_name:                       Albedo
    units:                           %

If I use (the now EOL) pynio loader (note that pynio loads the data as ALBDO_P0_L1_GST0 whereas cgrid loads as al)

 ds = xr.open_dataset('download_grib/CMC_hrdps_west_ALBDO_SFC_0_ps2.5km_2021021306_P001-00.grib2',engine='pynio')

ds.ALBDO_P0_L1_GST0
Out[5]:
<xarray.DataArray 'ALBDO_P0_L1_GST0' (ygrid_0: 485, xgrid_0: 685)>
[332225 values with dtype=float32]
Coordinates:
    gridlat_0  (ygrid_0, xgrid_0) float32 ...
    gridlon_0  (ygrid_0, xgrid_0) float32 ...
Dimensions without coordinates: ygrid_0, xgrid_0
Attributes:
    center:                                         Canadian Meteorological S...
    production_status:                              Operational products
    long_name:                                      Albedo
    units:                                          %
    grid_type:                                      Polar stereographic can b...
    parameter_discipline_and_category:              Meteorological products, ...
    parameter_template_discipline_category_number:  [ 0  0 19  1]
    level_type:                                     Ground or water surface
    forecast_time:                                  [60]
    forecast_time_units:                            minutes
    initial_time:                                   02/13/2021 (06:00)

I am unsure if I am doing something wrong or if there is a bug in the loader? I would also like to reiterate that being able to load all the attributes from the grib file without knowing a priori their names would be helpful!

shahramn commented 3 years ago

Doing a grib_dump on the product definition section of that file, I see:

% grib_dump -O -p section_4 CMC_hrdps_west_ALBDO_SFC_0_ps2.5km_2021021600_P000-00.grib2
======================   SECTION_4 ( length=34, padding=0 )    ======================
1-4       section4Length = 34
5         numberOfSection = 4
6-7       NV = 0
8-9       productDefinitionTemplateNumber = 0 [Analysis or forecast at a horizontal level or in a horizontal layer at a point in time (grib2/tables/4/4.0.table) ]
10        parameterCategory = 19 [Physical atmospheric properties (grib2/tables/4/4.1.0.table) ]
11        parameterNumber = 1 [Albedo  (%)  (grib2/tables/4/4.2.0.19.table) ]
12        typeOfGeneratingProcess = 2 [Forecast (grib2/tables/4/4.3.table) ]
13        backgroundProcess = 50
14        generatingProcessIdentifier = 50
15-16     hoursAfterDataCutoff = 0
17        minutesAfterDataCutoff = 0
18        indicatorOfUnitOfTimeRange = 0 [Minute (grib2/tables/4/4.4.table) ]
19-22     forecastTime = 0
23        typeOfFirstFixedSurface = 1 [Ground or water surface  (grib2/tables/4/4.5.table) ]
24        scaleFactorOfFirstFixedSurface = MISSING
25-28     scaledValueOfFirstFixedSurface = MISSING
29        typeOfSecondFixedSurface = 255 [Missing (grib2/tables/4/4.5.table) ]
30        scaleFactorOfSecondFixedSurface = MISSING
31-34     scaledValueOfSecondFixedSurface = MISSING

So the key that tells you the forecast time units is "indicatorOfUnitOfTimeRange" (0 means minute)

Chrismarsh commented 3 years ago

Thanks, that helps!

Is it possible to load all these attributes into the xarray .attrs?

ecmwf / cfgrib

Unable to read additional grib2 attributes #201